NGS分析入门:fastq文件的处理

拿到fastq文件的第一步是进行备份。使用gzip以及md5将文件压缩备份,并将其md5信息及实验信息写在readme文件中。

第二步去除barcodes。barcodes是在混合测序中运用到的一个区分测序样品的手段,使用不同的barcode,就可以将不同来源的样品区分出来。barcode是一段很短的oligo,比如ATCACG等。它很有可能是由一段酶切位点的一部分加上几个碱基组成。因为barcode半不是真实的序列,而是人为加上去的,所以需要把它从测序结果中去除。在去除的同时,需要你可以依据自己的实验设计,使用barcode信息将不同来源的测序结果分割开。

第三步是了解测序质量。可以使用fastqc工具。这一工具有用户界面,使用十分简单。

第四步是对fastq文件依据测序质量进行筛选。

以上这些步骤有很多可以使用的工具软件,这里推荐的是FASTX-Toolkit

以上这些步骤都是最基础的步骤,一开始是fastq文件,最终得到的依然是fastq文件。

接下来的工作,可能会用到,那就是将fastq文件转换成unique sequence count文件。
其中一种文件是tab间格的文本,格式为:

TTGCCTGCCTATCATTTTAGTGCCTGTGAGGTGGAGATGTGAGGATCAGT    2345
GATTTGTATGAAAGTATACAACTAAAACTGCAGGTGGATCAGAGTAAGTC    2241

这种文件的好处是文件小,适合网络传输,但其缺点是完成丢失了测序质量信息。可以使用miRanalyzer提供的工具完成。

还有一种是miRdeep使用的fa格式,它的每一条序列名都是由name_uniqueNumber_xNumer这种格式组成,比如:

>PAN_123456_x969696
ATACAATCTACTGTCTTTCCT

这种格式可以使用miRDeep2中的mapper.pl来完成。

