天天看點
傳回

轉錄組入門(6)： reads計數

要求

實作這個功能的軟體也很多，還是煩請大家先自己搜尋幾個教程，入門請統一用htseq-count，對每個樣本都會輸出一個表達量檔案。

需要用腳本合并所有的樣本為表達矩陣。參考：生信程式設計直播第四題：多個同樣的行列式檔案合并起來

對這個表達矩陣可以自己簡單在excel或者R裡面摸索，求平均值，方差。

看看一些生物學意義特殊的基因表現如何，比如GAPDH,β-ACTIN等等。

理論基礎

在上篇的比對中，我們需要糾結是否真的需要比對，如果你隻需要知道已知基因的表達情況，那麼可以選擇alignment-free工具（例如salmon, sailfish)，如果你需要找到noval isoforms，那麼就需要alignment-based工具（如HISAT2, STAR）。到了這一篇的基因（轉錄本）定量，需要考慮的因素就更加多了，以至于我不知道如何說清才能理清邏輯。

定量分為三個水準

  • 基因水準(gene-level)
  • 轉錄本水準(transcript-level)
  • 外顯子使用水準(exon-usage-level)。

在基因水準上，常用的軟體為HTSeq-count，featureCounts，BEDTools, Qualimap, Rsubread, GenomicRanges等。以常用的HTSeq-count為例，這些工具要解決的問題就是根據read和基因位置的overlap判斷這個read到底是誰家的孩子。值得注意的是不同工具對multimapping reads處理方式也是不同的，例如HTSeq-count就直接當它們不存在。而Qualimpa則是一人一份，平均配置設定。

轉錄組入門(6)： reads計數

_images/count_modes.png

對每個基因計數之後得到的count matrix再後續的分析中，要注意标準化的問題。如果你要比較同一個樣本(within-sample)不同基因之間的表達情況，你就需要考慮到轉錄本長度，因為轉錄本越長，那麼檢測的片段也會更多，直接比較等于讓小孩和大人進行賽跑。如果你是比較不同樣本（across sample）同一個基因的表達情況，雖然不必在意轉錄本長度，但是你要考慮到測序深度（sequence depth)，畢竟測序深度越高，檢測到的機率越大。除了這兩個因素外，你還需要考慮GC%所導緻的偏差，以及測序儀器的系統偏差。目前對read count标準化的算法有RPKM（SE）, FPKM（PE），TPM, TMM等，不同算法之間的差異與換算方法已經有文章進行整理和吐槽了。但是，有一些下遊分析的軟體會要求是輸入的count matrix是原始資料，未經标準化，比如說DESeq2，這個時候你需要注意你上一步所用軟體會不會進行标準化。

在轉錄本水準上，一般常用工具為Cufflinks和它的繼任者StringTie， eXpress。這些軟體要處理的難題就時轉錄本亞型（isoforms）之間通常是有重疊的，當二代測序讀長低于轉錄本長度時，如何進行區分？這些工具大多采用的都是expectation maximization（EM）。好在我們有三代測序。

上述軟體都是alignment-based，目前許多alignment-free軟體，如kallisto, silfish, salmon，能夠省去比對這一步，直接得到read count，在運作效率上更高。不過最近一篇文獻[1]指出這類方法在估計豐度時存在樣本特異性和讀長偏差。

在外顯子使用水準上，其實和基因水準的統計類似。但是值得注意的是為了更好的計數，我們需要提供無重疊的外顯子區域的gtf檔案[2]。用于分析差異外顯子使用的DEXSeq提供了一個Python腳本（dexseq_prepare_annotation.py）執行這個任務。

小結

計數分為三個水準： gene-level, transcript-level, exon-usage-level

标準化方法： FPKM RPKM TMM TPM

輸出表達矩陣

在RNA-Seq分析中，每一個基因就是一個feature（特征？），而基因被認為是它的所有外顯子的和集。在可變剪切分析中，可以單獨把每個外顯子當作一個feature。而在ChIP-Seq分析中，feature則是預先定義的結合域。但是确定一個read到底屬于哪一個feature有時會非常棘手。是以HTSeq提供了三種模式，示意圖見前一幅圖

  • the union of all the sets S(i) for mode union. This mode is recommended for most use cases.
  • the intersection of all the sets S(i) for mode intersection-strict.
  • the intersection of all non-empty sets S(i) for mode intersection-nonempty.

