计算机的角度来说,生物的序列属于一种字符串,也是一种文本,因此生物信息分析属于文本处理范畴。文本存储为固定格式文件,生物信息的工作就是各种文本文件之间格式的转换,例如通过序列拼接将fastq转换为fasta,通过短序列比对将fastq与fasta合并为bam,通过变异检测将bam中突变位点提取出来转换为vcf。因此,我们可以通过总结每一种生物数据文件格式的处理方法来学习生物信息,这样当拿到固定格式的文件之后,就知道该如何来处理了。
fastq格式文件处理大全(一)fastq格式文件处理大全(二)
fastq格式文件处理大全(三)
去除接头adapter
接头adapter主要是指illumina测序时加入的P7接头与P5接头,理论上来说测序从测序引物开始,是测序不到接头的,但是由于部分打断片段果断或者由于adapter空载,会导致adpter自身连接,以上两种情况都会导致测序reads中包含adapter序列。adapter序列非基因组本身的序列,会干扰分析,因此需要去除掉。主要是illumina测序中包含adapter。去除adapter主要是采用两种方式,一种是直接给定adapter序列,与reads比对,比对上的就把整条reads去掉,另外一种是测序完成之后,给定一个adater list文件,其中包含了含有adapter序列的reads ID列表,给定一个阈值将这些reads去除掉。cutadapt可以根据给定的adapter序列进行过滤。
cutadapt -g AACMGGATTAGATACCCKG -a GGAAGGTGGGGATGACGT -o output1.fastq -p output2.fastq SRR8651554_1.fastq.gz SRR8651554_2.fastq.gz
去除低质量
低质量主要是指去除测序质量错误率较高的位点,一般以Q20作为标准,如果一个碱基质量值低于Q20,则认为一个低质量,如果一条序列中低质量碱基达到一定比例,例如达到40%以上,则过滤掉此条序列。是过滤数据主要处理指标。seqtk工具seq功能通过-q,-n可以将低质量碱基进行标记,例如替换为小写字母或者其他字符,但是不进行过滤,有专门的数据过滤工具。
#将小于Q20的替换为小写字母
seqtk seq -q 20 SRR8651554_1.fastq.gz |less -S
去除N碱基
如果测序仪无法准确判断出测序碱基的类型,则选择输出N,N碱基无法进行各种分析,因此需要去除掉序列中包含过多N碱基的数据。去除N碱基并不是讲N碱基切除,而且去除包含N碱基过多的整条数据,例如N碱基含量达到10%以上,则过滤掉数据,有些程序按照连续N碱基比率进行过滤。
去除Duplication
duplication是指一对完全一样的测序数据,是由于打断不随机或者PCR过程中导致的,duplication会干扰序列拼接,还会影响变异检查,因此去要去除掉。但是RNAseq和宏基因组测序由于本身序列短,并且丰度不同,因此不能去除duplication。去除dupilication 数据可以只保留一对数据,去除多余一致的测序片段,但是在变异检测过程中采取的是在bam文件中对比对到同一位置的duplication片段进行标记的方法,称为Mark Duplication。因为比较reads是否为duplication比较消耗资源,而采用标记的方法更加快速。一般duplication与其他过滤条件一起过滤,或者采用比对之后标记的方式。fastx-toolkit工具中可以去除duplicantion,但只能处理单端,因此用处不大。
fastx_collapser -v -i BC54.trimmed.fa -o BC54.collapsed.fa
截取头尾
illumina测序一般头部有些波动,尾部质量较差,如果想取出尾部,或者截取部分区域,可以使用seqtk trim功能,例如去除头部15bp,尾部15bp,可以使用-b 15 -e15。b为begin的意思,e为end的意思。
seqtk trimfq -b 15 -e 15 -Q SRR8651554_1.fastq.gz | head
数据过滤
有很多软件可以一次性完成数据过滤工作,包括去除adapter,去除低质量,去除N碱基,去除duplication,截取reads,常用的包括fastp,trimmomatic,SOAPnuke等,不过这些软件都各有优缺点,SOAPnuke很好用,但是去除adapter需要提供一个adapter reads ID的列表,从网上下载的数据没有这个,fastp利用固定adapter序列,但是不能去除duplication,trimmomatic选项参数太长,而且也不是很好用。这里推荐使用fastp。
fastp -i SRR8651554_1.fastq.gz -I SRR8651554_2.fastq.gz -o clean.1.fq.gz -O clean.2.fq.gz -z 4 -q 20 -u 40 -n 10 -f 15 -t 15 -F 15 -T 15 -h fastp.html
---------- END ----------
(添加作者微信,请注明单位姓名)
您可能还会感兴趣的
生物信息暑期班(北京站)开始报名
基因学苑文章列表(201906)上传数据,直接分析,1T内存服务器来了手把手教你生信分析平台搭建专栏合集生物信息重要资源站点合集不会编程,如何进行批量操作一个人全基因组完整数据分析脚本一个细菌基因组完整分析脚本如何在Linux下优雅的装X