要求

實(shí)現(xiàn)這個(gè)功能的軟件也很多扫皱，還是煩請大家先自己搜索幾個(gè)教程足绅，入門請統(tǒng)一用htseq-count捷绑，對每個(gè)樣本都會(huì)輸出一個(gè)表達(dá)量文件。

需要用腳本合并所有的樣本為表達(dá)矩陣氢妈。參考：生信編程直播第四題：多個(gè)同樣的行列式文件合并起來

對這個(gè)表達(dá)矩陣可以自己簡單在excel或者R里面摸索粹污，求平均值，方差首量。

看看一些生物學(xué)意義特殊的基因表現(xiàn)如何壮吩，比如GAPDH,β-ACTIN等等。

理論基礎(chǔ)

在上篇的比對中加缘，我們需要糾結(jié)是否真的需要比對鸭叙，如果你只需要知道已知基因的表達(dá)情況，那么可以選擇alignment-free工具（例如salmon, sailfish)拣宏，如果你需要找到noval isoforms沈贝，那么就需要alignment-based工具（如HISAT2, STAR）。到了這一篇的基因（轉(zhuǎn)錄本）定量勋乾，需要考慮的因素就更加多了宋下，以至于我不知道如何說清才能理清邏輯仙逻。

定量分為三個(gè)水平

基因水平(gene-level)

轉(zhuǎn)錄本水平(transcript-level)

外顯子使用水平(exon-usage-level)伯病。

在基因水平上规阀，常用的軟件為HTSeq-count队伟，featureCounts葵礼，BEDTools, Qualimap, Rsubread, GenomicRanges等颂龙。以常用的HTSeq-count為例葵第，這些工具要解決的問題就是根據(jù)read和基因位置的overlap判斷這個(gè)read到底是誰家的孩子税弃。值得注意的是不同工具對multimapping reads處理方式也是不同的绅你，例如HTSeq-count就直接當(dāng)它們不存在伺帘。而Qualimpa則是一人一份，平均分配忌锯。

_images/count_modes.png

對每個(gè)基因計(jì)數(shù)之后得到的count matrix再后續(xù)的分析中伪嫁，要注意標(biāo)準(zhǔn)化的問題。如果你要比較同一個(gè)樣本(within-sample)不同基因之間的表達(dá)情況偶垮，你就需要考慮到轉(zhuǎn)錄本長度张咳，因?yàn)檗D(zhuǎn)錄本越長，那么檢測的片段也會(huì)更多似舵，直接比較等于讓小孩和大人進(jìn)行賽跑脚猾。如果你是比較不同樣本（across sample）同一個(gè)基因的表達(dá)情況，雖然不必在意轉(zhuǎn)錄本長度砚哗，但是你要考慮到測序深度（sequence depth)龙助，畢竟測序深度越高，檢測到的概率越大蛛芥。除了這兩個(gè)因素外提鸟，你還需要考慮GC%所導(dǎo)致的偏差军援，以及測序儀器的系統(tǒng)偏差。目前對read count標(biāo)準(zhǔn)化的算法有RPKM（SE）, FPKM（PE）称勋，TPM, TMM等胸哥，不同算法之間的差異與換算方法已經(jīng)有文章進(jìn)行整理和吐槽了。但是赡鲜，有一些下游分析的軟件會(huì)要求是輸入的count matrix是原始數(shù)據(jù)空厌，未經(jīng)標(biāo)準(zhǔn)化，比如說DESeq2银酬，這個(gè)時(shí)候你需要注意你上一步所用軟件會(huì)不會(huì)進(jìn)行標(biāo)準(zhǔn)化嘲更。

在轉(zhuǎn)錄本水平上，一般常用工具為Cufflinks和它的繼任者StringTie捡硅， eXpress哮内。這些軟件要處理的難題就時(shí)轉(zhuǎn)錄本亞型（isoforms）之間通常是有重疊的，當(dāng)二代測序讀長低于轉(zhuǎn)錄本長度時(shí)壮韭，如何進(jìn)行區(qū)分北发？這些工具大多采用的都是expectation maximization（EM）。好在我們有三代測序喷屋。

