NGS基础

FASTQ格式解释和质量评估

• FASTQ文件格式和命名
  高通量测序之后用于下游分析的数据一般存储在FASTQ文件中。为了节省空间，又不影响下游使用，也一般用gzip压缩的格式。

单端测序每个文库只返回一个FASTQ文件，
双端测序两个FASTQ文件，左端一般命名为_1或R1，右端命名为_2或R2。

image.png

image.png

第一行以@开头，后面是reads的ID以及其他信息。

第二行为read的序列。

第三行以+开头，一般后面没有内容；若有则为序列的名字，与第一行相同。

第四行代表reads的质量值。质量值的计算方式为 (e是错误率，estimated probability of the base call being wrong)。如果该碱基测序出错的概率(e)为0.001，则Q应该为30。现在测序数据多采用Phred33编码，那么30+33=63，那么63对应的ASCii码为?，则在第四行中该碱基对应的质量代表值即为?。一般地，碱基质量从0-40, 即ASCii码为从 !(0+33)到I(40+33)。

image.png

image.png

image.png

• FASTQ文件一般比较大，传输过程中可能会出现不完整情况，一般都会附带一个MD5校验值文件用于判断文件的完整性。如果文件内容不变，则MD5值不变，与是否压缩无关。

#获取MD5
md5sum test.fq.gz | tee >MD5.txt
#验证MD5
md5sum -c MD5.txt

测序量评估

• 测序深度：
  测序得到的碱基总量（bp）/参考基因组大小
  （基因组中每个碱基被测序到的平均次数）

• 测序覆盖度（coverage）
  测序获得的序列占整个基因组的比例

  
  
  image.png

• 单端测序
  数据量=reads长度 * reads个数 (reads长度很容易得知，reads个数等于测序所得到的fastq文件的总reads数)

• 双端测序
  数据量=单端reads长度 * 单端reads个数 * 2
  通常测序数据量的单位都是用“G"表示，例如1G

1个碱基=1 byte ；
1kb=10^3 byte 1M=10^6 byte 1G=10^9 byte。
：1G的数据量=10^9=10亿个碱基

测序reads数和测序碱基数。

zcat test.fq.gz
# 获取行数，然后除以4，得到reads数
 echo $(( `zcat test.fq.gz | wc -l` / 4 ))

## 获取碱基数
 zcat test.fq.gz | grep -A 1 '^@' | grep '^[ACGTN]' | wc -c

测序质量

FastQC软件

fastqc 1.fq.gz 2.fq.gz
FastQC质控报告解读

image.png

image.png

高通量测序原理

NGS基础 - 高通量测序原理 (qq.com)
涉及测序文库构建原理、链特异性文库的构建方式和识别方法、测序簇生成过程、双端测序过程、测序接头产生、PCR duplicate、测序通量选择标准等。

文库构建

文库构建即为测序片段添加接头即为测序片段添加接头

常用的试剂盒NEBNext®Ultra™ II DNA Library Prep Kit for Illumina®为例阐述二代测序文库构建的流程及其原理，具体如下所示：

从基因组提取DNA/RNA（RNA须逆转录为cDNA），打断（超声、酶切）。得到小片段后

1.末端修饰。（补平末端，加3'A，5'磷酸化）
目前很多PCR使用的高保真Pfu聚合酶产生的片段末端是平齐的（也即没有不配对的碱基）；鸟枪法产生的片段则是随机断裂，其末端可能是平齐的也可能是不平的。
因此，建库第一步是使用Taq聚合酶补齐不平的末端，并在两个末端添加突出的碱基A，从而产生粘性末端（若使用Taq酶扩增，则无需末端修饰），产生粘性末端的片段可以添加接头（Adaptor）。

2.添加接头。
接头具有突出的T尾，可以使用连接酶将接头添加到DNA片段两端
接头：
Rd1 SP：PE 测序引物
Rd2 SP：SE 测序引物
P5\P7:Flow cell 连接位点

3.磁珠纯化。
添加接头后的文库体系中含有聚合酶、连接酶等各种酶以及辅助物质，接头的添加也是过量的，而且由于末端的不稳定性，容易形成自连片段，鸟枪法打断的片段中也可能有大片段存在，所以需要特殊磁珠（AMPure XP Beads）纯化来去除大片段以及各种杂质，从而获得成功添加接头的文库片段。

• 原理为磁珠可以通过氢键等作用力来吸附DNA片段，磁珠本身不具有片段大小选择的能力，但其储存的buffer里面含有20%的PEG 8000，PEG浓度越大则可以吸附的DNA片段越小。因此磁珠纯化的时候要根据文库片段不同严格控制磁珠添加量（其实是PEG添加量）来实现片段选择。

