安装软件
二进制版本软件,直接找到官网下载解压即可使用。
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/