原文鏈接：HOPTOP轉(zhuǎn)錄組入門(6)： reads計數(shù)-轉(zhuǎn)錄組-生信技能樹
要求
實現(xiàn)這個功能的軟件也很多，還是煩請大家先自己搜索幾個教程，入門請統(tǒng)一用htseq-count，對每個樣本都會輸出一個表達量文件。
需要用腳本合并所有的樣本為表達矩陣。參考：生信編程直播第四題：多個同樣的行列式文件合并起來
對這個表達矩陣可以自己簡單在excel或者R里面摸索，求平均值，方差。
看看一些生物學意義特殊的基因表現(xiàn)如何，比如GAPDH,β-ACTIN等等。
理論基礎(chǔ)
在上篇的比對中，我們需要糾結(jié)是否真的需要比對，如果你只需要知道已知基因的表達情況，那么可以選擇alignment-free工具（例如salmon, sailfish)，如果你需要找到noval isoforms，那么就需要alignment-based工具（如HISAT2, STAR）。到了這一篇的基因（轉(zhuǎn)錄本）定量，需要考慮的因素就更加多了，以至于我不知道如何說清才能理清邏輯。
定量分為三個水平
基因水平(gene-level)
轉(zhuǎn)錄本水平(transcript-level)
外顯子使用水平(exon-usage-level)。
在基因水平上，常用的軟件為HTSeq-count，featureCounts，BEDTools, Qualimap, Rsubread, GenomicRanges等。以常用的HTSeq-count為例，這些工具要解決的問題就是根據(jù)read和基因位置的overlap判斷這個read到底是誰家的孩子。值得注意的是不同工具對multimapping reads處理方式也是不同的，例如HTSeq-count就直接當它們不存在。而Qualimpa則是一人一份，平均分配。
對每個基因計數(shù)之后得到的count matrix再后續(xù)的分析中，要注意標準化的問題。如果你要比較同一個樣本(within-sample)不同基因之間的表達情況，你就需要考慮到轉(zhuǎn)錄本長度，因為轉(zhuǎn)錄本越長，那么檢測的片段也會更多，直接比較等于讓小孩和大人進行賽跑。如果你是比較不同樣本（across sample）同一個基因的表達情況，雖然不必在意轉(zhuǎn)錄本長度，但是你要考慮到測序深度（sequence depth)，畢竟測序深度越高，檢測到的概率越大。除了這兩個因素外，你還需要考慮GC%所導致的偏差，以及測序儀器的系統(tǒng)偏差。目前對read count標準化的算法有RPKM（SE）, FPKM（PE），TPM, TMM等，不同算法之間的差異與換算方法已經(jīng)有文章進行整理和吐槽了。但是，有一些下游分析的軟件會要求是輸入的count matrix是原始數(shù)據(jù)，未經(jīng)標準化，比如說DESeq2，這個時候你需要注意你上一步所用軟件會不會進行標準化。
在轉(zhuǎn)錄本水平上，一般常用工具為Cufflinks和它的繼任者StringTie， eXpress。這些軟件要處理的難題就時轉(zhuǎn)錄本亞型（isoforms）之間通常是有重疊的，當二代測序讀長低于轉(zhuǎn)錄本長度時，如何進行區(qū)分？這些工具大多采用的都是expectation maximization（EM）。好在我們有三代測序。
上述軟件都是alignment-based，目前許多alignment-free軟件，如kallisto, silfish, salmon，能夠省去比對這一步，直接得到read count，在運行效率上更高。不過最近一篇文獻[1]指出這類方法在估計豐度時存在樣本特異性和讀長偏差。
在外顯子使用水平上，其實和基因水平的統(tǒng)計類似。但是值得注意的是為了更好的計數(shù)，我們需要提供無重疊的外顯子區(qū)域的gtf文件[2]。用于分析差異外顯子使用的DEXSeq提供了一個Python腳本（dexseq_prepare_annotation.py）執(zhí)行這個任務。
小結(jié)
計數(shù)分為三個水平： gene-level, transcript-level, exon-usage-level
標準化方法： FPKM RPKM TMM TPM
輸出表達矩陣
在RNA-Seq分析中，每一個基因就是一個feature（特征？），而基因被認為是它的所有外顯子的和集。在可變剪切分析中，可以單獨把每個外顯子當作一個feature。而在ChIP-Seq分析中，feature則是預先定義的結(jié)合域。但是確定一個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] htseq-count -r pos -f bam RNA-Seq/aligned/SRR3589957_sorted.bam reference/gencode.v26lift37.annotation.sorted.gtf > SRR3589957.count
用一個循環(huán)處理多個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}.logdone
運行的時間會比較久，所以可以去了解不同參數(shù)的用法了，其中比較常用的為：
-f bam/sam： 指定輸入文件格式，默認SAM
-r name/pos: 你需要利用samtool sort對數(shù)據(jù)根據(jù)read name或者位置進行排序，默認是name
-s yes/no/reverse: 數(shù)據(jù)是否來自于strand-specific assay。DNA是雙鏈的，所以需要判斷到底來自于哪條鏈。如果選擇了no， 那么每一條read都會跟正義鏈和反義鏈進行比較。默認的yes對于雙端測序表示第一個read都在同一個鏈上，第二個read則在另一條鏈上。
-a 最低質(zhì)量， 剔除低于閾值的read
-m 模式 union（默認）, intersection-strict and intersection-nonempty。一般而言就用默認的，作者也是這樣認為的。
-i id attribute: 在GTF文件的最后一欄里，會有這個基因的多個命名方式（如下）， RNA-Seq數(shù)據(jù)分析常用的是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測序的結(jié)果，但是我們其實只下載了3個sra文件。一般而言RNA-Seq數(shù)據(jù)分析都至少要有2個重復，所以必須要有4個sra文件才行。我在仔細讀完文章的方法這一段以后，發(fā)現(xiàn)他們有一批數(shù)據(jù)用的是其他課題組的： 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交流之后，他也承認自己只分析了小鼠的數(shù)據(jù)，而沒有分析人類的數(shù)據(jù)。所以我們需要根據(jù)文章提供的線索下載另外一份數(shù)據(jù)，才能進行下一步的分析。
這個時候就有一個經(jīng)常被問到的問題：不同來源的RNA-Seq數(shù)據(jù)能夠直接比較嗎？甚至說如果不同來源的RNA-seq數(shù)據(jù)的構(gòu)建文庫都不一樣該如何比較?不同來源的RNA-Seq結(jié)果之間比較需要考慮?批次效應（batch effect)?的影響。
處理批次效應，根據(jù)我搜索的結(jié)果，是不能使用FPKM/RPKM，關(guān)于這個標準化的吐槽，我在biostars上找到了如下觀點：
FPKM/RPKM 不是標準化的方法，它會引入文庫特異的協(xié)變量
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結(jié)果會存在很大差異。
結(jié)論： 可以問下原作者他們是如何處理數(shù)據(jù)的，居然有一個居然沒有重復的分析也能過審。改用小鼠數(shù)據(jù)進行分析?；蛘呤褂脽o重復的分析方法，或者模擬一份數(shù)據(jù)出來，先把流程走完。
合并表達矩陣
HTSeq-count輸出結(jié)果是一個個獨立的文件，后續(xù)分析需要把多個文件合并成一個行為基因名，列為樣本名，中間為count的行列式文件。肯定是不會用Excel手動處理（雖然可以寫一個VBA腳本，但是數(shù)據(jù)量過大不好處理了），這里使用的Python寫一個腳本。
基本邏輯：
讀取文件
建立一個字典，如果key不在字典中，新增key和value，如果key在字典中，追加value。
輸出
#!/usr/bin/python3import sysmydict = {}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] = valuefor key,value in mydict.items():? ? print(key + '\t' + value +'\n')
這幾行代碼寫了2個番茄鐘，但是debug花了我一個番茄鐘。問題出在str和list兩種數(shù)據(jù)格式的混亂使用。還有一個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讀取數(shù)據(jù)。
1.導入數(shù)據(jù)
options(stringsAsFactors =?FALSE)# import data if sample are smallcontrol <- 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.數(shù)據(jù)整合。gencode的注釋文件中的gene_id（如ENSG00000105298.13_3）在EBI是不能搜索到的，所以我就只保留ENSG00000105298這部分。
# merge data and delete the unuseful rowraw_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，所以可以刪掉節(jié)省體積。
summary(raw_count_filt)? ? control? ?? ?? ?? ???rep1? ?? ?? ?? ?? ?rep2? ?? ?? ? Min.? ?:0.0Min.? ?:0.0Min.? ?:0.01st Qu.:0.01st Qu.:0.01st Qu.:0.0Median :0.0Median :0.0Median :0.0Mean? ?:356.1Mean? ?:370.3Mean? ?:316.63rd Qu.:?15.03rd Qu.:?15.03rd Qu.:?10.0Max.? ?:161867.0Max.? ?:121902.0Max.? ?:105565.0
4.看幾個具體基因
在EBI上搜索GAPDH找到ID為ENSG00000111640。
GAPDH <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000111640",]? ?? ?? ?? ?? ?? ?? ?? ?? ???gene_id control??rep1??rep2ENSG00000111640 ENSG00000111640.14_2? ?41857 53902 55302
文章研究的AKAP95（ENSG00000105127）的表達量在KD中都降低了
> AKAP95 <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000105127",]> AKAP95? ?? ?? ?? ?? ?? ?? ?? ?? ? gene_id control rep1 rep2ENSG00000105127 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)
[size=0em]?
