使用featureCount进行alignment-based的定量

featureCount是subread套件的一个模块，最大的优点就是速度非常快，使用全部overlap的reads计数，灵活考虑多比对的reads的计数。

subread 软件官网：http://subread.sourceforge.net/

官方文档：http://bioinf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf

安装subread

使用bioconda安装subread

conda install -y subread

从官网下载安装

wget -c https://sourceforge.net/projects/subread/files/latest/download -P ~/software/
cd ~/software/
# 解压
mv download subread-1.6.4-source.tar.gz
tar zxvf subread-1.6.4-source.tar.gz
# 安装
cd subread-1.6.4-source/
make -f Makefile.Linux
# 添加系统环境变量
echo 'PATH=$PATH:~/software/TrimGalore-0.6.2/bin' >> ~/.bashrc
source ~/.bashrc
# 测试

featureCounts的用法

Usage: featureCounts [options] -a <annotation_file> -o <output_file> input_file1 [input_file2] ...

featureCounts的选项和参数：

# 必须参数:
  -a <string>         输入GTF/GFF格式的。-F选项提供更多的格式信息。软件安装目录'annotation'中内建的SAF格式的注释文件也可以使用。Gzip压缩的文件也支持
  -o <string>         read计数结果的输出文件名。额外还有输出结果统计信息，保存在 ('<string>.summary')文件中。所有输出文件都以制表符分隔

  input_file1 [input_file2] ...   多个SAM 或 BAM格式输入文件。按照名称或位置排序。如果没有提供输入文件，则应输入<stdin>。 位置排序的双末端read将自动按read名称顺序读取

# 可选参数:
## 注释
  -F <string>         指定输入的注释文件的格式。
  										可接受的格式包括'GTF'/'GFF'/'SAF'。
  										默认是'GTF'
  -t <string>         指定GTF注释文件中的feature类型。默认：'exon'
  										用来计数的feature根据提供的参数从注释文件中提取
  -g <string>         指定GTF注释文件中属性的类型。默认：'gene_id'
  										从注释文件中提取Meta-features信息用于read计数
  --extraAttributes   从提供的GTF注释文件文件中提取额外的属性，将他们包含在计数结果中。 这些属性不会用于对特征分组，多个属性值使用逗号分隔

  -A <string>         包含注释文件中染色体名和read中染色体名的对应关系的文本文件，该文件是以逗号为分隔符的两列的表格。第一列包含注释文件的染色体名，第二列是reads文件中的染色体名，根据实际情况修改，不要表头
  
## summarization的水平
  -f                  在feature水平计数 (如对外显子计数)

