按：RNA-seq 流程非常多樣化，每一個步驟可選工具都非常多。我的選擇是 StringTie 進(jìn)行組裝取得 read counts 數(shù)值覆山，DESeq2 分析差異基因在张。這里把我方法跟大家共享。一來給新手一個指引爆土，如果你毫無相關(guān)經(jīng)驗可以照著我這個出結(jié)果；二來拋磚引玉，如果你有相關(guān)經(jīng)驗歡迎一起交流進(jìn)步慨飘。
總體流程
如流程圖所示，總體上有三大步驟译荞。第一瓤的，用 stringtie 組裝取得每個樣本結(jié)果 gtf 文件；第二吞歼，用 StringTie 提供的 prepDE.py 腳本從所有(需要)樣本提取 read counts 數(shù)值并以矩陣格式存儲到一個 csv 文件圈膏；第三，用 DESeq2 進(jìn)行差異基因分析浆熔。
StringTie
StringTie 的輸入文件是排序好的 BAM 格式文件本辐，所以比對得到 sam 文件后用 samtools 進(jìn)行 sam->bam 格式轉(zhuǎn)換，并對 bam 進(jìn)行排序和建立索引文件医增。組裝需要一個基因組注釋文件慎皱，推薦去 genecode 下載。網(wǎng)址 GENCODE - Human Release 33叶骨。如下圖所示下載 gtf/gff3 格式茫多。
示例 stringtie 命令如下：
stringtie /Example/Bam/PC3_NC_1_sorted.bam -G /Example/GRCh38_hg38/gencode.v32.chr_patch_hapl_scaff.annotation.gff3 -o /Example/GTF/PC3_NC_1/PC3_NC_1.gtf -A /Example/GTF/PC3_NC_1/PC3_NC_1_abundance.txt -p 4 -B -e
其中 -G 參數(shù)指定基因組注釋文件， -o 輸出的 gtf 路徑忽刽， -e 參數(shù)限定分析基因組注釋文件轉(zhuǎn)錄本天揖，也就是說只分析已知轉(zhuǎn)錄本，不分析新轉(zhuǎn)錄本跪帝。 -B 是產(chǎn)生文件給下游 Ballgown 使用今膊，我們不使用 Ballgown 的話可以不需要指定。得到 gtf 后用 prepDE.py 轉(zhuǎn)換到 read counts 矩陣伞剑。 prepDE.py 輸入是目錄路徑斑唬，或者提供一個文件，文件里是你需要的 gtf 文件路徑。像我當(dāng)前目錄包含了需要的 gtf 文件恕刘。那么輸入當(dāng)前目錄路徑給 prepDE.py 即可缤谎。
$find . -name "*.gtf"
./PC3_si730_2/PC3_si730_2.gtf
./PC3_si730_1/PC3_si730_1.gtf
./PC3_NC_2/PC3_NC_2.gtf
./PC3_NC_1/PC3_NC_1.gtf
./PC3_NC_3/PC3_NC_3.gtf
./PC3_si730_3/PC3_si730_3.gtf
假設(shè)就在當(dāng)前目錄運行，示例命令如下：
python prepDE.py -i . -l 150
其中 -l 是設(shè)置平均 reads 長度褐着。 prepDE.py 計算 read counts 代碼如下坷澡，read_len 就是 -l 參數(shù)『兀可見這個參數(shù)需要按實際實驗去設(shè)置频敛，默認(rèn)值是75.
int(ceil(coverage*transcript_len/read_len)
運行后得到 gene_count_matrix.csv 和 transcript_count_matrix.csv 兩個文件，其中 gene_count_matrix.csv 就是需要的基因 read counts 矩陣谴餐。查看一下姻政。
gene_id,PC3_NC_1,PC3_NC_2,PC3_NC_3,PC3_si730_1,PC3_si730_2,PC3_si730_3
ENSG00000255427.1|LINC02707,0,0,0,0,0,0
ENSG00000227189.2|AC092535.2,0,2,0,5,0,0
ENSG00000244171.4|PBX2P1,1782,1921,1717,1320,1829,1704
ENSG00000122188.13|LAX1,3,3,0,3,0,0
ENSG00000225415.2|ACKR4P1,0,0,0,0,0,0
ENSG00000141736.13|ERBB2,4397,4340,3813,3544,4978,4346
ENSG00000284348.1|AC245128.33,0,0,0,0,0,0
ENSG00000273728.1|Metazoa_SRP,0,0,0,0,0,0
可見其 gene_id 是混合了 ensembl 和其他基因名，用 | 分隔岂嗓，后續(xù) DESeq2 流程我們需要根據(jù)這個進(jìn)行相應(yīng)處理汁展。關(guān)于 ensembl id 可以查看我之前文章《Ensembl基因ID格式》了解。
DESeq2
DESeq2 是一個 R 包厌殉，每個人 R 代碼習(xí)慣不同食绿，所以下面每一步我都會展示結(jié)果/數(shù)據(jù)， 你不需要跟著我代碼寫公罕，只要能做到相同數(shù)據(jù)結(jié)構(gòu)就行器紧。 DESeq2 是非常復(fù)雜強大的包，建議去官網(wǎng)把教程至少看一遍楼眷。
第一步先導(dǎo)入需要的 R 包铲汪。其中 biomaRt 是用于基因ID轉(zhuǎn)換，包使用方法參見我以前文章《批量轉(zhuǎn)換基因ID》罐柳。
library(tidyverse, quietly = TRUE)
library(DESeq2, quietly = TRUE)
library(biomaRt, quietly = TRUE)
DESeq2 的輸入是 read counts 矩陣和樣品信息表掌腰，樣品信息主要是樣品分組信息。在我這里是個簡單的敲低實驗张吉，其中3個對照組3個實驗組齿梁，分析要求是取得 實驗組 vs 對照組 差異基因。我的樣品信息就存在這樣文件里：
$head SampleGroup.tsv
Group
PC3_NC_1 Control
PC3_NC_2 Control
PC3_NC_3 Control
PC3_si730_1 Treatment
PC3_si730_2 Treatment
PC3_si730_3 Treatment
使用下面代碼將表格讀取進(jìn) R：
sampleGroup <- read.table("/Example/SampleGroup.tsv", header = TRUE)
sampleGroup$Group <- factor(sampleGroup$Group, levels = c("Control", "Treatment"))
print(sampleGroup)
這里我特意把分組列轉(zhuǎn)換到因子肮蛹，并且設(shè)定第一個因子水平為 Control. 這是用 R 的好習(xí)慣勺择，在 R 語言里如果涉及到排序，往往默認(rèn)按因子水平排伦忠。如果你的列表不是因子省核，會自動 as.factor 進(jìn)行轉(zhuǎn)換，采用默認(rèn)的因子水平昆码。返回 DESeq2 這里芳撒，把 Control 設(shè)為第一個因子邓深，那么后面的結(jié)果就是 其他組 vs 對照組。在我這個例子里就是 Treatment vs Control 笔刹。假如說你有2組 A 和 C ， C 是你的對照組冬耿，你希望是 A vs C 舌菜。如果你不設(shè)定 C 為第一個因子，那 DESeq2 默認(rèn) A 為第一個因子亦镶，結(jié)果就是 C vs A 日月。
使用下面代碼讀取 read counts 矩陣。
readCount <- read_csv("/Example/GTF/gene_count_matrix.csv") %>% tidyr::separate(col = gene_id, into = c("ensembl_gene_id", "gene_name"), sep = "\\|", remove = TRUE) %>% dplyr::select(ensembl_gene_id, rownames(sampleGroup)) %>% dplyr::mutate(ensembl_gene_id = map_chr(ensembl_gene_id, geneID)) %>% as.data.frame()
rownames(readCount) <- readCount$ensembl_gene_id
readCount$ensembl_gene_id <- NULL
print(head(readCount, n = 3))
這里 StringTie 給的結(jié)果基因ID是 ENSG00000214428.3|NPM1P10 這樣的格式缤骨，其中豎線右邊基因名是雜亂的爱咬，為了方便處理我把這部分移除了，只保留 Ensembl ID: ENSG00000214428 绊起。后面再用 biomaRt 轉(zhuǎn)換到 hgnc_symbol 精拟。總而言之虱歪，處理后得到下面的數(shù)據(jù)框蜂绎，這就是 DESeq2 需要的輸入格式。要注意列名順序跟樣品信息 sampleGroup 的行名順序是一致的笋鄙。
PC3_NC_1 PC3_NC_2 PC3_NC_3 PC3_si730_1 PC3_si730_2 PC3_si730_3
ENSG00000255427 0 0 0 0 0 0
ENSG00000227189 0 2 0 5 0 0
ENSG00000244171 1663 1793 1602 1232 1707 1591
剛剛說了要用 biomaRt 進(jìn)行基因ID轉(zhuǎn)換师枣，下面是代碼和結(jié)果示例，詳情參見《批量轉(zhuǎn)換基因ID》文章萧落。
# 其中 is.na(hgnc_symbol) 或許應(yīng)該替換為 hgnc_symbol == ""
ensembl <- useMart("ensembl", dataset="hsapiens_gene_ensembl")
geneMap <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol", "entrezgene_id"), filters = "ensembl_gene_id", values = rownames(readCount), mart = ensembl) %>% dplyr::filter(!(is.na(hgnc_symbol) & is.na(entrezgene_id))) %>% dplyr::distinct(ensembl_gene_id, .keep_all = TRUE)
print(head(geneMap, n = 3))
# 得到如下表格
ensembl_gene_id hgnc_symbol entrezgene_id
1 ENSG00000010404 IDS 3423
2 ENSG00000013293 SLC7A14 57709
3 ENSG00000029534 ANK1 286
使用下面代碼將數(shù)據(jù)導(dǎo)入 DESeq2 践美。
dds <- DESeqDataSetFromMatrix(countData = readCount, colData = sampleGroup, design = ~ Group)
keep <- rowSums(counts(dds)) > 12
dds <- dds[keep, ]
dds <- DESeq(dds)
print(resultsNames(dds))
# 打印輸出
[1] "Intercept" "Group_Treatment_vs_Control"
函數(shù) DESeqDataSetFromMatrix 將 read counts 矩陣和樣品信息導(dǎo)入到 DESeq2 。其中 design 參數(shù)設(shè)定實驗設(shè)計找岖，也即那一列信息用于分組比較(也包括批次效應(yīng)或其他因素陨倡，在我這例子里沒有)。后面的 keep 是移除低 read counts 基因宣增，這里也沒有一個數(shù)值標(biāo)準(zhǔn)玫膀，看自己想法。我這里是要求每個基因平均 read counts > 2 爹脾，我有6樣本所以要求總 read counts > 12 帖旨。只保留符合條件的基因。像 edgeR 的話是推薦按照 cpm(count-per-million) 進(jìn)行篩選而不是直接對 read counts 篩選灵妨，比如 keep <- rowSums(cpm(dds)>1) >= 2 也可以參考一下這個方法解阅。
用 results 函數(shù)取的差異基因結(jié)果，這個結(jié)果包含所有基因泌霍。
res <- results(dds, name = "Group_Treatment_vs_Control") %>% as_tibble(rownames = "ensembl_gene_id") %>% dplyr::left_join(geneMap, by = "ensembl_gene_id") %>% dplyr::select(ensembl_gene_id, hgnc_symbol, entrezgene_id, everything())
print(head(res, n = 3))
用 name 參數(shù)選擇需要的比較組货抄，也可以用 contrast 參數(shù)設(shè)置比較組述召。results 函數(shù)有個 lfcThreshold 參數(shù)跟P值計算直接相關(guān)，默認(rèn)是0所以檢驗的是基因表達(dá)有無差異蟹地，如果你想檢驗基因表達(dá)是否差2倍以上积暖，可以設(shè)為1。同時我們把基因ID匹配上去怪与。得到差異基因結(jié)果如下：
# A tibble: 3 x 9
ensembl_gene_id hgnc_symbol entrezgene_id baseMean log2FoldChange lfcSE stat
<chr> <chr> <int> <dbl> <dbl> <dbl> <dbl>
1 ENSG00000244171 PBX2P1 NA 1583. -0.123 0.118 -1.04
2 ENSG00000141736 ERBB2 2064 3916. 0.0734 0.105 0.702
3 ENSG00000268324 LRRC2-AS1 NA 15.8 1.20 0.860 1.40
# … with 2 more variables: pvalue <dbl>, padj <dbl>
得到結(jié)果后可以根據(jù)自己需要設(shè)置條件篩選夺刑，比如：
degFilter <- dplyr::filter(res, abs(log2FoldChange) >= 1 & padj < 0.05)
DESeq2 提供了log2差異倍數(shù)(log2 fold change)收縮方法 lfcShrink ，跟原始的 results 取的結(jié)果相比會把差異倍數(shù)進(jìn)行一定縮小分别。 lfcShrink 官方文檔如下：
Beginning with the first row, all shrinkage methods provided by DESeq2 are good for ranking genes by “effect size”, that is the log2 fold change (LFC) across groups, or associated with an interaction term. It is useful to contrast ranking by effect size with ranking by a p-value or adjusted p-value associated with a null hypothesis: while increasing the number of samples will tend to decrease the associated p-value for a gene that is differentially expressed, the estimated effect size or LFC becomes more precise. Also, a gene can have a small p-value although the change in expression is not great, as long as the standard error associated with the estimated LFC is small.
總體而言縮放數(shù)據(jù)更適合去可視化以及進(jìn)行排序等遍愿，像我們后面需要進(jìn)行GSEA分析，那么可以采用這個數(shù)據(jù)耘斩。有3種縮放方法沼填，官方給出下面總結(jié)：
示例代碼：
deg <- lfcShrink(dds, coef = "Group_Treatment_vs_Control", type = "apeglm") %>% as_tibble(rownames = "ensembl_gene_id") %>% dplyr::left_join(geneMap, by = "ensembl_gene_id") %>% dplyr::select(ensembl_gene_id, hgnc_symbol, entrezgene_id, everything())
print(head(deg, n = 3))
# 打印輸出
# A tibble: 3 x 8
ensembl_gene_id hgnc_symbol entrezgene_id baseMean log2FoldChange lfcSE
<chr> <chr> <int> <dbl> <dbl> <dbl>
1 ENSG00000244171 PBX2P1 NA 1583. -0.0000311 0.00144
2 ENSG00000141736 ERBB2 2064 3916. -0.00000352 0.00144
3 ENSG00000268324 LRRC2-AS1 NA 15.8 0.00000157 0.00144
# … with 2 more variables: pvalue <dbl>, padj <dbl>
除了差異基因數(shù)據(jù)，往往我們還需要表達(dá)數(shù)據(jù)用于其他分析括授。 DESeq2 提供 rlog 函數(shù)將原始 read counts 轉(zhuǎn)換為 log2 后值并且進(jìn)行調(diào)整坞笙。得到的矩陣適合后續(xù)機(jī)器學(xué)習(xí)相關(guān)分析，比如聚類刽脖、熱圖展示等羞海。如果你樣本數(shù)很多，建議使用 vst 函數(shù)曲管，會快非常多却邓。代碼示例。
rld <- rlog(dds)
rldt <- assay(rld) %>% as_tibble(rownames = "ensembl_gene_id") %>% dplyr::left_join(geneMap, by = "ensembl_gene_id") %>% dplyr::select(ensembl_gene_id, hgnc_symbol, entrezgene_id, everything())
print(head(rldt, n = 3))
# 輸出
# A tibble: 3 x 9
ensembl_gene_id hgnc_symbol entrezgene_id PC3_NC_1 PC3_NC_2 PC3_NC_3
<chr> <chr> <int> <dbl> <dbl> <dbl>
1 ENSG00000244171 PBX2P1 NA 10.6 10.7 10.7
2 ENSG00000141736 ERBB2 2064 11.9 11.9 11.9
3 ENSG00000268324 LRRC2-AS1 NA 3.84 3.84 3.87
# … with 3 more variables: PC3_si730_1 <dbl>, PC3_si730_2 <dbl>,
# PC3_si730_3 <dbl>
使用 counts 函數(shù)取的 normalized 后 read counts 數(shù)據(jù),記得參數(shù) normalized = TRUE .
nReadCount <- counts(dds, normalized = TRUE) %>% as_tibble(rownames = "ensembl_gene_id") %>% dplyr::left_join(geneMap, by = "ensembl_gene_id") %>% dplyr::select(ensembl_gene_id, hgnc_symbol, entrezgene_id, everything())
print(head(nReadCount, n = 3))
# 輸出
# A tibble: 3 x 9
ensembl_gene_id hgnc_symbol entrezgene_id PC3_NC_1 PC3_NC_2 PC3_NC_3
<chr> <chr> <int> <dbl> <dbl> <dbl>
1 ENSG00000244171 PBX2P1 NA 1575. 1678. 1698.
2 ENSG00000141736 ERBB2 2064 3888. 3789. 3772.
3 ENSG00000268324 LRRC2-AS1 NA 8.53 8.42 11.7
# … with 3 more variables: PC3_si730_1 <dbl>, PC3_si730_2 <dbl>,
# PC3_si730_3 <dbl>
