17-有参转录组实战3-计算readcount和TPM表达量

发布时间 2023-11-30 17:29:19作者: 啊辉的科研

 

#本教程部分文件,来源B站15天入门生物信息教程

“https://www.bilibili.com/video/BV1K44y1B7Dg/?spm_id_from=333.337.search-card.all.click&vd_source=19eea84b7c1e944fd3c6b3ccca066ade”

#1,先将毛果杨的GFF格式的注释文件转换为GTF格式:

conda install gffread

gffread Ptri_genome.gff -T -o Ptri_genome.gtf

#2,在B站教程的评论区里面得到压缩包rnaseq-apple-training.zip,在其中找到run-featurecounts.R文件,将其放入工作文件夹中,生成批量的R脚本命令:

awk '{print "Rscript run-featurecounts.R --bam "$3".sort.bam --gtf Ptri_genome.gtf --output "$3" &"}' sample.txt >command_Rscript.sh

#3,检查命令:

 

#4,进入R环境,请务必并提前安装好R包(argparser,Rsubread,limma,edgeR),运行R脚本:

source activate R

sh command_Rscript.sh

#5,约20min,得到了count文件,我们需要将其合并,这里需要用到压缩包里的abundance_estimates_to_matrix.pl文件,和support_scripts文件夹,都弄到服务器里:

ls *.count > genes.quant_files.txt

perl abundance_estimates_to_matrix.pl --est_method featureCounts --quant_files genes.quant_files.txt --out_prefix genes

#6,这里有个报错,说support_scripts文件夹里面的文件无权限,赋予权限:

chmod 777 support_scripts/run_TMM_scale_matrix.pl

#7,再重新运行下,就会出来结果了。gene.count.matrix是readcount矩阵文件,genes.TPM.not_cross_norm是TPM文件,后续的分析不建议使用FPKM。本教程的脚本文件来源于B站的15天入门生物信息视频,若有疑惑请以视频为准。若侵请联系删除。

#8,步骤3可以自己尝试使用for循环。