4.PCR扩增。
使用与接头互补的引物来扩增。片段还需要添加用于区分不同文库的特异性index，以及与测序仪芯片互补的两种寡核苷酸序列（P5/P7）。

测序是以单链为单位的，建库完成后的每条DNA的单链均一端连有测序引物Read1 Sequencing Primer（Rd1SP）和P5，另一端为Rd2 SP、Index（Barcode）和P7。Index用来区分不同的文库，因为测序仪一个run产生数据量巨大，由于实际情况不同，一次上机常会进行多个文库测序，因此需要加上Index来区分。

在建库过程中，文库中每个DNA短片段的正链与反链都加上了P5与P7，因此建库后每个DNA片段都会扩增出两种结果（详见上面插图），如果全部上机，最终两条链都会有测序结果。因为上机测序起始是以DNA单链为单位，单链化的DNA片段进入测序仪流通池，会随机的结合在不同位置，且相互距离足够远以保证测序信号的独立读取。最终获得的测序结果会有重复的reads(反向互补也会有重复)，所以都会有去重步骤，而且测序量越大重复率会越高。

5.第二次磁珠纯化。
PCR后需要将产物DNA片段与聚合酶等杂质分离，因此再次进行磁珠纯化，之后进行质量检测，包括DNA浓度检测、琼脂糖凝胶电泳和片段长度检测，完成建库

上机测序

Illumina测序技术为基于基因芯片的边合成边测序，整个平台可解剖为三个系统：
1 温度控制系统，原理和普通PCR仪一样，来控制反应的进行；
2 酶控制系统，通过各种酶来控制DNA合成与剪切；
3 荧光信号收集系统

在Illumina测序平台的流通池（Flow cell）表面，通过基因芯片技术交错固定了无数条寡核苷酸链（即短核苷酸链），分别为P5'(P5互补)和P7，单链化的文库DNA片段进入流通池后，包含P5或P7'的单链可以与表面的寡核苷酸基于互补配对结合，从而进入测序过程。测序具体流程如下：

1. 首先以寡核苷酸为引物、文库片段为模板进行DNA复制（因为文库稀释后浓度足够低，可以认为文库片段均匀的结合在流通池表面，每个片段结合的位置相距足够远，这很重要，否则测序时会导致信号叠加而不能识别）。复制完成后解链，将文库片段洗去，留在流通池表面的为与文库模板互补的DNA链。

2. 因为单链DNA另一端为不同的接头序列，可以与相邻的另一种寡核苷酸互补结合，之后进行“桥”式扩增（假如第一次结合的为P7，则复制完成洗脱模板后顶端可以与相邻的P5互补结合形成“桥”，并以P5为引物进行复制，完成后再次解链并与相邻不同种接头结合来进行复制，如此类推）。25-28个循环完成后，原来散布在表面的单核苷酸序列变成散布的DNA簇，这一步主要是为后续测序做准备，因为测序时单分子产生的光信号很弱，难以检测。
   

   image.png

3.“桥”式扩增后一个DNA簇都是由最初的一个文库模板复制而来，但是这时候P7上的序列与P5上的序列是分别从两端开始的，测序要保证每个片段一致性（都是正向或都是反向），因此再次解链线性化，切割并洗去P5上的DNA链，只留P7上的DNA单链。Illumina巧妙地利用了甲酰胺基嘧啶糖苷酶Fpg对8-氧鸟嘌呤糖苷8-oxo-G的选择性切断作用，在合成的引物链上加入了一个8-oxo-G，用Fpg处理，就把带8-oxo-G基团切掉，并把DNA链切断，留下一带不完整糖基的磷酸基。这个磷酸基在接下来的过程中，起到了阻止P5延伸的作用。此后的双末端测序中需要恢复3'-OH，则用脱嘌呤嘧啶内切核酸酶AP-endonuclease把带不完整糖基的那个磷酸基切掉。

4.加入测序引物Read1 SP和修饰过的DNA聚合酶，则在测序引物3’端开始DNA复制。在流通池加入可逆终止荧光dNTP，其3'-OH被阻隔（糖基3'连接有叠氮基团，在链延伸时起到了阻止添加下一个dNTP作用，因此在除去阻隔前只能添加一个碱基），4种dNTP在碱基上分别连接有不同颜色的荧光基团（也可以相同颜色荧光标记，但是测序会更慢，每次只能添加一种碱基）。之后洗掉多余的dNTP，使用激光扫描，收集留在流通池表面的荧光信号。用巯基试剂去掉3’位阻断的叠氮基团，用TCEP(Tris(2-carboxyethyl)phosphine,三(2-羧乙基)膦)去掉荧光基团，进入下一个碱基的测序反应。因为每条DNA单链扩增形成的DNA簇均固定在表面，随着反应进行根据相同位置出现的荧光信号情况，就逐渐读出了改位点DNA链的序列。

