因为我是直接trinity开始的,后面面才开始补加数据前处理,毕竟都是练手所以没关系。但实际上这部分应该放在前面。
下载一个基因组很小的细菌做试炼:Pelagibacter phage Greip EXVC021P
1. 下载解压
nohup wget www.XXXX &
nohup fastq-dump -gzip -split-3 -A SRR11559267 &
gunzip SRR11559267_1.fastq.gz
2. 原始数据质控图#####fastqc######
安装fastqc
conda create -n fastqc
conda activate fastqc
conda install -c bioconda fastqc
fastqc --help
运行fastqc
fastqc -t 4 -o ./ SRR11559267_1.fastq SRR11559267_2.fastq
得到文件SRR11559267_1_fastqc.html SRR11559267_2_fastqc.html
浏览器查看质控文件,得出的结果不是很好
3. 对reads进行过滤
NGSQCToolkit官网打不开了,而且用的人不是很多。
安装NGSQCToolkit前先安装libgd和GD[gd库提供了一系列用来处理图片的API,使用GD库可以处理图片,或者生成图片。在网站上GD库通常用来生成缩略图或者用来对图片加水印或者对网站数据生成报表。也就是有了gd库,我们用php对图片的处理将会得心应手。]
#conda install libgd#
A. FASTX-Toolkit
在使用过程前,我们需要简单的判断下这个测序格式是Phred+33还是Phred+64。有=号的一般都是Phred+33,实际上最近几年测序结果一般都是Phred+33,网上下载的早期数据可能是Phred+64。有人也这么判断:
grep 2 rosalind_filt_1_dataset.txt #有结果
grep X rosalind_filt_1_dataset.txt # 无结果
# 基本上断定这个是Phred33
B. FASTQ/A Clipper 去接头序列
这里 -v可以显示输入和输出功能, -l 18是去掉长度小于18nt的reads,要善用fastx_clipper -h,这样就能挑选自己想要的参数。-Q 33 在Fastx Toolkit的应用中都要加,这个在-h中不显示,我暂时能找到的解释是这么说的-Q is an undocumented parameter to indicate that quality values use ASCII 33 encoding。结果显示如下:
fastx_clipper -Q 33 -l 18 -a TGGAATTCTCGGGTGCCAAGG -v -i SRR11559267_1.fastq -o SRR11559267_1_clipped.fastq
关于接头序列:自己的序列自然好说,在网上下载的一些序列,比如NCBI是找不到接头序列的,这个时候可以借助Fastqc工具,在结果里面的Adapter Content里,找到打红色叉的adapter content ,可以通过这几个序列进行一个统计,作为接头adapter。如果没有太多接头残留,这一步可以省略。
C. fastq_quality_filter 去低质量的reads
fastq_quality_filter -Q 33 -v -q 30 -p 80 -i SRR11559267_1.fastq -o SRR11559267_1_qualified.fastq
关于 -q 和 -p 下列这个图片解释的很清楚,-q 30 -p 80 所过滤掉的reads介于-q 20 -p 90和-q 20 -p 100之间。
D. trimmomatic 过滤低质量碱基
如果使用Fastqc发现序列前后几个碱基质量不太好时,我们可以使用trimmomatic过滤掉按照一定的阈值对read前后进行过滤。例如
java -classpath trimmomatic-0.22.jar org.usadellab.trimmomatic.TrimmomaticSE -phred33 data/s1.fq data/tmp.fq TRAILING:30 MINLEN:50
就是过滤前后低于30的碱基,然后删掉不足50的read。
#单端
trimmomatic SE -phred33 SRR11559267_2.fastq out1.fq LEADING:22 TRAILING:22
#双端
trimmomatic PE -threads 5 -phred33 SRR11559267_1.fastq SRR11559267_2.fastq -baseout SRR11559267.fastq SLIDINGWINDOW:4:5 LEADING:5 TRAILING:5 MINLEN:25
在学习过程中,笔者发现有很多软件都可以做质控过滤,比如提到的 FASTX-Toolkit、trimmomatic,还有sickle、seqtk。
有帖子比较了trimmomatic、sickle、seqtk
结果显示:
如果需要同时去除测序接头序列,那么建议使用Trimmomatic;
如果只需过滤低质量碱基或者低质量read,可以选择Trimmomatic或sickle,有时sickle会更快一些,
如果不愿意read被过滤掉而且质量值体系是phred33,那么可以选择seqtk。
在笔者自己的试炼中,接头信息很少,然而FASTX-Toolkit好像不能解决PE去除低质量的reads后导致的左右两端测序文件不平衡问题(如果能平衡请批评指正),所以又学习了sickle 。使用说明详情请参照下一条博文:sickle转录组数据过滤·使用案例
PS:如何利用conda安装已经下载好的软件包:
1. wget下载.tar.gz2文件 然后移入miniconda/pkgs/文件夹
2. pkgs里面找到urls.txt文件 并手动添加下载地址
3. install
几个比较好的组装案例:
1. 组装细菌基因组-Trimmomatic
2. 转录组分析学习笔记(持续补充)
3. phylogenomic_dataset_construction