转载的生信技能书的文章，仅用于自己学习记录！
文章原链接为：http://www.bio-info-trainee.com/2775.html ；可以直接去原链接学习。
史上最快的转录组流程-subread
安装软件
二进制版本软件，直接找到官网下载解压即可使用。
cd ~/biosoft
# http://bioinf.wehi.edu.au/featureCounts/
mkdir featureCounts && cd featureCounts
## 之前以为这个软件就是用来计算表达量的，所以把文件夹取名为 featureCounts
wget https://sourceforge.net/projects/subread/files/subread-1.5.3/subread-1.5.3-Linux-x86_64.tar.gz
tar zxvf subread-1.5.3-Linux-x86_64.tar.gz
建立索引
每个比对工具的算法不一样，所以每个工具都需要对参考基因组建立自己的索引。本身参考基因组占一篇空间就不小，索引之后更大!
需要自行从UCSC下载参考基因组，我放在了~/reference/genome/ 目录
buildindex=~/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/subread-buildindexcd/home/jianmingzeng/reference/index/subread/$buildindex-omm10 ~/reference/genome/mm10/mm10.fa$buildindex-ohg19 ~/reference/genome/hg19/hg19.fa$buildindex-ohg38 ~/reference/genome/hg38/hg38.fa
得到的索引文件如下：
749M Sep 15 17:37 hg19.00.b.array
4.9G Sep 15 17:37 hg19.00.b.tab
5.5K Sep 15 17:33 hg19.files
0 Sep 15 17:17 hg19.log
2.3K Sep 15 17:38 hg19.reads
766M Sep 15 18:01 hg38.00.b.array
5.0G Sep 15 18:01 hg38.00.b.tab
29K Sep 15 17:57 hg38.files
0 Sep 15 17:38 hg38.log
14K Sep 15 18:01 hg38.reads
652M Sep 15 17:17 mm10.00.b.array
4.4G Sep 15 17:17 mm10.00.b.tab
3.9K Sep 15 17:13 mm10.files
0 Sep 15 16:52 mm10.log
1.6K Sep 15 17:17 mm10.reads
批量比对
做好一个配置文件，就可以运行下面的脚本。
subjunc="/home/jianmingzeng/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/subjunc";
subjunc_mm10_index='/home/jianmingzeng/reference/index/subread/mm10';
cat $config |while read id
do
arr=($id)
fq1=${arr[1]}
fq2=${arr[2]}
sample=${arr[0]}
echo " start alignment for $sample" `date`
#$hisat -p 5 -x $mm10_index -1 $fq1 -2 $fq2 -S $sample.sam 2>$sample.hisat.log
#samtools sort -O bam -@ 5 -o $sample.bam $sample.sam
$subjunc -T 5 -i $subjunc_mm10_index -r $fq1 -R $fq2 -o ${sample}_subjunc.bam
echo " end alignment for $sample" `date`
done
配置文件就3列，第一列是样本名，第二列是该样本的fastq1，第二列是fastq2。多个样本的样本名不运行重复。
之前我以为hisat就很快了，换成了这个subjunc才知道没有最快，只有更快。
批量计算表达量
mm10_gtf='/home/jianmingzeng/reference/gtf/gencode/gencode.vM12.annotation.gtf';
featureCounts='/home/jianmingzeng/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/featureCounts';
$featureCounts -T 5 -p -t exon -g gene_id -a $mm10_gtf -o counts.txt *.bam
实在是没有想到这个软件居然会如此快，1M的reads耗时三五秒即可，甩之前的htseq-counts好几条街。
还有更多计算的模型和参数可以供挑选；http://bioinf.wehi.edu.au/featureCounts/