5.要保证测序的准确性，需要一个位点DNA簇的每条链同步复制，然而随着反应进行，不同链复制情况会出现差异，因此二代测序读长目前限制在300bp以内。Read1结束后，解链并洗掉测序中已经合成的部分，加入测序引物Index引物（也即Read2 SP互补的寡核苷酸），这时会继续在3’端进行复制，读出接头中Index序列，从而可以确定出每个位点的DNA属于哪个文库。

⑥为了增长测序长度，进行另一个方向测序，也即双末端测序。洗掉前面复制合成的片段，DNA单链继续在流通池表面形成桥式连接，这时要用脱嘌呤嘧啶内切核酸酶处理修复P5的3’-OH末端，加入聚合酶，则在P5末端开始DNA复制。十几个循环后，将P7上的DNA切割并洗掉。Illumina通过在P7核酸链中加入一个U碱基，用USER酶（Uracil Specific Excision Reagent，尿嘧啶链特定切断试剂）来切隔断链。这时只留下P5上的DNA链，与Read中方向相反。加入测序引物Read2 SP，进行另一端的序列读取。

测序结果

测序数据
一般我们接触到的测序数据为fastq格式的碱基序列，然而早期Illumina平台直接下机数据为bcl格式文件，其储存的是显微拍摄得到的荧光信号信息，

参考基因组和基因注释文件获取

通常测序生成的reads要与参考基因组或参考转录组进行比对，或Pseudo-alignment。所以首先需要获取参考基因组和参考转录组信息

Ensembl

http://www.ensembl.org/info/data/ftp/index.html 是常用的信息齐全的参考基因组和GTF文件下载网站。
Ensembl提供的参考基因组有2种组装形式和3种重复序列处理方式, 分别是primary,toplevel和unmasked(dna)、soft-masked(dna_sm)和masked (dna_rm)。

一般选择dna.primary或dna_sm.primary。

• 为什么选择Primary

Primary assembly contains all toplevel sequence regions excluding haplotypes and patches. This file is best used for performing sequence similarity searcheswhere patch and haplotype sequences would confuse analysis.

• 为什么不选择masked

Masked基因组是指所有重复区和低复杂区被N代替的基因组序列，比对时就不会有reads比对到这些区域。一般不推荐用masked的基因组，因为它造成了信息的丢失，由此带来的一个问题是uniquely比对到masked基因组上的reads实际上可能不是unique的。而且masked基因组还会带来比对错误，使得在允许错配的情况下，本来来自重复区的reads比对到基因组的其它位置。 另外检测重复区和低复杂区的软件不可能是完美的，这就造成遮盖住的重复序列和低复杂区并不一定是100%准确和敏感的。

soft-masked基因组是指把所有重复区和低复杂区的序列用小写字母标出的基因组，由于主要的比对软件，比如BWA、bowtie2等都忽略这些soft-mask，直接把小写字母当做大写字母比对，所以使用soft-masked基因组的比对效果和使用unmasked基因组的比对效果是相同的。

• ENSEMBL的基因注释文件与GeneCode（http://www.gencodegenes.org/）V26版本一致。

ENSEMBL中基因组和GTF文件中染色体的名字都没有添加chr，最好收到添加，以保持与UCSC或下游操作一致。

下载基因功能和结构注释信息

ENSEMBL数据库的BioMart http://www.ensembl.org/biomart/martview工具为下载基因的功能信息、序列信息、结构信息、ID的转换等提供了很大的便利。

注意在BioMart的Attribute选项里如果选择了蛋白相关的选项，得到的结果中只有蛋白编码基因的信息。如果要下载所有基因信息，请不要选择蛋白相关的选项。
示例：Ensembl Genes 89数据集
选择Ensembl Genes 89数据集>
以Human为例，选择Human genes (GRCh38.p10)>如果下载全部的基因信息，Filters部分可以略过不填。如果只想下载比如说某个GO通路的基因或给定列表的基因信息，可以在Filters中指定对应的GO ID。>Attribute中包含基因的名字、位置、注释、在不同数据库中的名字、GO注释、KEGG注释、功能域信息等，按需选择下载。>选择好后，点击Results，获取结果>Export al results to选择存储到文件中。如果特别大，而自己网速又比较慢，可以选择通过邮件发送下载链接

