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

史上最快的转录组流程-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-buildindex
cd /home/jianmingzeng/reference/index/subread/
$buildindex -o mm10  ~/reference/genome/mm10/mm10.fa
$buildindex -o hg19  ~/reference/genome/hg19/hg19.fa
$buildindex -o hg38  ~/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/

Comments are closed.