生信学习笔记：生物信息学测序分析基本流程入门笔记

文章目录

• • • 前言
    • 短序列比对软件
    • sam文件
    • insertsize
    • 基因差异表达计算
    • 变异检测
    • 物种组成与丰富度计算
    • kmer估计基因组大小
    • 序列拼接
    • Pregraph
    • 常用序列拼接软件
    • 基因组污染分析
    • RNA-seq与meta序列拼接
    • 基因功能注释
    • 非编码RNA
    • 小RNA
    • 共线性分析
    • 在线序列分析
    • 序列比对
    • 数据下载

前言

根据B站视频教程生物信息快速入门边自学边随手记的笔记，省略了开头几节测序原理以及数据质控，在之前的文章中有粗略提过该部分。

目前我自己理解的大致流程基本就是

测序——质控——比对拼接（contig、scaffol、mappingd）——注释、差异分析和富集分析等

短序列比对软件

bwa、华大的soap、bwotie2

sam文件

主要用于测序序列mapping到基因组上的结果表示，当然也可以表示任意的多重比对结果

基本是短序列比对默认的标准格式

12列 tab分割

1：序列名字 reads的id

2：序列标记

-1序列是一对序列中的一个

-2比对结果是一个pair-end比对的未端

-4没有找到位点

-8这个序列是pair中的一个但是没有找到位点

-16在这个比对上的位点，序列与参考序列反向互补

-32这个序列在pair-end中的的mate序列与参考序列反向互补

-64序列是mate1

-128序列是mate2

3：参考序列的名字

4：在参考序列上的位置

5：mapping的质量值，越高表示越独特

6：CIGAR字符串

如37M1D2M1137个匹配，1个参考序列上的删除，2个匹配，1个参考序列上的插入

M代表的是alignment match（可以是错配）I代表Insert D代表Deletion

7：mate序列在参考序列上的名称

8：mate序列在参考序列上的位置

9：估计的片段长度，当mate序列位于本序列上游时，该值为负值

10：reads的序列

11：ASCII码格式的reads的质量值

12：12个用空格分隔的可选区域

-AS:i匹配的得分

-XS:i第二好的匹配的得分

-YS:i mate 序列匹配的得分

-XN:i在参考序列上模糊碱基的个数

-XM:i错配的个数

-XO:i gap open的个数

-XG:i gap延伸的个数

-NM:i经过编辑的序列

-YF:i说明为什么这个序列被过滤的字符串

-MD:Z代表序列和参考序列错配的字符串

根据1、10、11即reads部分，看这些能否匹配到目标序列上

insertsize

插入片段insertsize大小，也就是文库片段的长度，同样也是两条测序reads的物理距离

基因差异表达计算

有reads覆盖，则证明该基因有表达，对于一个基因测到的reads越多，则说明该基因表达量越大

基因表达量的计算使用RPKM法（Reads Per Kb per Million reads）

（Mortazavi，2008），其计算公式为：

RPKM=-106C / (NL/103)

设RPKM（A）为基因A的表达量，则C为唯一比对到基因A的reads数，N为唯一比对到参考基因的总reads数，L为基因A编码区的碱基数。RPKM法能消除基因长度和测序量差异对计算基因表达的影响，计算得到的基因表达量可直接用于比较不同样品间的基因表达差异。

如果一个基因存在多个转录本，则用该基因的最长转录本计算其测序覆盖度和表达量。

不适用于发生可变剪切的RNA-seq数据

FPKM

FPKM表示，每1百万个map上的reads中map到外显子的每1K个碱基上的reads个数。

FPKM与RPKM计算方法基本一致。公式如下：

FPKM = total exon Fragments / (mapped reads（Millions） * exon length（KB）)

RPKM计算的是reads，FPKM计算的是片段，片段包含的意义更广，可以是pair-end的一个fragment，也可以是一个reads

比较基因的差异表达

一般选取两个标准：

1、FoldChange

同一个基因表达水平的变化倍数，也就是两个基因FPKM或RPKM的差值，差值倍数越大，差异也越大

2、通过FDR校正

FDR=E（V/R）

FDR（False Discovery Rate）control 就是多重检验校正中一种校正p-value的统计学方法。E为期望值。在实际中，FDR是错误发现率的期望。假设挑选了R个差异表达基因，其中S个是真正有差异表达的基因，另外V个是其实没有差异表达的基因，为假阳性结果。

FDR越小，差异表达倍数越大，则表明差异表达越显著

变异检测

一比一完全匹配、单端read匹配、某个碱基不匹配（可能变异位点）

物种组成与丰富度计算

通过测定16s的高变区确定物种类别

reads组成tags，tags组成OUT

OTU：operational taxonomic unit 分类单元

kmer估计基因组大小