也可以通过Biomart提取基因结构信息，比如5’ UTR、3’ UTR、外显子、内含子的坐标等。

Biomart下载很方便，但一个点击也比较麻烦，可以看到截图中存在XML按钮，点击打开看到选择的下载信息都记录在了这个文件中。

wget -O result.txt 'http://www.ensembl.org/biomart/martservice?query= + XML中的内容 (调整为一行，并且行尾加一个单引号)即可反复使用。如果想换一个物种，只需修改对应的Dataset name

GTF/GFF文件格式解读和转换

NGS基础 - GTF/GFF文件格式解读和转换

GFF 文件GFF全称为general feature format，这种格式主要是用来注释基因组。

从 Ensembl 导出的GFF文件示例：

X    Ensembl    Repeat    2419108    2419128    42    .    .    hid=trf; hstart=1; hend=21
X    Ensembl    Repeat    2419108    2419410    2502    -    .    hid=AluSx; hstart=1; hend=303
X    Ensembl    Repeat    2419108    2419128    0    .    .    hid=dust; hstart=2419108; hend=2419128
X    Ensembl    Pred.trans.    2416676    2418760    450.19    -    2    genscan=GENSCAN00000019335
X    Ensembl    Variation    2413425    2413425    .    +    .    
X    Ensembl    Variation    2413805    2413805    .    +    .

GFF文件是以tab键分割的9列组成，
以下为每一列的对应信息：

• seq_id：序列的编号，一般为chr或者scanfold编号；

• source: 注释的来源，一般为数据库或者注释的机构，如果未知，则用点“.”代替

• type: 注释信息的类型，比如Gene、cDNA、mRNA、CDS等;

• start: 该基因或转录本在参考序列上的起始位置；(从1开始，包含);

• end: 该基因或转录本在参考序列上的终止位置；(从1开始，包含);

• score: 得分，数字，是注释信息可能性的说明，可以是序列相似性比对时的E-values值或者基因预测是的P-values值，.表示为空;

• strand: 该基因或转录本位于参考序列的正链(+)或负链(-)上;

• phase: 仅对注释类型为“CDS”有效，表示起始编码的位置，有效值为0、12. (对于编码蛋白质的CDS来说，本列指定下一个密码子开始的位置。每3个核苷酸翻译一个氨基酸，从0开始，CDS的起始位置，除以3，余数就是这个值，，表示到达下一个密码子需要跳过的碱基个数。该编码区第一个密码子的位置，取值0,1,2。0表示该编码框的第一个密码子第一个碱基位于其5’末端；1表示该编码框的第一个密码子的第一个碱基位于该编码区外；2表示该编码框的第一个密码子的第一、二个碱基位于该编码区外；如果Feature为CDS时，必须指明具体值。)；

• attributes:
  一个包含众多属性的列表，格式为“标签＝值”(tag=value)，以多个键值对组成的注释信息描述，键与值之间用“=”，不同的键值用“；”隔开，一个键可以有多个值，不同值用“,”分割。注意如果描述中包括tab键以及“，= ；”，要用URL转义规则进行转义，如tab键用 代替。键是区分大小写的，以大写字母开头的键是预先定义好的，在后面可能被其他注释信息所调用。

预先定义的键主要包括：

ID：注释信息的编号，在一个GFF文件中必须唯一；

name：注释信息的名称，可以重复；Alias：别名；Parent > >

Indicates：该注释所属的注释，值为注释信息的编号，比如外显子所属的转录组编号，转录组所属的基因的编号。

Parent指明feature所从属的上一级ID，用于将exons聚集成transcript，将transripts聚集成gene，值可以为多个；

Target 指明比对的目标区域，一般用于表明序列的比对结果。格式为 “target_idstart end [strand] ,其中strand是可选的(“+”或”-”),target_id中如果包含空格，则要转换成’ ‘。

Gap：T比对结果的gap信息，和Target一起，用于表明序列的比对结果。Derives_from：Note:备注;Dbxref:数据库索引。

GTF 文件

GTF全称为gene transfer format，主要是用来对基因进行注释。

从 Ensembl 导出的 GTF 文件示例：

image.png

GTF格式大部分与GFF相同，但有两个硬性标准：

feature types是必须注明的;

第9列必须以gene_id以及transcript_id开头。而且GTF文件的第9列同GFF文件不同，虽然同样是标签与值配对的情况，但标签与值之间以空格分开，且每个特征之后都要有分号;(包括最后一个特征);

gene_id “geneA”;transcript_id “geneA.1”;database_id “0012”;modified_by “Damian”;duplicates 0;

每个区域的行首增加chr标签，并去掉#开头的行。

