2021-05-25史上最快的转录组流程-subread

转载的生信技能书的文章，仅用于自己学习记录！

文章原链接为：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/