生物信息学/featureCounts
< 生物信息学
使用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数目