0.前期准备
先在工作目录下创建以下几个目录:
01.raw_data #用于存放原始数据
02.fastq #用于存放fastq格式数据
03.fastqc #用于存放QC结果
04.trinity_result #用于存放trinity结果
数据处理的顺序:01.raw_data > 02.fastq > 03.fastqc > 04.trinity_result
1.下载SRA
先cd
到01.raw_data
目录下,执行以下命令(二选一):
nohup axel -q ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR130/SRR1300434/SRR1300434.sra &
axel是一种多线程下载工具,比wget快;-q
参数能让axel进入静默模式(选用)。
nohup wget -c ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR130/SRR1300434/SRR1300434.sra &
用wget的时候记得加-c参数,支持断点续传,就会方便很多。
用脚本批量下载数据的代码参考附录
里的内容。
2.fastq-dump
fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files file.sra --split-3 -O|--outdir
这个是trinity给的命令,如果后面要用trinity进行无参组装的就需要加上--defline-seq '@$sn[_$rn]/$ri'
这串看不懂的“乱码”,就要用这个命令,否则trinity会报错说无法识别这个格式。
经过改良后的命令:
nohup fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-3 file.sra -o ../02.fastq &
默认输出在上一层目录的02.fastq
这个文件夹里。
-o
:重定向输出结果的位置。将结果输出到上一级目录下的02.fastq
文件夹内。
可选用的参数:
--gzip, --bzip2
: 以压缩文件的方式输出结果。有利于减少文件的占用空间。
3.质控
不管什么数据,质控一下还是很有必要的,多看fastqc结果,即是积累经验,也是及时发现数据可能出现的问题。首先cd
到../02.fastq
,之后执行:
fastqc -t 8 *.fq -o ../03.fastqc
-t
: 使用的线程数。现在用的服务器是20 CPU 40线程,以后这个值可以改大一点(range from 1~40,但是一般取10,如果文件比较多可以取20。)
-o
: 重定向输出结果的位置。放在03.fastqc
里面。
用filezilla
登陆服务器找到对应文件夹将生成的html文件下载到本地进行查看。如无问题则可进行下一步Trinity的操作,如有问题则进行质控部分余下的操作。
因为篇幅原因,将单独写一个关于如何查看fastqc结果的教程。敬请期待
需要去除接头的情况:
- cutadapt
- trimmomatic
需要对双端测序数据进行过滤低质量reads操作的情况:
- sickle
当你想改善你的数据的问题但是你又新入门啥软件都不会的时候(推荐使用):
- fastp
fastp -i name.fastq -o name.fastp.fastq #单端
fastp -i name_1.fastq -o name_1.fastp.fastq -I name_2.fastq -O name_2.fastp.fastq #双端
4.Trinity
cd
到../fastq
文件夹后,执行以下命令:
双端测序命令:
nohup Trinity --seqType fq\fa --left ? --right ? --CPU 20 --max_memory 50G -output ../trinity_result/?.trinity &
left和right不能随意设置,left必须是_1,right必须是 _2.
一般type是fq
将
?
的位置用具体需要操作的文件名代入即可。
单端测序命令:
nohup Trinity --seqType fq\fa --single ? --CPU 20 --max_memory 50G -output ../trinity_result/?.trinity &
ps:输出结果的名字里必须含有trinity这个关键词,否则无法正常运行。
5.TransDecoder
cd
到../trinity_result/filename.trinity
文件夹下后,找到一个叫做Trinity.fasta
的文件,执行:
#step 1: 提取最长的开放阅读框
TransDecoder.LongOrfs -t ./Trinity.fasta
#step 2: (可选),BlastP搜索和Pfam搜索
#BlastP搜索:蛋白库搜索, Swissprot (快) or Uniref90 (慢 but more comprehensive)
blastp -query transdecoder_dir/longest_orfs.pep -db uniprot_sprot.fasta -max_target_seqs 1 -outfmt 6 -evalue 1e-5 -num_threads 10 > blastp.outfmt6
#Pfam搜索:肽或蛋白域预测,需要安装hmmer3和Pfam数据库
hmmscan --cpu 8 --domtblout pfam.domtblout /path/to/Pfam-A.hmm transdecoder_dir/longest_orfs.pep
#step 3: 将Blast和Pfam搜索结果整合到编码区域选择这个步骤在尝试中发现比较费时间,建议用nohup命令进行操作。
nohup TransDecoder.Predict -t target_transcripts.fasta --retain_pfam_hits pfam.domtblout --retain_blastp_hits blastp.outfmt6 &
nohup TransDecoder.Predict -t trinity.fasta & #简化版
通过TransDecoder这一步将能获得cds、pep、gff3、bed文件。用filezilla将以下五个文件拿到本地:
Trinity.fasta
、Trinity.fasta.transdecoder.cds
、Trinity.fasta.transdecoder.pep
、Trinity.fasta.transdecoder.gff
、-
Trinity.fasta.transdecoder.bed
将这几个文件以可分辨的名字进行命名(如改成SRR1300215.fasta等)。
6.出个统计报告
还是在本目录下,对Trnity.fasta文件进行操作:
TrinityStats.pl trinity_out_dir/Trinity.fasta >> stat.txt
会生成一个统计报告,并将内容重定向到stat.txt这个文件中。这个文件也记得保存下来。
结果展示。
以SRR1300215为例:
################################
## Counts of transcripts, etc.
################################
Total trinity 'genes': 1524
Total trinity transcripts: 1793
Percent GC: 54.43
########################################
Stats based on ALL transcript contigs:
########################################
Contig N10: 1008
Contig N20: 757
Contig N30: 597
Contig N40: 487
Contig N50: 397
Median contig length: 291
Average contig: 379.04
Total assembled bases: 679625
#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################
Contig N10: 965
Contig N20: 715
Contig N30: 559
Contig N40: 461
Contig N50: 382
Median contig length: 287
Average contig: 369.46
Total assembled bases: 563055</pre>
7.附录
批量下载sra文件的脚本模板(mac版本,仅供参考)
# 在终端运行
#用axel下载会快一些
brew install axel
for i in {56..62}
do
# 也可以用wget 下载
axel ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP075/SRP075747/SRR35899$i/SRR35899$i.sra
done
批量处理fastq-dump的脚本模板:(从人家的教程里摘下来的,仅供参考)
#!/usr/bin/env bash
for i in {56..62}
do
fastq-dump --gzip --split-3 -O /Users/chengkai/Desktop/zhuanlu_files -A SRR35899${i}.sra
done
8.参考来源
文章:
-
转录组入门二(mac版):读文献下数据
Author:thinkando
-
Fastq-dump: 一个神奇的软件
Author:hoptop
-
转录组入门三(mac 版):了解fastqc测序数据
Author:thinkando
-
使用Trimmomatic过滤低质量序列
Author:To_2019_1_4
-
TRANSDECODER寻找和预测ORF
Author:冰冻三丈
微信公众号:
生信菜鸟团
生信者言
9.版权信息
Author:陈俊浩 Hanschen
QQ:1020607557