史上最快的转录组定量流程

前面我比较过各种alignment的工具，5款流行比对工具大比拼 提到了subjunc，它其实是subread套件的一个小功能，而很久以前就有朋友建议我用featureCounts替换之前教程的htseq-counts，直到现在才有空写成教程分享给大家。

安装软件

二进制版本软件，直接找到官网下载解压即可使用。

  1. cd ~/biosoft

  2. # http://bioinf.wehi.edu.au/featureCounts/

  3. mkdir featureCounts &&  cd featureCounts

  4. ## 之前以为这个软件就是用来计算表达量的，所以把文件夹取名为 featureCounts

  5. wget https://sourceforge.net/projects/subread/files/subread-1.5.3/subread-1.5.3-Linux-x86_64.tar.gz

  6. tar zxvf subread-1.5.3-Linux-x86_64.tar.gz

建立索引

每个比对工具的算法不一样，所以每个工具都需要对参考基因组建立自己的索引。本身参考基因组占一篇空间就不小，索引之后更大!

需要自行从UCSC下载参考基因组，我放在了 ~/reference/genome/ 目录

  1. buildindex=~/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/subread-buildindex

  2. cd /home/jianmingzeng/reference/index/subread/

  3. $buildindex -o mm10  ~/reference/genome/mm10/mm10.fa

  4. $buildindex -o hg19  ~/reference/genome/hg19/hg19.fa

  5. $buildindex -o hg38  ~/reference/genome/hg38/hg38.fa

得到的索引文件如下：

  1. 749M Sep 15 17:37 hg19.00.b.array

  2. 4.9G Sep 15 17:37 hg19.00.b.tab

  3. 5.5K Sep 15 17:33 hg19.files

  4.   0 Sep 15 17:17 hg19.log

  5. 2.3K Sep 15 17:38 hg19.reads

  6. 766M Sep 15 18:01 hg38.00.b.array

  7. 5.0G Sep 15 18:01 hg38.00.b.tab

  8. 29K Sep 15 17:57 hg38.files

  9.   0 Sep 15 17:38 hg38.log

  10. 14K Sep 15 18:01 hg38.reads

  11. 652M Sep 15 17:17 mm10.00.b.array

  12. 4.4G Sep 15 17:17 mm10.00.b.tab

  13. 3.9K Sep 15 17:13 mm10.files

  14.   0 Sep 15 16:52 mm10.log

  15. 1.6K Sep 15 17:17 mm10.reads

批量比对

做好一个配置文件，就可以运行下面的脚本。

  1. subjunc="/home/jianmingzeng/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/subjunc";

  2. subjunc_mm10_index='/home/jianmingzeng/reference/index/subread/mm10';

  3. cat $config |while read id

  4. do

  5.    arr=($id)

  6.    fq1=${arr[1]}

  7.    fq2=${arr[2]}

  8.    sample=${arr[0]}

  9.    echo "  start alignment for $sample" `date`

  10.    #$hisat -p 5 -x $mm10_index -1 $fq1 -2 $fq2 -S $sample.sam 2>$sample.hisat.log

  11.    #samtools sort -O bam -@ 5  -o $sample.bam $sample.sam

  12.    $subjunc -T 5  -i $subjunc_mm10_index -r $fq1  -R $fq2 -o ${sample}_subjunc.bam

  13.    echo "  end alignment for $sample" `date`

  14. done

配置文件就3列，第一列是样本名，第二列是该样本的fastq1，第二列是fastq2。多个样本的样本名不运行重复。

之前我以为hisat就很快了，换成了这个subjunc才知道没有最快，只有更快。

批量计算表达量

  1. mm10_gtf='/home/jianmingzeng/reference/gtf/gencode/gencode.vM12.annotation.gtf';

  2. featureCounts='/home/jianmingzeng/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/featureCounts';

  3. $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/

赞 (0)

相关推荐