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.

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