RNA-seq上游分析的最后一站：featureCounts

官网：http://bioinf.wehi.edu.au/featureCounts/。featureCounts与subread、subjunc软件一样同属subread包。

featureCounts 需要两个输入文件：

1.比对产生的BAM/ SAM文件（教程中用bam文件，因为bam文件占用空间小）
2.区间注释文件（GTF格式， SAF格式）

subread包

featureCounts常用参数

-a # 输入GTF/GFF基因组注释文件(annotation file)
-F # 指定区间注释文件的格式，默认是GTF
-o # 输出文件:可输出raw counts的txt文本及raw counts的summary文本
-p # 针对paired-end数据
-g # 指定统计的meta-feature名称，默认是gene_id
-t # 指定统计的feature名称，默认是exon
-T # 指定线程数

featureCounts的结果解析

Geneid：基因的ensemble基因号；
Chr：多个feature所在的染色体编号;
Start：多个feature起始位点,与前面一一对应
End：多个feature终止位点,与前面一一对应
Strand：正负链
Length：基因长度
sampleID：一列代表一个样本，数值表示比对到该基因上的read数目

MultiQC分析featureCounts结果

multiqc all.id.txt.summary  #对featureCounts的结果进行整合，html文件可视化

featurecounts的运行

featureCounts #打开featureCounts的帮助文档
#每次使用生信软件分析前，都必须要激活创建的小环境rna
conda activate rna

workdir=$HOME/project/airway   #建立目录变量
# 进入06.counts文件夹
cd ${workdir}/06.counts
# 将05.mapping文件夹下面的所有sort.bam文件链接到06.counts文件夹下。
ln -s ${workdir}/05.mapping/*.sort.bam ./

# 设置featureCounts所需要的gtf文件位置（注释文件）
gtf=/teach/database/gtf/gencode.v25.annotation.gtf.gz

# 运行featureCounts对所有样本进行基因水平定量
featureCounts -T 1 -p \
-t exon -g gene_id \
-a $gtf -o all.id.txt *.sort.bam
#解释：-T 1：线程数为1；-p：表示数据为paired-end，双末端测序数据；-t exon表示 feature名称为exon；-g gene_id表示meta-feature名称为gene_id(ensembl名称)；-a $gtf 表示输入的GTF基因组注释文件；-o all.id.txt 设置输出文件类型；*.sort.bam是被分析的bam文件。featureCounts支持通配符*

# 使用multiqc对featureCounts统计结果进行可视化。
multiqc all.id.txt.summary

# 只保留all.id.txt文件的【基因名】和【样本counts】
cat all.id.txt | cut -f1,7- > counts.txt