生信學習筆記：生物資訊學測序分析基本流程入門筆記

文章目錄

• • • 前言
    • 短序列比對軟體
    • 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