fastq&fasta 处理

fastq格式文件处理大全

fasta格式文件处理大全

FASTQ

FASTQ

完整性校验

md5sum
保证文件在传输过程中保持完整,没有丢失内容,一般采用md5校验方式,目前测序公司给定的测序数据都带有md5文件,使用md5sum -c命令检测文件,如果返回OK,说明文件完整

md5sum -c SRR8651554_1.fastq.md5
md5sum -c SRR8651554_5.fastq.md5

文件统计

统计序列条数,碱基总数,reads读长分布等,可以使用seqkit
统计每条序列ATCG四种碱基组成以及质量值分布,可以使用seqtk comp
统计每个位点所有序列ATCG以及质量值分布,seqtk fqchk命令。结果可以绘制碱基质量以及含量分布图。

合并

可以使用seqtk mergerpe进行合并多个fastq文件,其实cat或者zcat也可以合并,不过seqtk的合并方式有一些差别,cat是将一个文件追加到另一个文件结尾,seqtk mergerpe是每次取文件一个单位合并。

过滤短序列

Ion Torrent,pacbio,nanopore测序的fastq文件序列长度并不相同,通常需要过滤较短的序列,例如过滤掉长度小于150bp的序列。可以使用seqtk seq或者seqkit seq进行操作。

转换为列表

可以使用seqkit fx2tb
转换为列表方便根据ID进行处理。
将四行数据转换为一行三列,可以使用awk列表处理程序来进行处理
使用tab2fx将处理好的列表转为fastq格式。

质量值转换

目前测序得到的fastq文件,都采用phred+33的格式,但是如果处理之前的文件,还有可能遇见phred+64的模式,一般软件中包含--phred33或者--phred64选项,当然也可以直接在两种质量值之间进行转换。

QC

fastqc绘制碱基含量分布图与碱基质量分布图,通过这两个图来判断fastq文件质量好坏。如果样品太多可以使用multiqc合并多个结果。

过滤

去adapter

接头adapter主要是指illumina测序时加入的P7接头与P5接头,理论上来说测序从测序引物开始,是测序不到接头的,但是由于部分打断片段或者由于adapter空载,会导致adpter自身连接,以上两种情况都会导致测序reads中包含adapter序列。adapter序列非基因组本身的序列,会干扰分析,因此需要去除掉。一种是直接给定adapter序列,与reads比对,比对上的就把整条reads去掉,另外一种是测序完成之后,给定一个adater list文件,其中包含了含有adapter序列的reads ID列表,给定一个阈值将这些reads去除掉。cutadapt可以根据给定的adapter序列进行过滤。

去除低质量

低质量主要是指去除测序质量错误率较高的位点,一般以Q20作为标准,如果一个碱基质量值低于Q20,则认为一个低质量,如果一条序列中低质量碱基达到一定比例,例如达到40%以上,则过滤掉此条序列。seqtk工具seq功能通过-q,-n可以将低质量碱基进行标记,例如替换为小写字母或者其他字符,但是不进行过滤,有专门的数据过滤工具。

去除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,但只能处理单端,因此用处不大。

截取头尾

illumina测序一般头部有些波动,尾部质量较差,如果想截取部分区域,可以使用seqtk trim功能

数据过滤

有很多软件可以一次性完成数据过滤工作,包括去除adapter,去除低质量,去除N碱基,去除duplication,截取reads,常用的包括fastp,trimmomatic,SOAPnuke等,SOAPnuke很好用,但是去除adapter需要提供一个adapter reads ID的列表,从网上下载的数据没有这个
fastp利用固定adapter序列,但是不能去除duplication
trimmomatic选项参数太长,而且也不是很好用。
这里推荐使用fastp。

排序

如果想对fastq格式文件进行排序,可以使用seqkit sort功能

抽样

有时候需要从全部文件中抽取一部分进行分析,因为测序出来的数据本身就是随机分布的,因此,即使从头到尾开始取数据,出来的也是随机的。当然还是可以随机抽样的。seqtk和seqkit工具都提供了抽样功能。

拆分数据

有时候需要将fastq文件拆成多份,或者根据固定模式进行拆分
例如测序时同一个lane的数据根据index进行拆分,
16S测序中,同一个文件中序列根据barcode进行拆分等。
seqtk split与split2可以用来拆分文件

转换为fasta

一些软件只支持fasta格式,例如只有fasta格式才能进行blast比对,将fastq转换为fasta有多种方法,awk都可以,这里使用seqtk和seqkit分别演示一下。

提取序列

fastq格式文件需要使用bgzip压缩,一般的fastq都是采用gzip压缩,需要重新处理文件才行。

合并两条序列

双末端测序的reads来自一条片段的两侧,如果文库大小比较小,测序读长比较长,例如pairend 300,insertsize 500,就有一些片段中间具有overlap区域,因此可以将两条reads进行拼接,连成更长的片段。有利于进行序列拼接,也具有更好的唯一性,例如16S序列分析中常用此步骤。有很多工具例如cope,flash,fastq-join等。

kmer计数

kmer是一段序列,一般将fastq文件按照固定长度(奇数),例如15,从头到尾按照步移数为1开始截取序列,例如1-15,2-16,3-17……然后对kmer频率进行计数,kmer分析可以用来估计基因组大小等,通过kmer频数分布也可以反应测序错误,或者杂合位点的分布。一般认为kmer频数为1的序列包含测序错误位点。可以使用jellyfish对fastq文件进行kmer计数。

fastq到bam

很多分析的第一步就是通过短序列比对将fastq文件转换为bam,包括变异检测,RNAseq等。一些测序仪直接输出bam格式文件,例如Ion Torrent,属于uBam格式。可以将fastq进行短序列比对的工具很多,这里以bwa为例。除了需要fastq格式,还需要fasta格式作为参考序列。

FASTA

image.png

统计

统计序列条数
grep wc
seqkit的统计功能

格式化

调整,例如一行ID一行序列,或者让每行显示固定长度内容。

逐条统计

统计每条序列长度
seqtk grep awk
统计长度并按照长度计算频数
seqtk grep awk sort uniq

成分统计

计算每条序列的长度以及ATCG碱基的组成,seqtk comp

提取序列

根据基因ID,提取序列。
方法一:利用samtools为fasta建立索引,然后提取
方法二:利用seqtk工具,给定一个ID列表

截取序列

如果想截取某条序列固定区域,需要给定一个bed格式ID列表

seqtk subseq

排序

seqkit sort

按照长度过滤

使用seqtk或者seqkit的seq

反向互补

seqtk是一步完成反向互补操作,
seqkit可以单独取反向序列,也可以单独取互补序列。
-r -p

转化大小写

seq -lu

查找重复序列

RepeatMasker

串联重复序列

trf

blast比对

blastn

建立bwa索引

fasta序列可以作为BWA比对的参考序列,比对之前同样创建索引。

翻译成氨基酸

可以使用一些在线工具例如ExPaSy,EMBOSS等。

多序列比对 多序列比对可以用于构建系统发育树,可以直接使用muscle,clustalw,mafft,megacc。

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

推荐阅读更多精彩内容