上述軟件都是alignment-based琳拨，目前許多alignment-free軟件，如kallisto, silfish, salmon屯曹，能夠省去比對這一步狱庇，直接得到read count，在運(yùn)行效率上更高恶耽。不過最近一篇文獻(xiàn)[1]指出這類方法在估計(jì)豐度時(shí)存在樣本特異性和讀長偏差密任。

在外顯子使用水平上，其實(shí)和基因水平的統(tǒng)計(jì)類似偷俭。但是值得注意的是為了更好的計(jì)數(shù)浪讳，我們需要提供無重疊的外顯子區(qū)域的gtf文件[2]。用于分析差異外顯子使用的DEXSeq提供了一個(gè)Python腳本（dexseq_prepare_annotation.py）執(zhí)行這個(gè)任務(wù)涌萤。

小結(jié)

計(jì)數(shù)分為三個(gè)水平： gene-level, transcript-level, exon-usage-level

標(biāo)準(zhǔn)化方法： FPKM RPKM TMM TPM

輸出表達(dá)矩陣

在RNA-Seq分析中淹遵，每一個(gè)基因就是一個(gè)feature（特征？）负溪，而基因被認(rèn)為是它的所有外顯子的和集透揣。在可變剪切分析中，可以單獨(dú)把每個(gè)外顯子當(dāng)作一個(gè)feature川抡。而在ChIP-Seq分析中辐真，feature則是預(yù)先定義的結(jié)合域。但是確定一個(gè)read到底屬于哪一個(gè)feature有時(shí)會(huì)非常棘手。因此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

用一個(gè)循環(huán)處理多個(gè)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

運(yùn)行的時(shí)間會(huì)比較久恨闪，所以可以去了解不同參數(shù)的用法了倘感，其中比較常用的為：

-f bam/sam： 指定輸入文件格式放坏，默認(rèn)SAM

-r name/pos: 你需要利用samtool sort對數(shù)據(jù)根據(jù)read name或者位置進(jìn)行排序，默認(rèn)是name

-s yes/no/reverse: 數(shù)據(jù)是否來自于strand-specific assay老玛。DNA是雙鏈的淤年，所以需要判斷到底來自于哪條鏈。如果選擇了no蜡豹， 那么每一條read都會(huì)跟正義鏈和反義鏈進(jìn)行比較麸粮。默認(rèn)的yes對于雙端測序表示第一個(gè)read都在同一個(gè)鏈上，第二個(gè)read則在另一條鏈上镜廉。

-a 最低質(zhì)量弄诲， 剔除低于閾值的read

-m 模式 union（默認(rèn)）, intersection-strict and intersection-nonempty。一般而言就用默認(rèn)的娇唯，作者也是這樣認(rèn)為的齐遵。

-i id attribute: 在GTF文件的最后一欄里，會(huì)有這個(gè)基因的多個(gè)命名方式（如下）塔插， RNA-Seq數(shù)據(jù)分析常用的是gene_id梗摇， 當(dāng)然你可以寫一個(gè)腳本替換成其他命名方式。

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é)果想许，但是我們其實(shí)只下載了3個(gè)sra文件伶授。一般而言RNA-Seq數(shù)據(jù)分析都至少要有2個(gè)重復(fù)，所以必須要有4個(gè)sra文件才行流纹。我在仔細(xì)讀完文章的方法這一段以后糜烹，發(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交流之后漱凝，他也承認(rèn)自己只分析了小鼠的數(shù)據(jù)疮蹦，而沒有分析人類的數(shù)據(jù)。所以我們需要根據(jù)文章提供的線索下載另外一份數(shù)據(jù)碉哑，才能進(jìn)行下一步的分析挚币。

這個(gè)時(shí)候就有一個(gè)經(jīng)常被問到的問題：不同來源的RNA-Seq數(shù)據(jù)能夠直接比較嗎？甚至說如果不同來源的RNA-seq數(shù)據(jù)的構(gòu)建文庫都不一樣該如何比較?不同來源的RNA-Seq結(jié)果之間比較需要考慮批次效應(yīng)（batch effect)的影響扣典。

處理批次效應(yīng)妆毕，根據(jù)我搜索的結(jié)果，是不能使用FPKM/RPKM贮尖，關(guān)于這個(gè)標(biāo)準(zhǔn)化的吐槽笛粘，我在biostars上找到了如下觀點(diǎn)：

