RSEM软件
##step1 准备参考文件:共三个参数:
rsem-prepare-reference --gtf gtf文件 \
参考基因组文件 \
位置/命名
##step2 定量计算
rsem-calculate-expression \
--paired-end \ (双端测序)
--no-bam-output \ (只定量计算,不输出bam文件结果)
--alignments \ (输入文件本身就是STAR比对输出的结果,这个参数代表输入文件是比对好的bam,不是原始fastq文件)
-p 15 \ (线程数)
-q 03align_out/sample1_1Aligned.toTranscriptome.out.bam \ (转换为转录本后的bam文件)一定要是转换为转录本后的bam文件
arab_RSEM/arab_rsem \ (索引)
04rsem_out/sample1_1rsem (输出目录/输出文件)
会生成基于基因的定量结果_rsem.gene.results和基于转录本的定量结果_rsem.isoform.results,isoforms.results中包含了转录本ID,基因ID,转录本长度,有效长度,expected_count,TPM,FPKM和IsoPct(该转录本表达量占基因总表达量的百分比)。genes.results中的内容与之类似,只是少了IsoPct
##step3:整合结果,构建表达矩阵
#将n个rsem.genes.results整合到1个矩阵里面
rsem-generate-data-matrix *rsem.genes.results > output_matrix
##step4:删除在所有样本里表达量为0的基因
awk 'BEGIN{printf "geneid\ta1\ta2\tb1\tb2\n"}{if($2+$3+$4+$5>0)print $0}' output_matrix > edgeR_input.txt
#查看删除前后文件的行数
wc -l output_matrix
wc -l edgeR_input.txt
featurecounts软件
sudo apt install subread#安装 featureCounts是subread软件包中的一种工具
featureCounts --help #可以使用
featureCounts \
-p \ (双端测序)
-a \ (gtf文件)
-o \ (输出文件)
-T \ (线程)
-t \ 想要统计的feature 的名称, 取值范围是gtf 文件中的第3列的值,默认是exon
-g \ 想要统计的meta-feature 的名称,取值范围参考gtf 第9列注释信息,
gtf 的第9列为 key=value 的格式, -g 参数可能的取值就是所有的key, 默认值是gene_id
_Aligned.sortedByCoord.out.bam \ bam文件
2020-08-15总结:STAR+featureCounts和STAR+RSEM软件时平行的,前者是基于基因组比对,后者是基于转录本比对,目前用STAR+featureCounts即可。
2020.08.30
在GSE46224文件夹里,对STAR输出的bam文件用featurecounts进行定量
对SRR830973-SRR830988进行定量,写脚本:
for i in {830973,830974,830975,830976,830977,830978,830979,830980,830981,830982,830983,830984,830985,830987,830988}
do
featureCounts -p -a /home/sxw/HF/Homorefer/humangtf/genecode38.gtf -o SRR${i}.counts.txt \
-T 5 -t exon -g gene_id SRR${i}_Aligned.sortedByCoord.out.bam
done
结果:
53.7%的比对成功率,可以接受。