SUBREAD FEATURECOUNTS¶
FeatureCounts assign mapped reads or fragments (paired-end data) to genomic features such as genes, exons and promoters. For more information please see featureCounts tutorial, documentation of subread and commandline help.
URL:
Example¶
This wrapper can be used in the following way:
rule feature_counts:
input:
sam="{sample}.bam", # list of sam or bam files
annotation="annotation.gtf",
# optional input
# chr_names="", # implicitly sets the -A flag
# fasta="genome.fasta" # implicitly sets the -G flag
output:
multiext("results/{sample}",
".featureCounts",
".featureCounts.summary",
".featureCounts.jcounts")
threads:
2
params:
tmp_dir="", # implicitly sets the --tmpDir flag
r_path="", # implicitly sets the --Rpath flag
extra="-O --fracOverlap 0.2"
log:
"logs/{sample}.log"
wrapper:
"0.85.0/bio/subread/featurecounts"
Note that input, output and log file paths can be chosen freely.
When running with
snakemake --use-conda
the software dependencies will be automatically deployed into an isolated environment before execution.
Software dependencies¶
subread=2.0
Input/Output¶
Input:
- a list of .sam or .bam files
- GTF, GFF or SAF annotation file
- optional a tab separating file that determines the sorting order and contains the chromosome names in the first column
- optional a fasta index file
Output:
- .featureCounts file including read counts (tab separated)
- .featureCounts.summary file including summary statistics (tab separated)
- .featureCounts.jcounts file including count number of reads supporting each exon-exon junction (tab separated)
Authors¶
Code¶
__author__ = "Antonie Vietor"
__copyright__ = "Copyright 2020, Antonie Vietor"
__email__ = "antonie.v@gmx.de"
__license__ = "MIT"
from snakemake.shell import shell
log = snakemake.log_fmt_shell(stdout=True, stderr=True)
extra = snakemake.params.get("extra", "")
# optional input files and directories
fasta = snakemake.input.get("fasta", "")
chr_names = snakemake.input.get("chr_names", "")
tmp_dir = snakemake.params.get("tmp_dir", "")
r_path = snakemake.params.get("r_path", "")
if fasta:
extra += " -G {}".format(fasta)
if chr_names:
extra += " -A {}".format(chr_names)
if tmp_dir:
extra += " --tmpDir {}".format(tmp_dir)
if r_path:
extra += " --Rpath {}".format(r_path)
shell(
"(featureCounts"
" {extra}"
" -T {snakemake.threads}"
" -J"
" -a {snakemake.input.annotation}"
" -o {snakemake.output[0]}"
" {snakemake.input.sam})"
" {log}"
)