参考:有参转录组实战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#报错说无权限的话运行这个,然后再重新运行一下上面步骤
差不多就到这里