kmer：所谓k-mer，即为一段短的DNA片段（将reads切割为固定长度的小片段）。K为一个奇数，k等于几，就称为几mer。比如17bp大小的片段，则称17mer

绘制整个测序中kmer的频数分布，kmer的峰值为测序深度，通过深度来估计基因组的大小

估计基因组大小：G=S/D；基因组长度=碱基总数/平均深度

​ 我们设G为实际基因组大小，k为kmer长度，|为reads长度，nk为所有kmer的个数，nb为所有碱基的个数，dk为kmer期望深度，db为碱基期望深度，其中nk和dk皆可以统计出来。因为Kmer深度符合汇松分布，所以d k即为Kmer深度的平均值，而整个基因组约可产生G个Kmer，这里面忽略边界效应。

​ 最终我们就得到了如下的公式。

​ G=n_k/d_k=n_b/d_b；

​ d_b/d_k=l/(l-k+1);

​ G=d_k*1/(l-k+1)

​ 得到kmer的深度dk就可以推测出基因组的大小

序列拼接

两条路：短序列比对、序列拼接

小片段文库中称piarend，大片段（>1000bp）称为matepair

contig：reads之间通过overlap连接成contig（准确的说是kmer连成contig），即重叠群

scaffold：原意是脚手架，contig连接成的更长的片段

Pregraph

短序列拼接的步骤：

去除低频率的kmer

1、 构图

De bruijin图，kmer之间的overlap关系

图中每一个kmer都叫做一个node（即特异的kmer数目）；

点与点的连接形成edge，即边，越长越好；

边上还有很多断头，无法形成完整的圈，称为tip

两个kmer从同一个起点node出发经过两条不同的路又回到同一个终点node，称为bubble，可能是基因位点杂合造成

2、构建contig

合并nodes、去tips、合并bubble、解repeat（去除基因组重复区）、输出contig

3、构建scaffold（包括map）

mapping：

reads落在同一contig用于计算insertsize

reads落在不同contig用于构件scaffold

4、补洞

由N碱基构成的gap称为“洞”

关于N区域：

1、N碱基数目不完全准确

2、N区域来自基因组复杂区域，不容易拼接出来

补洞方法：

利用从洞两侧设计引物PCR扩增获得序列或者利用pacbio测序长片段补洞；

利用pairend关系reads补洞；

逐步延伸方法

常用序列拼接软件

SOAPdenovo（华大开发）

软件特点：

操作简单，适合短序列拼接，尤其适合illumina数据，构建Scaffold比较优秀。需要配置文件，一次可以配置多个文库的数据，可以加入454等长reads数据辅助补洞，配有数据细错和补洞程序

使用说明：

1、书写配置文件

2、运行程序

3、最终的scafSeq为最终的拼接结果

注意事项：

1、内控控制

内存的大小与输入数据量有关，其次与数据特点有关，错误率越高，内存消耗越大2、在内存允许范围内，可以尝试使用多个大小kmer，选择一个最优的拼接结果

velvet

软件特点：

操作简单，适短序列拼接，龙基适合illumina数据，支持fastq，fasta，sam，bam等多种格式，局部拼接效果好，不善于连接Scaffold。

使用说明：

1、velveth 构建hash表

2、利用velvetg来拼接

注意事项：

1、不能利用mate-pair文库

2、输入数据需要指定文件格式，数据类型等，多种数据类型时比较混乱3、输出结果命名不能自定义

SPAdes（适合小基因组的拼接）

软件特点：

适合多种测序数据的拼接，既可以利用illumina数据，也可以利用lon Torrent数据或者454数据等。同时可以加入sanger，pacbio nanopore等测序数据，非常适合多种数据之间混合拼接。同时有专门的针对二倍体杂合基因组的版本。可以选择最适合的kmer来拼接。

使用说明：

1、直接在命令行导入数据即可运行

注意事项：

1、需要指定数据类型以及文库类型

2、可以做数据纠错

3、最终的contigs.fasta或者scaffold.fasta即为最终的结果

基因组污染分析

污染来源：样品提取

基因组污染的判断：有两个GC%峰

RNA-seq与meta序列拼接

DNA测序与RNAseq比较：

1、DNA一般为全基因组测序，而RNA测序的是转录出来的转录本，都是独立断开的片段；

2、DNA测序一般是均匀测序，基因组上的区域被均匀覆盖，而RNA由于存在表达丰度的差异；

3、DNA全基因组测序中存在很多重复序列，干扰序列的拼接，而转录组中这个问题影响会小一些。

RNA序列拼接软件：

Trinity

软件特点：

主要应用于转录组数据从头拼接。适用于illumina测序数据的转录组拼接。

使用说明：