grep -v '^#' GRCh38.gtf | sed 's/^/chr/' >GRCh38.new.gtf

GFF 文件与 GTF 文件相互转换

• 使用Cufflinks里面的工具gffread

#gff2gtf
gffread my.gff3 -T -o my.gtf
#gtf2gff
gffread merged.gtf -o- > merged.gff3

统计GTF文件中染色体数目

cut -f1 GRCh38.new.gtf | uniq -c >chrCount.txt
awk '{print $2"\t"$1}' chrCount.txt 
 awk '/chr[0-9XYM]/' chrCount.txt 
awk '/chr[0-9XYM]/' chrCount.txt | sed 's/  *\([0-9]*\) \(chr.*\)/\2\t\1/' 

统计GTF文件中基因数目

time cut -f 3 GRCh38.new.gtf | sort | uniq -c
time awk '{a[$3]+=1}END{for(i in a) print i,a[i];}' GRCh38.new.gtf

awk '{if(a[$3]=="") a[$3]=1; else a[$3]=a[$3]+1;}END{for(i in a) print i,a[i];}' GRCh38.new.gtf

计算GTF中外显子总长度

# 这个是冗余的外显子
awk '{if($3=="exon") sum+=$5-$4+1;}END\
    {print "Total redundant exon length", sum;}' GRCh38.new.gtf

计算GTF文件中每个基因的转录本数目

# 第一个办法：基因和对应的转录本是排序好的，直接判断计算就可以

awk 'BEGIN{OFS=FS="\t";}{if($3=="gene" && count>0) {print count; count=0;} else \
    {if($3=="transcript") count+=1;}}END{print count}' GRCh38.new.gtf

# 第二个方法：取出所有基因和转录本名字
sed 's/"/\t/g' GRCh38.new.gtf | awk '$3=="transcript"' | cut -f 10,14 | cut -f 1 | uniq -c

# 第三个方法：与第二个类似，但使用了groupBy

sed 's/"/\t/g' GRCh38.new.gtf | awk '$3=="transcript"' | cut -f 10,14 | \
    bedtools groupby -g 1 -c 1,2 -o count,collapse | head

sed 's/"/\t/g' GRCh38.new.gtf | awk '$3=="transcript"' | cut -f 10,14 | \
    bedtools groupby -g 1 -c 1,2 -o count,collapse >geneTrCount.txt

计算GTF文件中基因所拥有的平均转录本数目

awk 'BEGIN{OFS=FS="\t"}{sum+=$2}END{print sum/NR;}' geneTrCount.txt

GTF 文件中提取转录本序列(.fa)

• Cufflink中的gffread

gffread transcripts.gtf –g genome.fa –w transcripts.output.fa
# 获取CDS序列
gffread transcripts.gtf –g genome.fa -x cds.output.fa
# 获取蛋白序列
gffread transcripts.gtf –g genome.fa -y protein.output.fa

• Tophat中的gtf_to_fasta

gtf_to_fasta transcripts.gtf genome.fa out_file

获取启动子区序列

image.png

test.fa中的序列全转成大写

sed -i '/^[^>]/ s/.*/\U&/' test.fa

计算多行FASTA文件test.fa中每条序列长度

#输出类似genome.txt格式的文件(文件有两列，第一列为序列ID，第二列为序列长度)
# 计算一个输出一个
awk 'BEGIN{OFS="\t"; size=0;}{if($0~/>/) {if(size>0) print geneName,size; \
    geneName=$0; sub(">","",geneName); size=0;} else \
    {size+=length}}END{print geneName,size}' test.fa

# 全部计算完存储起来再输出
awk 'BEGIN{OFS="\t";}{if($0~/>/) {geneName=$0; sub(">","",geneName); size[geneName]=0;} \
    else {size[geneName]+=length($0)}}END\
    {for (geneName in size) print geneName,size[geneName]}' test.fa

多行FASTA转单行FASTA序列

# conditions?true_value:false_value 三目运算符，条件为真时，返回冒号前结果，否则冒号后结果
# 对于非第一行的>，输出前先输出一个换行
awk '/^>/&&NR>1{print "";}{printf "%s",/^>/?$0"\n":$0}' test.fa >singleLine.fa

取出单行FASTA文件中序列长度大于40的序列的名字

awk 'BEGIN{OFS="\t";}{if($0~/>/) {geneName=$0; sub(">","",geneName); } else \
    {if (length($0)>40) print geneName;}}' singleLine.fa

定一个BAM文件，怎么计算有多少基因组区域被测到了？平均测序深度是多少？

bedtools genomecov -ibam ../bio/map.sortP.bam -bga