## read和features的overlap
  -O                  统计所有重叠的的meta-features的reads 
  										(或-f参数指定的feature）
  --minOverlap <int>  最小重叠（overlapping）碱基数，默认: 1。 双末端read的重叠碱基数是两个read重叠碱基的总数。如果该值为负，允许read和feature之间存在该值绝对值长度的gap
  --fracOverlap <float> read最小的重叠分值。该值在[0,1]区间内，默认: 0
  											对于双末端序列重叠碱基数为一对read的重叠碱基数之和。
  											该选项和 '--minOverlap' 选项都必须满足read分布
  --fracOverlapFeature <float> 根据feature统计的read重叠碱基的最小分值, 
  														 该值在[0,1]之间, 默认：0
  --largestOverlap    计算meta-feature/feature的read重叠碱基的最大值。
  --nonOverlap <int>  根据feature计算的非read（read pair）重叠的最大值。
  										默认无限制
  --nonOverlapFeature <int> 用于read评估的feature中非overlap的重叠碱基最大数
  													默认无限制
  --readExtension5 <int> Read向 5'端上游扩展 <int> 个碱基
  --readExtension3 <int> Read向 3'端上游扩展 <int> 个碱基
  --read2pos <5:3>    	 减少read 5'端或 3'端的大多数碱基。
  											 通过read减少的单碱基进行计数

## 多mapping read
  -M                  对所有多比对read进行计数。通过BAM/SAM格式文件中 'NH' 标签来检测是否是多比对read

## 计数分数

  --fraction          feature统计的计数分数。该选项必须和 '-M' 或 '-O' 或both一起使用。
  和 '-M' 一起使用时, 多比对（通过'NH'标签来确定）的报告中会计算记一个计数分值1/x, 来代替1 (one), x是比对报告中相同read的总数
  和 '-O' 一起使用时, 每个重叠的feature都有一个计数分值1/y, y是read重叠的features总数
  和'-M' 和 '-O' 一起使用时, 每个比对都会得到计数分值1/(x*y)

## Read过滤
  -Q <int>            比对上的read满足最小质量值才被计数。
  										对于双末端read，至少有一个read满足标准。默认：0
  --splitOnly         只对split比对计数(如，比对结果中CIGAR字符串中包含'N')。												RNA-seq数据的exon-spanning read
  										就是split比对的一个实例。
  --nonSplitOnly      对非split比对计数（比对结果中CIGAR字符串中不包含'N'）。											其他比对结果会被忽略。
  --primary           只统计主要比对结果。通过BAM/SAM格式中FLAG列中包含 Ox100 二进制标识符来确定
  --ignoreDup         忽略duplicate reads计数。Duplicate reads通过BAM/SAM格式中FLAG列中包含 Ox400 二进制标识符来确定。如果双末端read的一个是duplicate read，则该read-pair被忽略

## Strandness
  -s <int or string>  执行链特异性read计数。单个整数（应用于所有输入）或以逗号分隔的字符串（应用于对应的输入文件），数量要与输入文件一致。可能的取值包括: 0 (unstranded), 1 (stranded), 2 (reversely stranded)。默认：0 (针对所有输入文件进行非链特异read计数)

## Exon-exon junctions
  -J                  支持exon-exon junction的read计数。Junctions是从输入中的 exon-spanning eads产生(CIGAR字符串的 'N' )。 
  										储存计数结果的文件名： '<output_file>.jcounts'
  -G <string>         提供生成SAM/BAM格式文件的FASTA格式的参考序列。
  										该选项和'-J'结合提高junctions的read计数

## 双末端read参数
  -p                  设置此参数后，使用fragments (或 templates)代替
  										reads计数。只对paired-end reads有效
  -B                  只对都比对上的read pairs计数
  -P                  检查paired-end距离使用 -d 和 -D 来设置阈值
  -d <int>            最小fragment/template长度，默认：50
  -D <int>            最大fragment/template长度，默认：600
  -C                  对于比对到不同染色体或相同染色体的不同链的read对，
  										不进行计数
  --donotsort         不对输入BAM/SAM文件的reads排序。请注意，来自同一对的read必须在输入中彼此相邻

# CPU线程数
  -T <int>            使用的线程数。默认1

# Read分组
  --byReadGroup       通过read group计数。 "RG" 标签必须要存在输入的BAM/SAM 文件中。
                      
# 长read
  -L                  对Nanopore 和 PacBio 的长reads计数。长read计数只能使用单核，只能对reads (不能是read-pairs)计数。此参数对长read计数的CIGAR字符串中'M'的数量无限制

# 每个read的Assignment结果
  -R <format>         输出每个read或read-pair的详细分析结果。 结果保存到以下格式之一的文件中：CORE，SAM和BAM。
  --Rpath <string>    指定目录以保存详细的分析结果。 如果未指定，则使用保存计数结果的目录。

# 其他参数
  --tmpDir <string>   保存运行过程中结果的文件夹（运行完成后删除）。默认保存在 '-o' 选项的参数设置的路径中。
  --maxMOp <int>      CIGAR字符串中'M'操作符的最大值。默认：10。CIGAR字符串中，'X' 和 '=' 被当作 'M'，临近的 'M' 操作被合并
  --verbose           输出debugging的信息, 未匹配的chromosome/contig name
  -v                  输出程序版本

featureCounts的使用示例

# 每次使用生信软件分析前，注意要激活创建的小环境rna
conda activate rna
# 将5.mapping目录下的bam文件链接到6.featureCounts
ln -s ~/5.mapping/*.bam ~/6.featureCounts/

### 运行featureCounts计数
featureCounts -T 5 -p -t exon -g gene_id -a /teach/database/gtf/gencode.v25.annotation.gtf.gz -o ~/6.featureCounts/all.id.txt *sort.bam
# 简化结果，去除特征列，只保留基因计数的列
cat all.id.txt | cut -f1,7- > counts.txt

featureCounts的结果解析

all.id.txt文件的表达矩阵：

# Program:featureCounts v1.6.4; Command:"featureCounts" "-T" "5" "-p" "-t" "exon" "-g" "gene_id" "-a" "/teach/database/gt
Geneid  Chr     Start   End     Strand  Length  SRR1039510.sort.bam	SRR1039511.sort.bam     SRR1039512.sort.bam
ENSG00000278267.1       chr1    17369   17436   -       68      0       0       0
ENSG00000268020.3       chr1    52473   53312   +       840     0       0       0
ENSG00000240361.1       chr1    62948   63887   +       940     0       0       0
ENSG00000186092.4       chr1    69091   70008   +       918     0       0       0

从左到右依次：

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