有参转录组定量

参考:有参转录组实战1-批量质控-fastp - 简书 (jianshu.com)
有参转录组实战2-将批量转录组比对到基因组上 - 简书 (jianshu.com)
有参转录组实战3-计算readcount和TPM表达量 - 简书 (jianshu.com)

首先下载转录组,新建一个SRR_Acc_List.txt,大概长这样,每行是SRR号码


image.png

然后在linux中新建一个文件夹,把SRR_Acc_List.txt放进去,然后用sratoolkit下载sra文件

/路径/sratoolkit.3.0.0-centos_linux64/bin/prefetch --option-file  SRR_Acc_List.txt

要下很久,接下来把sra文件转换成fastq.gz

find * -name '*.sra' -exec mv {} /路径/  \;#把所有sra文件挪到一起。\;是exec结束
awk '{print "parallel-fastq-dump -t 12 --outdir /路径/ --split-3 --gzip -s /路径/"$1".sra -T /路径/tmp/ &"}' SRR_Acc_List.txt >command_para.sh
tr -d '\r' < command_para.sh > command_Rscript_fixed.sh
mv command_Rscript_fixed.sh command_para.sh#生成批量命令,转换成fastq.gz,并用tr -d解决换行不正确的问题,中间那个sh文件无所谓名字,就是一个中转文件
sh command_para.sh#批量转换,parallel超快的呢

接下来批量质控fastp,先重命名一个sample.txt,大概长这样,每行后面有个空格,最后一行是空行,


1720868519098.png
awk '{print "fastp -i "$1".fastq.gz -I "$2".fastq.gz -o "$1".clean.fq.gz -O "$2".clean.fq.gz -h "$3".html &"}' sample.txt>command_fastp.sh#生成批量质控的文件,同样用刚才的tr -d解决一下换行的问题
sh command_fastp.sh

然后看一下html文件,Q30>85就行好像,留下clean.fastq.gz,其余的删掉
接下来把转录组比对到基因组上

hisat2-build genome.fa genome#构建索引
awk '{print "nohup hisat2 -x genome -1 "$1".clean.fq.gz -2 "$2".clean.fq.gz -S "$3".sam > "$3".temp.txt 2>&1 &"}' sample.txt >command_hisat2.sh#批量命令,同样tr -d一下
sh command_hisat2.sh

要好久,好了以后再用samtools把sam文件变成bam文件

awk '{print "samtools view -bS "$3".sam > "$3".bam &"}' sample.txt > command_samtobam.sh#批量,同样tr -d
sh command_samtobam.sh

弄完以后再用samtools产生sort.bam文件,排序方便后续分析好像

awk '{print "samtools sort -o "$3".sort.bam "$3".bam &"}' sample.txt > command_sort.sh
sh command_sort.sh#同样的批量,同样的tr -d

结束了,留下sort.bam文件,其他都删掉
接下来是定量部分,从“https://www.bilibili.com/video/BV1K44y1B7Dg/?spm_id_from=333.337.search-card.all.click&vd_source=19eea84b7c1e944fd3c6b3ccca066ade”评论找到rnaseq-apple-training.zip,找到run-featurecounts.R和support_scripts文件夹,放到目录里

awk '{print "Rscript run-featurecounts.R --bam "$3".sort.bam --gtf genomic.gtf --output "$3" &"}' sample.txt >command_Rscript.sh#先生成批量文件,然后进入R幻境
source activate R
sh command_Rscript.sh#运行定量命令
ls *.count > genes.quant_files.txt#把所有的count文件放到一起
perl abundance_estimates_to_matrix.pl --est_method featureCounts --quant_files genes.quant_files.txt --out_prefix genes#定量
chmod 777 support_scripts/run_TMM_scale_matrix.pl#报错说无权限的话运行这个,然后再重新运行一下上面步骤

差不多就到这里

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 219,366评论 6 508
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,521评论 3 395
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 165,689评论 0 356
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,925评论 1 295
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,942评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,727评论 1 305
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,447评论 3 420
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,349评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,820评论 1 317
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,990评论 3 337
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,127评论 1 351
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,812评论 5 346
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,471评论 3 331
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,017评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,142评论 1 272
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,388评论 3 373
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 45,066评论 2 355

推荐阅读更多精彩内容