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.80.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}"
)