第一步，Inchworm：将RNA-seq的原始reads数据组装成Unique序列；第二步，Chrysalis：将上一步生成的contigs聚类，然后对每个类构建Bruijn图；最后，Butterfly：处理这些Bruijn图，依据图中reads和成对的reads来寻找路径，从而得到具有可变剪接的全长转录子，同时将旁系同源基因的转录子分开。

注意事项：

1、java 1.72、Bowtie13、不支持压缩格式数据

4、当数据量比较大的时候，trinity运行的时间会很长，同时，内存不够等情况出现的时候有可能程序运行崩溃。最好是分步运行。

1、拼接结果中获取Unigene（选取最长转录本）

2、拼接时要去除duplication

3、表达定量时不能去除duplication

这里可以将数据保存两份，一份用于RNA-seq拼接，一份用于分析表达量

oases

SOAPdenovo-trans

宏基因组拼接：宏基因组（Metagenome），也称微生物环境基因组Microbial Environmental Genome，或元基因组。

如何改善拼接效果：

1、覆盖基因组所有位点；

2、重新调整insertsize；

3、去除insertsize偏差大的pairend reads；

4、调整kmer大小以及软件选项参数阈值；

5、混合拼接

如何选择kmer的大小？

1、reads准确度越高，选择大kmer较好

2、reads错误率高，选择小kmer，reads利用率高

3、基因组本身特点，重复率，测序深度等因素，都会对kmer取值造成一定影响

kmer只能取奇数

基因功能注释

将氨基酸序列与数据库进行blast比对。同源基因功能相似

**U在打分软件中自动被替换成X：**因为U无法被识别

COG

COG数据库，是Clusters of Orthologous Groups of proteins的简称。，蛋白相邻类的聚簇。该数据库是对细菌、藻类和真核生物的21个完整基因组的编码蛋白，根据系统进化关系分类构建而成的。对于预测单个蛋白质的功能和整个新基因组中的蛋白质的功能非常有用。

GO

一个基因注释是对基因产物的描述

有特定的分子功能（molecular function）

涉及到特定的生物过程（biological process）

作用在特定的细胞组分（cellular component）

KEGG

单个基因的作用、相互关系

pathway

非编码RNA

Noncodeing RNA：ncRNA

rnammer预测rRNA

小RNA

miRNA的作用机制

• 抑制或降解

  取决于miRNA与靶mRNA种子区域的互补程度

  种子区域

  通常指miRNA5’端第二位到第八位的核苷酸序列

• 两者完全互补

  降解

• 两者不完全互补

  抑制翻译

共线性分析

展示不同基因间的亲缘关系，大的差异能通过图直观的一眼看出

所谓共线性，顾名思义，表示二者在一条直线上。基因组的共线性分析也类似，主要是用一种线性图的方式来比较两个或者多个基因组是否具有较好的同源性。

软件：mummer+mummerplot 和 lastz 和 BRIG

在线序列分析

blast

GLIMMER 和 PRODIGAL：在线进行基因预测

RNAmmer：核糖体RNA预测

tRNAscan-SE Search Server：tRNA预测

DAVIDBioinformatics Resources：GO功能分析

POLYAH：预测分析转录终止信号

SignalP 4.1 Server：预测信号肽

Phage_Finder：预测噬菌体

RAST（Rapid Annotation using Subsystem Technology）：利用subsystem技术快速注释工具

可以用来给细菌的完成图和草图来做分析，包括基因预测、ncRNA分析、基因功能注释等，并可以利用这些信息构建代谢网络。免费使用。

分析meta数据可用MG-RAST（国内不能用，且需要火狐浏览器）

序列比对

序列两两比对：

1、两个碱基类型完全一致，没有发生变化；

2、碱基发生了变化，这种变化包括替换、插入和删除。

E-value:E值综合考虑了真实匹配与随机匹配的相对概率，序列的长度，数据库大小和序列组分的偏向性等数据而计算出来的

Bootstrap，即自展值，是用来检验计算的进化树分支可信度的。简单地讲就是把序列的位点都重排，重排后的序列再用相同的办法构树，如果原来树的分枝在重排后构的树中也出现了，就给这个分枝打上一分，如果没出现就给0分，这样经过你给定的repetitions次（至少1000次）重排构树打分后，每个分枝就都得出分值，计算机会给你换算成bootstrap values值。

数据下载

MyBioSoftware：相关生信分析软件

Nucleic Acids Research总结的数据库

NCBI

ftp下载：FTP 是File Transfer Protocol（文件传输协议）的英文简称，用于Internet上的控制文件的双向传输。同时，它也是一个应用程序（Application）。

NCBI中批量下载：Batch Entrez（将需要下载的序列号存入文本文件上传）

Linux下一直想保持最新数据：数据同步rsync