56 thoughts on “NGS分析入门:fastq文件的处理

  1. 版主,你好!我是个菜鸟,看到你的BLOG受益匪浅,我有问题想请教一下:我在ENA上下载了一些fastq文件数据,我现在有一段200bp序列,我想拿这段序列和 fastq文件的数据blast,请问这种想法是否可行,如果可行的话 ,能否帮忙提供一些解决思路。

    1. 当然是可行的。这就是序列比对的基本操作。你可试着搜索一下align two or more sequence,比如说NCBI的工具http://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastSearch&PROG_DEF=blastn&BLAST_PROG_DEF=megaBlast&SHOW_DEFAULTS=on&BLAST_SPEC=blast2seq&LINK_LOC=align2seq, 或者clustal omega, http://www.ebi.ac.uk/Tools/msa/clustalo/, 又或者专业人员喜欢的phylip,http://evolution.genetics.washington.edu/phylip.html

  2. 版主您好,我最近在尝试对miRNA-seq结果进行分析。我现在想将我的原始数据fastq文件转换成unique sequence count文件。这种文件的好处是文件小,适合网络传输,但其缺点是完成丢失了测序质量信息。可以使用miRanalyzer提供的工具完成。我尝试打开版主您给的链接,但页面无法打开,我又登陆miRanalyzer主页寻找相应的工具,也没有找到╮(╯_╰)╭。麻烦版主可不可以把相应工具的地址再公布一下,或者有没有其他的工具可以完成fastq文件格式的转换~~O(∩_∩)O谢谢!

    1. 我这边可以打开。你可以再试试:https://www.mdc-berlin.de/8551903/en/research/research_teams/systems_biology_of_gene_regulatory_elements/projects/miRDeep

  3. 斑竹您好,请问这NGS的流程就是关于illumina的么?有没有Ion proton的RNA-seq数据处理的介绍?还是说和illumina的差不多?

  4. 谢谢斑竹。我是新手,请教一下。 Ion不是有自己的一套数据分析系统软件么?如Torrent Suite 什么的? 是不是说illumina关于NGS的流程 如tophat-cufflink-cuffdiff不太适合Ion平台的数据?否则为何它们自己单独开发了一套?谢谢

    1. 说句实在话,对于Ion的系统,我没做过。不过,最有可能的是它的mapping可能会不一样,其余的步骤应该不受影响。但以我的理解,即使是mapping,也应该差不多才是,毕竟输出的都是reads.但是,如果你是做的原核生物的话,那应该会和真核生物很不一样。

  5. 您好,我想请问一下数据下机之后除了通过barcode区分样本信息及去除Barcodes之外还有其他要做的吗?
    1、barcode序列是不是包含在read的序列里?
    2、barcode 和 insert size 之间不是还有接头序列吗,接头序列有没有测到read中,有的话是不是也要去掉?
    3、如果不进行质量过滤,我们拿去比对的fq文件是不是已经把barcode和adapter除去掉了?

    1. 1. barcode是否包含在read的序列里?是的,都会包含在read里,除非你错误的使用了测序方法。
      2. barcode 和 insert size 之间不是还有接头序列吗,接头序列有没有测到read中,有的话是不是也要去掉?这与你测序的方法有关。通常状况下,人们并不是非常需要超出正常测序长度的reads。比如可以测50nt的测序仪通常不会提供超过50nt的reads。但是如果你将barcode整合至adapter中的话,需要仪器逆向测序出barcode,这会非常费钱,几乎没有人这么做了。正常的测序是adapter-barcode-insert-adapter这样。通常而言,第二个adapter是测不到的。仪器在测序时,因为最开始的碱基是测不到的,之后会测到adapter的一部分,但这一部分会在输出前就被切除,只报告barcode-insert–的部分,除非你的barcode包含在adapter中。通常的测序长度小于片段长度,所以第二的adapter测不出来。但是如果你是很小的insert并做了环化的话,有可能会测出第二个adapter,这时需要将其去除,但这种情况非常少见。所以,普通的测序不需要去除adapter。
      3. 如果不进行质量过滤,我们拿去比对的fq文件是不是已经把barcode和adapter除去掉了?当去除barcode之后(通常使用demultiplex这个词来代表这个过程),得到的reads就是可以拿来比对到基因组上的reads了。

  6. 在服务器上安装FASTX-Toolkit时,自己不是root,安装不了,有其他不需要是root的软件可以进行质量删选么?

    1. 我想你不用安装,可以直接下载相应系统下的二进制文件运行就可以了。如果不行,你可以下载原代码自己编译一下,具体方法你可以在网上搜索一下。

  7. 现在的Illumina RNA-seq数据还需要去除barcode吗?因为根据FASTX-Toolkit 的 –bcfile选项来看,还需要设置mybarcode.file 这个从哪里来呢?

    1. 这取决于你的实验。如果你的实验加入了barcode,你就需要去除。如果没有,那就不需要。对于barcode文件的准备,基本上就是类似这样的一个文本文件:
      BC1 GATCT
      BC2 ATCGT
      BC3 GTGAT
      BC4 TGTCT
      每行一个barcode,前面是名称,后面是barcode序列,中间使用tab符号间隔。

  8. 我从ncbi下载了一些chip-seq数据,illumina的,reads是35bp。我问我怎么找这个数据的barcode和adaptor呢?

    1. 你需要从原始文献中去查找。ncbi中会要求上传一些资料,但是没有要求上传barcode和adaptor。另外,你也可以email作者,一般会告诉你。如果都拿不到,只能靠猜了。

      1. 谢谢回复,我已经给作者发邮件了,如果有问题想请教您能给我留个邮箱吗?
        如果不方便公开可以通过邮箱告诉我吗?

  9. 你好!博主,我根据你的方法正在做转录组测序数据分析,但是,我用fastx-toolkits软件时,发现fastq_clipper无法处理PE数据,因为pe=pairend 数据有两个adapter,而fastq_clipper只能一次行处理一个adapter,所以,我该怎么办?我查了资料发现fastx-toolkits较好地处理50bp一下的数据,超过100bp的数据,其结果就不是很好了(不知道真假:https://www.biostars.org/p/74111/)博主,请查证!请推荐合适我的软件!谢谢!!

  10. 博主,您好!
    我手中现在拿到了单端测序的原始数据,没有经过处理
    reads中应该是包含了index与adapter
    请问,我如果使用fastx toolkit来剔除index,adapter序列,可不可以直接把index-adapter序列的反义链作为参数使用?(index序列紧接着adapter序列,或者说我的理解错误了,这种情况实际是index包含在adapter里面了?)
    谢谢指教

    1. 你需要了解的是测序是如何做的。因为index和adapter有两种可能的组合的方式,见http://www.umassmed.edu/uploadedFiles/nemo/Landing_Pages/Indexing%20and%20Barcoding%20for%20Illumina%20NextGen%20Sequencing.pdf。我考虑以后有机会写一下这个方面的中文介绍。

  11. 老师您好,看了您的博客非常有帮助!请问我得到了几个fastq文件合并后的fastq文件,有什么软件可以看新的fastq的>Q30%和PF reads%吗?得出来类似做完bcl2fastq后得到的report那样?非常感谢!

  12. 你好,感谢你的博客!我有个问题是,对序列质量进行过滤以后,还需要对含“N”数目的序列进行过滤么?这两者有对应关系么?

    1. 没有这个必要了。有时候个别的N并不影响测序的下游分析。你需要查看文件确认一下N是否出现在5’端,如果在,你需要对5’端进行截断处理。

  13. 博主您好,我有个问题想要请教您。在我对paired-end数据质量过滤(quality score > 20)之后,发现包含”N”的reads有不少,而且有些reads是连续几个N,有的出现在reads中间,有的在中间偏后的位置,如果是这种情况,是否需要在mapping之前过滤掉这样的reads?

    1. 不需要。出现在末尾说明是片段小于测序的长度。出现在中间,可能是测序操作的问题。但一般出现在中间的都只有一两个碱基,不过它们是很影响mapping的。

  14. 老师您好,刚刚接触NGS很多东西还不太明白,想请教一下您,从NCBI sra上下载的数据转换成fastq格式之后如何查看是否有接头序列?

  15. 版主,我在使用trim-galore进行adaptor去处后发现fastqc文件的Overrepresented sequences选项中序列经比对都是一些ribosomal RNA片段。但是当时准备文库时候使用的 都是poly-T oligo纯化的polyA mRNA,而且最后pcr后 都是经过bioanalyzer验证过的没有问题。不知道这些ribosomal RNA哪里来的,导致fastqc的结果很差(例如,Per base sequence content [FAIL] Per base GC contentSequence Duplication Levels [FAIL] Overrepresented sequences[WARNING] Kmer Content)。

    1. 这种情况我遇到过不少。其实就是rRNA污染。rRNA污染并不象你想象中那么容易去除。你是不切片段的时候切得太宽了,或者片段控制得太小了?

      1. 是的,构建文库时片段化过程的时间设置为1min(default是8min),期望片段较长便于覆盖跟拼接!我听其他人建议说如果后期分析着重于差异比较转录组分析,rRNA可以相互抵消,是这样么?

        1. 其实rRNA对最后的影响并不是很大,但是很消耗运算量,所以人们希望先过滤掉它们。而且,rRNA/tRNA等因为量大的话,会占用很多reads,客观地减少了coverage。

      2. 另外,版主,面对rRNA污染,在制备文库时有什么比较好的解决办法么?(提取RNA时选择mRNA提取方法而不是totalRNA?)还有,比如像我这样样品中已经混有rRNA,如何在后续分析时进行软件去除?(blast?)。。。问题有点儿多哈 :)

  16. 博主 您好,我现在打算用miRanalyzer 做microRNA的定量分析,当我用groupReads.pl将fastq文件转换成unique read count(rc)文件后,将其作为input file运行miranalyzer后,但是显示错误,用mfa格式也是。我尝试将dapter去除和不去除两种情况,但是都是错误的。想请问这是什么原因呢?谢谢啦
    rc file(未去adpter):
    TCCTGTACTGAGCTGCCCCGATGGAATTCTAGGGTGCCAAGGACTCCAGT 1
    TATTGCACTTGTCCCGGCCTGTGGAATTCTCGGGTGCCAAAGAACTCCAG 1
    ATATAGATATCATCAGAGCTTGGAATTCTCGGGTGCCAAGGAACTCCAGT 1
    TCCTGTACTGAGCTGCCCCGATGGAATTCTCGGGTAAGGAACTCCAGTCA 6
    rc file(去除adpter):
    TCCTGTACTGAGCTGCCCCCGG 1
    TAGGCACCGTGGCTCACGCCT 1
    CGATCTTCTCACCTTGGCCTT 1
    AACAGATCCATCTCTAGACATG 1
    CAAAGTGCTGTTCGTGCAGGTAT 1
    AAACTGAAGAGTTTGATCCT 1
    ATCCTGTACTGAGCTGCCCCGAAA 12

      1. 是的 我还尝试了自己用linux命令建立unique read count的txt file, 也是没办法运行,命令如下: awk ‘(NR%4==2){print}’ reads.fastq|sort|uniq -c|awk ‘{print$2″\t”$1}’ > read_for_miranalyzer.txt

        显示 An error occured with you job: 62512359 Please report it indicating the jobID错误报告

  17. 楼主,你好, 请问一下,分析一个 50X-60x的 人 全基因组测序原始fastq,到 vcf, 一般要多长时间呢?我要用10来天,除了拆分文件,并行分析,还有什么其他方法可以缩短时间呢?

  18. 版主你好。我是新手,剛接觸NGS對於檔案格式不是很明白,想請問.fq跟.fastq file 都是raw data 的形式,那麼它們的差別是什麼?CLC genomic workbench 是否只支援.fastq檔?那我該如何處理.fq檔?

  19. 版主你好。我是新手,剛接觸NGS對於檔案格式不是很明白,想請問.fq跟.fastq 都是raw data 的形式,那麼它們的差別是什麼?CLC genomic workbench 是否只支援.fastq檔?那我該如何處理.fq檔?

    1. .fq和.fastq只是文件扩展名不同,其格式要求是一致的。如果软件只认.fastq的话,你改个名字就好了。有人会混淆.fa和.fq,.fa没有测序质量,所以每条记录比fastq少两行。fastq是每条记录四行。你可以打开文件看一眼。

      1. 好的,感謝!!
        其實原本是擔心資料太大,個人筆記型電腦容量不夠大,才沒有解壓縮打開文件。

  20. 欧老师,有两个问题想请教一下,需要处理一些PE的原始数据。看过之后发现,有很多序列中出现成片的N,如下(R2),这种事什么原因呢?需要继续做过滤么?还是直接认为有问题呢?
    @HWI-ST754:95:D0GEWACXX:3:1101:6314:2002 2:Y:0:CAGATC
    GNANTGGATNNNNNNNNNNNNNNNNNNNNNNNNGNNNNNNNNNNNNNNNAGCATGACCTCACGCCTTNTATAAAAGCCCTCCCTTAGTCTTTTGATTATCC
    +
    C#1#4ADDH########################0###############–5?CHEEFFFFDDEDDD#,5-,9??CAAB@6?8;>:ACDCCC9:44@@?@###############
    @HWI-ST754:95:D0GEWACXX:3:1101:6285:2007 2:N:0:CAGATC
    ATTATAAAAANNNNNNNNNNNNNNNNNNNNNNAANNNNNNNNNNNNCNNATAAAAAATATGTGTGGATTATTAAACAAAAACTACATATTATTCATTAGAA
    +
    CCCFFFFFHH######################00############-##.-BGHHHFFFFE>CEEDDDDDDEFFDDDDDDDDDDDDDDEEDEEEEEEEDDD

    第二个问题是,看到很多相连的文件中,N都出现在第三个碱基位,这可能是什么原因造成的呢?如下:
    @HWI-ST754:95:D0GEWACXX:3:1104:15160:39682 1:N:0:CAGATC
    AGNACCACGGTGTTAATTACGTATGTGAAATAACGGCAATGGGAATTTGCTTGGATTGCAGATACGCTAGAACCACGGTGTTAATTACGTATGTGAAATAA
    +
    @@#4ADDFGFFHFIJJIIIJIFHIHHIJJIEHJEHIGIJJIIIIEHHIIFJIIIIIJIEGGGHHHHFEFDCEEDD?BBBDDDDEEEEDDBCDDBDEDACDD
    @HWI-ST754:95:D0GEWACXX:3:1104:15135:39688 1:Y:0:CAGATC
    CTNTATCCTTAGAAAACAATTGTTTGCACTGGATGAGAAAACACAAATGAATAATAAAACGTAAAAATAGAGAAAAGAACAGAAATTCTAAAAAAGAACAC
    +
    ;<#-).=3=??5)0639>?:.;)../-??@;?AACC>CBB(9<ACCCC:?>AA>@CCDAC4:>>>ACCC@DACB
    @HWI-ST754:95:D0GEWACXX:3:1104:15047:39694 1:N:0:CAGATC
    GCNTGATTCTAGGAGAATGCTGGCGGACGCGAACACTTTATCAAGACAATTTCTGCTGTTTATACTATGGCGCGGGCGTCCGCCAGCATTCTCCTAGACTG
    +
    ??#44ADD??<4++<A:E@DEA>,:;:>>((5(055&&&)5-0052)009(:>:A>AAAAAA>

    多谢!

    1. 对于这些情况,我通常不做处理。如果mapping效率太低的话,对于第二个问题中出现的情况,可以对所有的reads都删除前三个碱基来提高mapping的效率。

  21. 老师您好,我是一个新手,dna测序结果是几个g的fq文件,请问该用什么软件去处理这些文件,让我可以分析里面这些序列中的相同序列和相似序列呢?我用notepad++打开了fq文件,可是因为数据太多不能复制到doc文件中进一步处理,看了你上面的文章尝试去用你推荐的软件可是没能成功,还不会用。

发表评论

电子邮件地址不会被公开。 必填项已用*标注