|
读而思deepreader 关注微信号每天收听我们的消息 读而思为您推送精品阅读 读而思deepreader 关注微信号每天收听我们的消息 读而思为您推送精品阅读 读而思deepreader 关注微信号每天收听我们的消息 读而思为您推送精品阅读 读而思deepreader 关注微信号每天收听我们的消息 读而思为您推送精品阅读 生信媛 关注微信号, 获取我们不定期的高质量文章推送 写在前面今天徐爷提了一个问题：
天地本无心认为：
然后徐爷问了一下公司，得到如下回复 这个时候生物女博士（本文作者）默默给了一张截图 于是我决定重新把这篇文章发出来，让小伙伴重新复习一下（PS，我是日更） 1 为什么需要进行QC ?测序质量问题包括[1]：
这些问题的存在会影响reads mapping到参考基因组、组装、表达量估计等。有的问题可以通过过滤、剪切、错误或偏差纠正来解决。而有的问题我们是没办法解决的。但至少，我们在解释分析结果的时候，能意识到这些问题的存在[1]。 导致质量问题的原因：
网页链接：BGI常见问题解答：http://www./science/cat-91-95-409.html 质控并没有一个确定标准。在质控过程中，会去除一些低质量reads、切除质量差的碱基、在切除质量差的碱基后还可能会去除一些长度过低的reads……而这些处理，其实相当于我们损失了一些信息。而由于信息的损失，也可能会造成最后结果有偏差。 质控是信息量和信息准确度之间的一个平衡。所以，质控并见得是越严格越好。应该根据自己实际情况的需要来选择质控的严格程度。一般来说，对基因组知道的越少，质控应该越严格[2]。 质控分两个阶段[2]：
2. 序列质量检测及可视化软件
2.1 FastQC的安装和使用FastQC简介[1]FastQC是一个Java程序。既有图形操作界面（graphical user interface，GUI），也可以使用命令行操作。图形界面的高通量数据分析软件Galaxy[3]和Chipster[4]也将其整合其中。 FastQC运行速度比较快。支持的输入文件格式有[5]：
FastQC安装以下安装流程仅供参考，读者可根据自己情况选择别的安装方式。 # 进入你软件源码存放的文件夹。cd ~/software# 下载软件curl -O http://www.bioinformatics./projects/fastqc/fastqc_v0.11.5.zip# 或使用wget：wget http://www.bioinformatics./projects/fastqc/fastqc_v0.11.5.zip# 用Linux的: 用sudo apt-get 安装和解压。unzip fastqc_v0.11.5.zipcd FastQCls -l# 使fastqc执行。chmod +x fastqc# 把可执行的fastqc文件链接到入~/bin目录下ln -s ~/software/FastQC/fastqc ~/bin/fastqc# 或者把可执行的fastqc文件直接复制到到入~/bin目录下# cp ~/software/FastQC/fastqc ~/bin/fastqc# 测试软件是否安装成功fastqc -h# 出现帮助文档则为安装成功。在下边这一步中，需要注意的是，你的 # 把可执行的fastqc文件放入~/bin目录下ln -s ~/software/FastQC/fastqc ~/bin/fastqc如果没有添加，可以按下边命令进行添加 # 将~/bin 文件夹加到PATH# Mac下:echo 'export PATH=~/bin:$PATH' >> ~/.profile# 使设置生效source ~/.profile# 在Linux 下这么做:echo 'export PATH=~/bin:$PATH' >> ~/.bashrcsource ~/.bashrc添加完了之后再按照上边的命令进行软件安装。 此外，也可以直接含有可执行命令的目录 当然，关于软件安装还有许多其他的方法，比如使用conda等工具，这里不赘述。 关于如何使用conda，请参考洲更同学的 FastQC命令使用使用下边的命令对reads.fastq.gz这个压缩格式的fastq文件进行质量检测，并且对每个测序文件生成两个结果文件: html格式的报告文件以及一个zip压缩文件。 # 使用prefetch从NCBI下载测试数据集。下载的数据会放在目录~/ncbi/public/sra/prefetch SRR1553607# 创建文件夹并把sra文件拆分为两个测序文件（双端测序）mkdir QCcd QCfastq-dump --split-files SRR1553607# 获得2个数据文件 SRR1553607_1.fastq 和 SRR1553607_2.fastqls # 查看文件# 使用fastqc进行质量检查fastqc *fastq# 此时会在当前文件夹生成结果文件。# 假如你想把所有的结果文件统一放到一个文件夹：mkdir resultsfastqc -o results *fastq可以从百度云下载该测试数据：链接: http://pan.baidu.com/s/1o8r8Lz8 密码:xm9v 关于从NCBI上下载数据，更多信息请参考我们公众号的文章《从NCBI下载测序数据 | 也许是目前最详细的版本》 。 FastQC给出的报告结果包括以下几个内容：
并且给出一个判断：通过（pass）、警告（warn）和不通过（fail）。不过，软件设置的判断阈值未必适合你的数据。 所以即便有一些项目出现警告或者不通过也不用过于担心。 比如，RNAseq的数据在序列重复水平（Sequence Duplication Levels）这一项不通过就很正常。 而RNAseq数据由于使用了随机六聚体引物，在reads开头会出现碱基成分偏离，所以Per base sequence content一项也会出现警告或者不通过的信号，而且在5’端会产生比较多的K-mer。这是一个真正的技术偏差,也无法通过修剪reads来纠正,不过在大多数情况下似乎不会对下游分析产生负面影响[15]。 更多fastqc结果文件的解读可以参见微信公众号“生信菜鸟团”的文章《测序数据质量控制之FastQC》或者FastQC官方帮助文档[14]。这里不再赘述。 2.2 PRINSEQ的安装和使用PRINSEQ有网页版本, 也有命令行版本。Chipster GUI中也有PRINSEQ。 PRINSEQ输出的报告内容包括[1]：
PRINSEQ可以处理的文件类型为非压缩的FASTQ、FASTA以及QUAL[1]。不支持压缩格式的fastq。 安装# 进入你软件源码存放的文件夹。cd ~/software# 下载软件curl -OL http://downloads./project/prinseq/standalone/prinseq-lite-0.20.4.tar.gz# 解压tar zxvf prinseq-lite-0.20.4.tar.gz# 进入文件夹cd prinseq-lite-0.20.4ls # 查看目录里的内容# 把该目录下的三个perl脚本(prinseq-lite.pl、prinseq-graphs.pl和prinseq-graphs-noPCA.pl )变得可执行。chmod +x *.pl# 把目录 ~/software/prinseq-lite-0.20.4/ 添加到环境变量# Mac下:echo 'export PATH=$PATH:~/software/prinseq-lite-0.20.4/' >> ~/.profile# 使设置生效source ~/.profile# 在Linux 下这么做:echo 'export PATH=$PATH:~/software/prinseq-lite-0.20.4/' >> ~/.bashrcsource ~/.bashrc# 测试软件是否安装成功prinseq-lite.pl -h# 出现帮助文档则为安装成功。通过上面的安装，prinseq-lite.pl是可以正常使用的。但是由于画图还需要别的perl模块，所以需要进行模块安装才能正常使用。 非root权限下的perl模块安装：查看软件目录下的README文档可以知道：
prinseq-lite.pl无需安装perl模块。
非root用户运行下面的代码获取自己的私人cpan下载器cpanm。 wget -O- http:// | perl - -l ~/perl5 App::cpanminus local::lib# 不知为何，在这一步我需要运行两次才能安装好。eval `perl -I ~/perl5/lib/perl5 -Mlocal::lib`echo 'eval `perl -I ~/perl5/lib/perl5 -Mlocal::lib`' >> ~/.bashrcecho 'export MANPATH=$HOME/perl5/man:$MANPATH' >> ~/.bashrcecho 'export PATH=$PATH:~/perl5/bin' >> ~/.bashrcsource ~/.bashrc# ~/.bashrc可能需要更改为.bash_profile, .profile, etc等等，取决于你的linux系统！ # 然后你直接运行cpanm Module::Name，就跟root用户一样的可以下载模块啦！注：安装cpanm的代码根据perl模块安装大全略有修改。 用cpanm基本能安装上大部分的模块。但是prinseq需要的模块，有的安装起来真心心累。比如在安装Cairo的时候，需要系统安装了cairo这个画图的工具库。而cairo画图工具库的安装需要依赖另外两个库：libpng和Pixman。没有root权限，安装很难。而且看质量情况一般用FastQC也比较足够，所以果断弃坑。 命令使用PRINSEQ的质控报告可以是文本格式或者html格式，如果想输出html格式，需要把前面提到的模块安装好。 输入html格式需要分为两步[1]：
# 一张张的png图格式prinseq-graphs.pl -i graph -png_all -o QCreport# html格式prinseq-graphs.pl -i graph -html_all -o QCreport3 reads数据的质控这部分内容主要来自：RNAseq Data Analysis: A Pratical Approach 3.1 碱基质量碱基质量一般用Phred值表示。碱基错误的可能性，用log~10~进行变换，再乘以-10。举例来说，某个碱基有1%的概率是错误的，那么它的质量值为q = −10*log~10~(0.01) = 20。质量值的范围一般是0到40。在fastq文件里，为了节省存储，质量值是用ASCII字符表示的。目前的fastq文件都用的Sanger编码方式：第33个ASCCII字符表示质量为0。但是在1.8版本之前的fastq文件是以第64个ASCCII字符表示质量为0。如果你不知道你的数据是哪种编码系统，用FastQC可以检测出来，如图。而Trimmonatic软件可以帮你把fastq的质量编码方式进行转换[1]。 关于测序质量编码方式的讲解，可以查看《小白生信学习记3》。
3.2 低质量碱基的reads的处理方式含有低质量碱基的reads的处理方式包括过滤和切除。 如果你决定要过滤或者切除双端测序的reads，则需要选择一个能保持双端序列顺序对应的工具。因为当把reads mapping到参考序列时，比对软件会以相同的顺序在两个测序文件里寻找配对的reads。注意，不单是过滤，切除也可能会对顺序造成影响：当一些reads被完全切除或者切除得太短而被过滤掉，从而导致顺序的变化[1]。 3.2.1. 过滤：去除整条readsTrimmomatic, FastX,以及PRINSEQ都可以通过质量过滤reads。
Trimmomatic的安装安装方法参考了《Biostar课程7、8-二代测序质控：fastqc、prinseq、trimmomatic、seqtk》 # 安装trimmomaticcd ~/softwarecurl -O http://www./cms/uploads/supplementary/Trimmomatic/Trimmomatic-Src-0.36.zipunzip Trimmomatic-Src-0.36.zipcd Trimmomatic-0.36# 可以到这查看一下手册：http://www./cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf# 运行的命令比较长java -jar trimmomatic-0.36.jar# 我们写一个脚本来替代我们运行。# 你也可以通过命令行来达成echo '#!/bin/bash' > ~/bin/trimmomaticecho 'java -jar ~/software/Trimmomatic-0.36/trimmomatic-0.36.jar $@' >> ~/bin/trimmomaticchmod +x ~/bin/trimmomatic# 安装成功。trimmomatic# 也可以使用conda进行安装。conda install trimmomaticTrimmomatic[1]Trimmomatic是一个多功能的基于JAVA的reads处理工具。Trimmomatic既有命令版本，也被Galaxy 和 Chipster图形界面软件整合了。Trimmomatic的输入和输出文件都是FASTQ（压缩或者非压缩的都可以）Trimmomatic是多线程的，所以运行起来很快。 Trimmomatic的功能：
Trimmomatic可以一条命令，按着你给的顺序把多个步骤完成。 下边这条命令是用Trimmomatic过滤双端reads，设定read的碱基平均质量阈值为20 (AVGQUAL:20)，输入和输出文件都为压缩的[1]。 trimmomatic PE -phred33 SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz paired1.fq.gz unpaired1.fq.gz paired2.fq.gz unpaired2.fq.gz AVGQUAL:20# 如果没有按照上边那样写一个脚本进行命令名称替代，则按下边写:java -jar trimmomatic-0.36.jar -phred33 SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz paired1.fq.gz unpaired1.fq.gz paired2.fq.gz unpaired2.fq.gz AVGQUAL:20-phred33 表质量编码体系是示Illumina 1.8之后的版本。质量编码体系是Illumina 1.8之前的版本，需要改为-phred64。不过，从0.32版本，软件会自动识别是33还是64。所以 -phred参数可以不加： trimmomatic PE SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz paired1.fq.gz unpaired1.fq.gz paired2.fq.gz unpaired2.fq.gz AVGQUAL:20用PRINSEQ进行同样的质量过滤[1] ： prinseq-lite.pl -fastq SRR1553607_1.fastq -fastq2 SRR1553607_2.fastq -min_qual_mean 20 -out_good qual_filtered -out_bad null –no_qual_header –log -verbose# 同样，如果质量编码体系是Illumina 1.8之前的版本需要加上-phred64：prinseq-lite.pl -fastq SRR1553607_1.fastq -fastq2 SRR1553607_2.fastq -phred64 -min_qual_mean 20 -out_good qual_filtered -out_bad null –no_qual_header –log -verbose
3.2.2. 切除：从read一端切除低质量碱基切除方案1：设定切除长度或者是切剩余的长度。假如在reads末端发现低质量碱基，最简单的方法就是从某端开始切除，并设定切除reads的剩余长度或者给一个需要切除的长度。 但是这个方法会把质量好的那部分序列也去除了。如何减少高质量序列的损失？ 切除方案2：切除低于质量阈值的碱基（从3’或者5’端开始，逐个切除低于你设定质量阈值的碱基。）
用Trimmomatic从3’端开始切除低质量碱基（TRAILING表示3’端开始，如果想从5’开始，则改为LEADING），碱基质量阈值为20 （TRAILING:20），过滤切除后长度小于50个碱基的reads（MINLEN:50）。切除和过滤参数需要按正确顺序排。 trimmomatic PE SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz trim_paired1.fq.gz trim_unpaired1.fq.gz trim_paired2.fq.gz trim_unpaired2.fq.gz TRAILING:20 MINLEN:50切除方案3：使用滑窗的方法（设定一个滑窗，将窗口的碱基质量与设定的阈值进行比较。）
注意：
由于有的reads有可能会出现中间质量很低，而两端比较高的情况。而测序质量一般都是从5’高而3’低，所以使用5’端开始滑窗的方法通常会导致质控后产生的reads更短。导致低于read长度阈值的序列更多，使得被过滤掉reads也随之更多。而这个问题对于双端测序影响会更大，因为其中一端测序的read被过滤，在另一端测序文件里对应配对的序列也会被去除。 除了可以选择reads扫描方向，PRINSEQ还有其他灵活的设置：Trimmomatic允许用户设置窗口大小，但只能使用窗口内碱基质量的平均值作为窗口的质量值，而PRINSEQ还允许设定滑窗的步长、允许选择使用平均值或最小值来跟阈值作比较。 下边的命令是使用Trimmomatic软件，设置窗口大小为3碱基，从5’开始切除reads，read的平均质量低于20（SLIDINGWINDOW:3:20），开始切除后边的碱基序列。过滤长度小于50碱基的reads（MINLEN:50）。使用更大的窗口可以减轻切除的力度。 trimmomatic PE SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz win_paired1.fq.gz win_unpaired1.fq.gz win_paired2.fq.gz win_unpaired2.fq.gz SLIDINGWINDOW:3:20 MINLEN:50下边的命令是使用PRINSEQ软件，设置窗口大小为3碱基，从3’开始切除平均质量低于20的序列（高于20时停止）。同样过滤长度小于50碱基的reads。 prinseq-lite.pl -trim_qual_window 3 -trim_qual_type mean -trim_qual_right 20 -trim_qual_rule lt -fastq SRR1553607_1.fastq -fastq2 SRR1553607_2.fastq -out_good window -out_bad null -verbose -min_len 50 -no_qual_header切除方案4：另一种滑窗的方法是“running sum”，最初是在BWA比对软件中采用，所以通常也被称为“BWA quality trimming” 。从右端（3’端）开始扫描reads，将每个碱基的质量值与给定阈值进行比较，并计算累计差异的加和。在累计“badness”最高的read位置切除。这个方法被Cutadapt软件采用。 最后，Trimmomatic提供了一种适应性的切除方法，叫MAXINFO。其目标是在尽可能在保留reads长度和去除错误碱基序列之间取得一个平衡。 它包括两个参数：
可以通过严格度参数来控制这个平衡，严格度是从0到1的值，数值越高越严格。参数值低 (<0.2)有利于保留更长的reads，而高(>0.8) 则保证reads的准确性。 MAXINFO过滤的原理参见：
下边的MAXINFO切除命令是设置目标长度为50，严格度为0.7[1]。 trimmomatic PE SRR1553607_1.fastq.gz SRR1553607_2.fastq.gz MI_paired1.fq.gz MI_unpaired1.fq.gz MI_paired2.fq.gz MI_unpaired2.fq.gz MAXINFO:50:0.7 MINLEN:503.2.3 不确定碱基（Ambiguous Bases）在测序时，如果一个碱基无法判定是哪种，则记为N。不同的组装软件和比对软件对于这种不确定碱基的处理方式不一样：有的把N用4种碱基进行随机替代，有的则是用固定某个碱基替代N。由于Ns会导致错误的组装和mapping，含N多的reads需要被去掉。 PRINSEQ的质控报告会给出一个Ns含量的图，可以看reads含N比例的一个频率分布情况。
使用PRINSEQ过滤含Ns过多的reads： prinseq-lite.pl -fastq SRR1553607_1.fastq -fastq2 SRR1553607_2.fastq -ns_max_n 2 -out_good nfiltered -out_bad null -no_qual_header -log -verbose
保留的reads放在 3.2.4 接头（Adapters）测序时候，一般会使用测序接头序列，而接头序列的存在可能会导致基因组组装和转录本组装出现问题，所以需要在分析数据之前去除掉。此外，其他的tags，比如multiplexing identifiers以及引物也需要被去除[1,2]。 去接头是很有挑战性的[1]：
fastqc支持你把接头序列加到adapter文件中，进行接头分析。 查看接头序列。 cat ~/software/FastQC/Configuration/adapter_list.txt假如接头序列未知，可以用TagCleaner软件预测。而FastQC的k-mer overrepresentation图和PRINSEQ的Tag Sequence Check图也可以检测出是否有接头。 Trimmomatic, FastX, TagCleaner以及Cutadapt都可以去除接头序列。 其中Trimmomatic, TagCleaner和 Cutadapt：
注意，如果你要把general、quality以及adapter trimming 同时放到Trimmomatic的一个命令中，你要把切除接头放到第一步，因为识别完整的接头序列比识别部分接头序列要容易得多。 我们也可以自己创建接头序列文件[2]，或使用Trimmomatic自带的： 创建接头序列文件: echo '>illumina' > adapter.faecho 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC' >> adapter.fa使用trimmomatic去除另一个sra文件SRR519926的接头: # 单端去除接头trimmomatic SE SRR519926_1.fastq output.fq ILLUMINACLIP:adapter.fa:2:30:5# 查看质量情况fastqc SRR519926_1.fastqfastqc output.fq去除接头后的reads存放在 左边的图示去除接头之前，右边的图是去除接头之后。
使用bbduk进行切除接头: bbduk.sh in=SRR519926_1.fastq ref=adapter.fa out=bbduk.fq大部分工具都允许多个步骤放到同一条命令运行。 # 双端模式trimmomatic PE SRR519926_1.fastq SRR519926_2.fastq trimmed_1.fq unpaired_1.fq trimmed_2.fq unpaired_2.fq SLIDINGWINDOW:4:30 TRAILING:30 ILLUMINACLIP:adapter.fa:2:30:5先后顺序一般来说很重要，因为后面的处理依赖于前面的结果。 先以质量值切除，再切adapter，还是反过来更好？由于adapters有重叠区，如果我们先以质量切除，那么有可能因为切除了部分碱基而导致无法识别出adapter，另一方面，如果碱基质量很差，也可能无法比对上adapter而导致adapter无法被检测出来。通常，这个顺序不会导致特别大的差异。不过，对于有问题的数据，最好还是先试验一下[2]。 3.2.5 随机六聚体引物导致的sequence-Specific Bias以及错配在文库构建时，RNA是片段化的，随机引物六聚体会被作为引物来反转录产生cDNA，然后再从末端进行测序。所以，理论上说，我们期待reads从转录本的随机位置开始，因此reads中不应该会有碱基组成偏离的情况。但是，实际情况表明，随机六聚体会在RNAseq的reads开头引入碱基组成偏离[7]。这种序列特异的偏差会影响基因以及转录本表达量的估算，因为测序reads对转录本的覆盖是不均一的。在不同碱基位置的严重偏差也可能是未切除接头或者是文库测序过度的信号[1]。 FastQC会给一个reads每个位置的碱基组成图（如下图），理论上应该每种碱基都是一条直线，该值应近似等于物种中该碱基的含量。在任一碱基位置上，如果A和T的差异或者C和G的差异大于10%，则FastQC会显示警告[1]。
序列特异的偏离无法通过过滤或者切除来去除，不过 Cuffinks 和 eXpress包提供了一个校正方法，这个方法学习从数据中挑选出来的序列并把这些信息包含到丰度估计中[1]。 如果只是前面几个碱基出现碱基成分偏差，也可以直接切除。 注意，除了造成序列特异性偏差，随机六聚体引物还会导致Illumina RNAseq reads开头的错配[8]。错配率最高的是第一个核苷酸，虽然有高达7个核苷酸会被影响。错配的碱基通常有很好的质量值，因此不同通过质量来检测它。为了避免这种情况，第一个碱基应当用Trimmomatic 或 PRINSEQ给切除。 3.2.6 GC 含量reads的GC含量通常服从正态分布，中间的GC值应该近似为物种的GC含量。异常的分布性状或跟源基因组的GC含量有很大偏差，可能意味着文库污染。 FastQC和PRINSEQ都会将reads的平均GC含量的频率分布给出来。
注意，标准的文库构建流程会使用PCR扩增，而PCR对于GC含量过高或过低的区域扩增比较难。这些GC含量是样本特异的，从而导致差异表达分析的复杂度增加。GC偏差无法用reads预处理的方式消除。不过，在后续分析阶段，各种各样的标准化和校正方法会考虑这一点[9,10]。 3.2.7 重复序列（Duplicates）duplicates有两种可能：
duplicates检测方法：
是否该去除duplicates？
duplicates的检测和去除FastQC 和 PRINSEQ 都能分析duplicates。FastQC只检测精确的duplicates，而PRINSEQ还会检测5’和3’端的duplicates。
这些工具用于原始reads（raw reads），基于序列识别duplicates。注意，实际中，PCR同样片段产生的duplicates的测序reads未必会完全一致，因为有测序错误。所以，基于序列的方法会低估了duplicates的数目。 当reads被比对到参考基因组，duplicates可以通过mapping位置的一致性（而非序列内容）来识别。 使用PRINSEQ去除原始reads中的duplicates。下边的命令是去除出现超过100次(-derep _min 101) prinseq-lite.pl -fastq reads1.fastq -fastq2 reads2.fastq -derep 1 -derep_min 101 -log -verbose -out_good dupfiltered -out_bad null -no_qual_headerFastQC的duplicate图怎么看？[2,11]看懂这个图其实有点难。首先，红色圈里的值很重要，她表示有多少比例的reads是不一样的。
有两种duplicate的计算方式：
大神Istvan Albert在Tutorial: Revisiting the FastQC read duplication report（https://www./p/107402/）对这作了更详细的讲解： 首先，定义什么是distinct sequences： distinct sequences = singletons（重复值为1的序列）的数目 + doubles（重复值为2的序列）的数目 + triplets（重复值为3的序列）的数目+ … 这里，每种序列不论重复了几次，都只算一种distinct sequences。 图标题上的值 = distinct/总序列数 * 100% 这就是蓝线和红线想要展示的：
举例来说，有两种情况，都含有20条reads。 A B C D E F G H I J K K L L M M N N O OCase 2: 10 unique reads + 1种重复了10次的read A B C D E F G H I J K K K K K K K K K KCase 1 （在上面的图）展示了15种distinct reads，即比例为15/20=75%对于all sequences的计算方法（蓝线）：总数为20的总reads数中：
所以蓝线是一个平台，因为二者比例都为50%。 对于distinct sequences的计算方法（红线）：对于总共为15种distinct sequences有：
所以红线是一条斜线。 Case 2 （在下面的图）则有11种distinct reads，所以比例为11/20=55%。对于all sequences的计算方法（蓝线）：总数为20的总reads数中：
所以，蓝线在重复值为1和10的地方，比例都为50%。 对于distinct sequences的计算方法（红线）：对于总共为11种distinct sequences有：
所以红线在重复值为1时，比例为91%，而在重复值为10时，比例为9%。
3.2.8 序列污染（Sequence Contamination）有的时候reads可能会被别的生物体或者载体的序列污染。这会在GC含量中体现。PRINSEQ的dinucleotide frequency 图也可以给出一些样品污染的线索。但也许最直接的方法就是把reads比对到可能的污染源序列上。FastQ Screen工具可以把reads mapping到用户怀疑的污染源序列上，用的Bowtie比对工具，并且会把结果总结成文本和图形的形式[12]。或者，你可以随机抽取一部分reads，使用BLAST对常规核苷酸数据库进行比对[1]。 3.2.9 低复杂度序列和polyA尾巴（Low-Complexity Sequences and PolyA Tails）低复杂度的序列低复杂度的序列由于缺乏信息量，所以很难可靠地map到参考序列上。比如，单一（homopolymer）重复序列 比如：AAAAAAAAAA，二核苷酸（dinucleotide）重复序列如： CACACACACA，三核苷酸（ trinucleotide）重复序列如：CATCATCATCAT。 PRINSEQ通过两种方法计算reads的序列复杂度：DUST 和 Entropy。
你可以通过选项 PolyA/T尾巴PolyA/T尾巴是reads末端的A或T重复序列。 PRINSEQ会报告有多少reads含有大于等于5个的碱基长度的polyA/T尾巴，以及尾巴长度的频率分布。你可以通过设置最小尾巴长度选项 4 Reads过滤和切除软件推荐Reads过滤和切除软件那么多，应该如何选择？每个人都有自己的偏好。这里笔者选了几个好友以及书籍《Biostar Handbook》和《RNAseq Data Analysis: A Practical Approach》的推荐罗列如下供大家参考： 4.1 健明兄
4.2 天地本无心三种软件组合使用：
4.3 Viswanathan Satheesh
4.4 Biostar Handbook (January 2017)[2]在Biostar Handbook里罗列了许多QC软件，如下：
用于QC的R（Bioconductor）包: PIQA、ShortRead 其中，作者推荐使用以下几种：
4.5 RNAseq Data Analysis: A Practical Approach (June 2014)[1]在这本书里，作者主要介绍了如下几个软件：
最后质控软件太多，这里不可能一一予以介绍。以后如果有时间和经历可能会把软件安装和使用继续写写。 注：本文大部分内容参考了《RNAseq Data Analysis: A Pratical Approach》一书。参考资料
|
|
|