基本用法非常的簡單：

# 安裝
conda install htseq
# 使用
# htseq-count [options] <alignment_file> <gtf_file>
htseq-count -r pos -f bam RNA-Seq/aligned/SRR3589957_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > SRR3589957.count

用一個循環處理多個BAM檔案(在/mnt/f/Data目錄下)

mkdir -p RNA-Seq/matrix/
for i in `seq 56 58`
do
    htseq-count -s no -r pos -f bam RNA-Seq/aligned/SRR35899${i}_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > RNA-Seq/matrix/SRR35899${i}.count 2> RNA-Seq/matrix/SRR35899${i}.log
done

運作的時間會比較久，是以可以去了解不同參數的用法了，其中比較常用的為：

  • -f bam/sam： 指定輸入檔案格式，預設SAM
  • -r name/pos: 你需要利用samtool sort對資料根據read name或者位置進行排序，預設是name
  • -s yes/no/reverse: 資料是否來自于strand-specific assay。DNA是雙鍊的，是以需要判斷到底來自于哪條鍊。如果選擇了no， 那麼每一條read都會跟正義鍊和反義鍊進行比較。預設的yes對于雙端測序表示第一個read都在同一個鍊上，第二個read則在另一條鍊上。
  • -a 最低品質， 剔除低于門檻值的read
  • -m 模式 union（預設）, intersection-strict and intersection-nonempty。一般而言就用預設的，作者也是這樣認為的。
  • -i id attribute: 在GTF檔案的最後一欄裡，會有這個基因的多個命名方式（如下）， RNA-Seq資料分析常用的是gene_id， 當然你可以寫一個腳本替換成其他命名方式。
gene_id "ENSG00000223972.5_2"; transcript_id "ENST00000456328.2_1"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-002"; exon_number 2; exon_id "ENSE00003582793.1_1"; level 2;...

Jimmy的伏筆

我們這次分析是人類mRNA-Seq測序的結果，但是我們其實隻下載下傳了3個sra檔案。一般而言RNA-Seq資料分析都至少要有2個重複，是以必須要有4個sra檔案才行。我在仔細讀完文章的方法這一段以後，發現他們有一批資料用的是其他課題組的： For 293 cells, the mRNA-seq results of the control samples include (1) those from the doxycycline-treated parental Flp-In T-REx 293 cells by us and (2) those from the doxycycline-treated control Flp-In T-REx 293 cells performed by another group unrelated to us (sample GSM1095127 in GSE44976)。 然後和Jimmy交流之後，他也承認自己隻分析了小鼠的資料，而沒有分析人類的資料。是以我們需要根據文章提供的線索下載下傳另外一份資料，才能進行下一步的分析。

這個時候就有一個經常被問到的問題：不同來源的RNA-Seq資料能夠直接比較嗎？甚至說如果不同來源的RNA-seq資料的建構文庫都不一樣該如何比較?不同來源的RNA-Seq結果之間比較需要考慮 批次效應（batch effect) 的影響。

處理批次效應，根據我搜尋的結果，是不能使用FPKM/RPKM，關于這個标準化的吐槽，我在biostars上找到了如下觀點：

  • FPKM/RPKM 不是标準化的方法，它會引入文庫特異的協變量
  • FPKM/RPKM has never been peer-reviewed, it has been introduced as an ad-hoc measure in a supplementary 沒有同行評審
  • One of the authors of this paper states, that it should not be used because of faulty arithmetic 作者說算法有問題
  • All reviews so far have shown it to be an inferior scale for DE analysis of genes Length normalization is mostly dispensable imo in DE analysis because gene length is constant

有人建議使用一個Bioconductor包

http://www.bioconductor.org/packages/devel/bioc/html/sva.html

我沒有具體了解，有生之年去了解補充。

還有人引用了一篇文獻 IVT-seq reveals extreme bias in RNA-sequencing 證明不同文庫的RNA-Seq結果會存在很大差異。