FPKM/RPKM 不是標(biāo)準(zhǔn)化的方法，它會(huì)引入文庫特異的協(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

有人建議使用一個(gè)Bioconductor包http://www.bioconductor.org/packages/devel/bioc/html/sva.html我沒有具體了解，有生之年去了解補(bǔ)充薪前。

還有人引用了一篇文獻(xiàn) IVT-seq reveals extreme bias in RNA-sequencing 證明不同文庫的RNA-Seq結(jié)果會(huì)存在很大差異润努。

結(jié)論： 可以問下原作者他們是如何處理數(shù)據(jù)的，居然有一個(gè)居然沒有重復(fù)的分析也能過審示括。改用小鼠數(shù)據(jù)進(jìn)行分析铺浇。或者使用無重復(fù)的分析方法垛膝，或者模擬一份數(shù)據(jù)出來鳍侣，先把流程走完。

合并表達(dá)矩陣

HTSeq-count輸出結(jié)果是一個(gè)個(gè)獨(dú)立的文件吼拥，后續(xù)分析需要把多個(gè)文件合并成一個(gè)行為基因名倚聚，列為樣本名，中間為count的行列式文件凿可』笳郏肯定是不會(huì)用Excel手動(dòng)處理（雖然可以寫一個(gè)VBA腳本，但是數(shù)據(jù)量過大不好處理了）枯跑，這里使用的Python寫一個(gè)腳本惨驶。

基本邏輯：

讀取文件

建立一個(gè)字典，如果key不在字典中全肮，新增key和value敞咧，如果key在字典中，追加value辜腺。

輸出

#!/usr/bin/python3importsysmydict = {}forfileinsys.argv[1:]:forlineinopen(file,'r'):? ? ? ? key,value = line.strip().split('\t')ifkeyinmydict:? ? ? ? ? ? mydict[key] = mydict[key] +'\t'+ valueelse:? ? ? ? ? ? mydict[key] = valueforkey,valueinmydict.items():? ? print(key +'\t'+ value +'\n')

這幾行代碼寫了2個(gè)番茄鐘休建，但是debug花了我一個(gè)番茄鐘。問題出在str和list兩種數(shù)據(jù)格式的混亂使用评疗。還有一個(gè)bug： 由于詞典是無序的测砂，所以原本代表樣本來源的第一行，會(huì)跑到其他行百匆。

在論壇上找到一個(gè)更加簡潔的代碼（要求基因名順序排列）砌些，用到paste,awk,printf這三個(gè)shell命令。

paste *.txt | awk'{printf $1 "\t";for(i=2;i<=NF;i+=2) printf $i"\t";printf $i}'

保存為countCombiner.py加匈，輸入文件為count, 輸出為標(biāo)準(zhǔn)輸出存璃，需要重定向。

簡單分析

這一步需要用到R語言或者是Excel讀取數(shù)據(jù)雕拼。

1.導(dǎo)入數(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)controlrep1rep2Min.? :? ? 0.0Min.? :? ? 0.0Min.? :? ? 0.01stQu.:? ? 0.01stQu.:? ? 0.01stQu.:? ? 0.0Median:? ? 0.0Median:? ? 0.0Median:? ? 0.0Mean:? 356.1Mean:? 370.3Mean:? 316.63rdQu.:? ? 15.03rdQu.:? ? 15.03rdQu.:? ? 10.0Max.:161867.0Max.:121902.0Max.:105565.0

4.看幾個(gè)具體基因

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

GAPDH<-raw_count_filt[rownames(raw_count_filt)=="ENSG00000111640",]gene_idcontrolrep1rep2ENSG00000111640ENSG00000111640.14_241857 53902 55302

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

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

下面的差異基因表達(dá)衰絮，讓我想想袍冷，該如何收拾Jimmy挖的坑。

參考文獻(xiàn)

[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)

作者：hoptop

鏈接：http://www.reibang.com/p/e9742bbf83b9

來源：簡書

簡書著作權(quán)歸作者所有猫牡，任何形式的轉(zhuǎn)載都請聯(lián)系作者獲得授權(quán)并注明出處胡诗。