結論： 可以問下原作者他們是如何處理資料的，居然有一個居然沒有重複的分析也能過審。改用小鼠資料進行分析。或者使用無重複的分析方法，或者模拟一份資料出來，先把流程走完。

合并表達矩陣

HTSeq-count輸出結果是一個個獨立的檔案，後續分析需要把多個檔案合并成一個行為基因名，列為樣本名，中間為count的行列式檔案。肯定是不會用Excel手動處理（雖然可以寫一個VBA腳本，但是資料量過大不好處理了），這裡使用的Python寫一個腳本。

基本邏輯：

  1. 讀取檔案
  2. 建立一個字典，如果key不在字典中，新增key和value，如果key在字典中，追加value。
  3. 輸出 
#!/usr/bin/python3

import sys
mydict = {}
for file in sys.argv[1:]:
    for line in open(file, 'r'):
        key,value = line.strip().split('\t')
        if key in mydict:
            mydict[key] = mydict[key] + '\t' + value
        else:
            mydict[key] = value

for key,value in mydict.items():
    print(key + '\t' + value +'\n')

這幾行代碼寫了2個番茄鐘，但是debug花了我一個番茄鐘。問題出在str和list兩種資料格式的混亂使用。還有一個bug： 由于詞典是無序的，是以原本代表樣本來源的第一行，會跑到其他行。

在論壇上找到一個更加簡潔的代碼（要求基因名順序排列），用到

paste

 ,

awk

printf

 這三個shell指令。 
paste *.txt | awk '{printf $1 "\t";for(i=2;i<=NF;i+=2) printf $i"\t";printf $i}'

儲存為countCombiner.py，輸入檔案為count, 輸出為标準輸出，需要重定向。

簡單分析

這一步需要用到R語言或者是Excel讀取資料。

1.導入資料

options(stringsAsFactors = FALSE)
# import data if sample are small
control <- read.table("F:/Data/RNA-Seq/matrix/SRR3589956.count",
                       sep="\t", col.names = c("gene_id","control"))
rep1 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589957.count",
                    sep="\t", col.names = c("gene_id","rep1"))
rep2 <- read.table("F:/Data/RNA-Seq/matrix/SRR3589958.count",
                    sep="\t",col.names = c("gene_id","rep2"))

2.資料整合。gencode的注釋檔案中的gene_id（如ENSG00000105298.13_3）在EBI是不能搜尋到的，是以我就隻保留ENSG00000105298這部分。

# merge data and delete the unuseful row
raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")
raw_count_filt <- raw_count[-1:-5,]
ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)
row.names(raw_count_filt) <- ENSEMBL

3.總體情況, 大部分基因都為0，是以可以删掉節省體積。

summary(raw_count_filt)
    control              rep1               rep2         
 Min.   :     0.0   Min.   :     0.0   Min.   :     0.0  
 1st Qu.:     0.0   1st Qu.:     0.0   1st Qu.:     0.0  
 Median :     0.0   Median :     0.0   Median :     0.0  
 Mean   :   356.1   Mean   :   370.3   Mean   :   316.6  
 3rd Qu.:    15.0   3rd Qu.:    15.0   3rd Qu.:    10.0  
 Max.   :161867.0   Max.   :121902.0   Max.   :105565.0

4.看幾個具體基因

在EBI上搜尋GAPDH找到ID為ENSG00000111640。

GAPDH <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000111640",]
                             gene_id control  rep1  rep2
ENSG00000111640 ENSG00000111640.14_2   41857 53902 55302

文章研究的AKAP95（ENSG00000105127）的表達量在KD中都降低了

> AKAP95 <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000105127",]
> AKAP95
                            gene_id control rep1 rep2
ENSG00000105127 ENSG00000105127.8_2    1168  539  506

下面的差異基因表達，讓我想想，該如何收拾Jimmy挖的坑。

參考文獻

[1] Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis

[2] Detecting differential usage of exons from RNA-seq data

[3] RNA-seq Data Analysis-A Practical Approach(2015)

算法 Python python算法 python算法問題 python算法選擇 python算法通過 算法python3
上一篇: Python和R的異同(一）R的向量矩陣數組與Python的清單
下一篇: R語言的循環

繼續閱讀