120分的轉錄組試題和答案
120分的轉錄組試題和答案
這個答案之前出過三份，最近整理了一份文字版，方便觀看，還請大家多多補充。
一、理論題目
1、說出至少5種高通量測序技術的應用，並簡要描述其在科研中的用途？ (2分) (選答)
-
de novo sequencing：獲得該物種的參考序列，為後續研究奠定基礎。
-
基因組重測序：進行全基因組重測序，在全基因組水平上掃描並檢測突變位點，發現個體差異的分子基礎。
-
轉錄組測序：基因表達定量，尋找差異基因；發現新的轉錄本、可變剪下、基因融合等。(
轉錄組分析的正確姿勢 ) -
ChIP-Seq: 結合ChIP和高通量測序技術的優勢，能夠在全基因組範圍內高效地研究目的蛋白的結合位點。( ChIP-seq基本分析流程 )
-
DNA甲基化：基因組水平高精度的甲基化檢測成為現實。甲基化DNA免疫共沉澱測序、甲基結合蛋白測序、亞硫酸氫鹽測序。
-
測序發展史：150年的風雨歷程 也看一下。
2、轉錄組測序可以解決的問題有哪些？(2分)
-
可以從轉錄組水平整體評估外界條件引起的基因表達響應。
-
對於檢測低丰度轉錄本和發現新轉錄本具有獨特優勢。
-
可以發現基因融合、可變剪下及結構變異。
3、轉錄組專案開展前需要考慮哪些問題？ (6分)
-
明確解決的生物學問題，選擇合適的建庫方法和測序方法。
-
設計合理的對照樣品和生物學重複。
-
儘量一次性完成樣本的處理和測序。不然就得處理 批次效應 。
-
經濟條件，樣本獲取難易成度，失敗風險。
4、轉錄組測序同一樣品不同重複之間相關性多少合適？分別敘述不同的需求可用reads數多少合適？ (3分)
-
相同基因型生物學重複Spearman係數>0.9；不同基因型生物學重複Spearman係數>0.8。
-
檢測新轉錄本，可變剪下： 50-100 million可用reads.
-
Long RNA-Seq library: 30 million aligned reads.
-
RAMPAGE library: 30 million aligned reads.
-
Small RNA-Seq library: 30 million aligned reads.
5、測序reads數和測序鹼基數之間如何換算?(1分)
- Reads數=鹼基數/讀長；例子：如果讀長150bp 鹼基數4.5G，則對應的reads為30M。即：30million reads=4.5G/150bp
6、描述下常規轉錄組測序的過程，從樣本準備到獲取測序資料中經歷的步驟？ (6分)
- 獲取樣本：培養細胞或從組織上進行解剖，必要時可以用流式細胞儀分離細胞，注意對照組和生物學重複。
- RNA提取：根據不同的提取試劑盒提total RNA，去除rRNA或富集mRNA，反轉錄為cDNA.
- 建庫：將cDNA片段化(超聲剪下或酶切)，末端修復+A，連結index接頭，一般需要再進行5-8個迴圈的文庫擴增，磁珠純化進行片段篩選。
- 上機測序：把index不同的樣品混合到一起，把雙鏈文庫變性成單鏈，與flowcell上的oligo互補結合，進行橋式擴增，形成cluster，進行邊合成邊測序 (SBS)。具體見 NGS基礎 - 高通量測序原理 。
- 獲取資料：按照接頭拆分下機後的bcl格式文件並轉為fastq格式文件。
7、描述下常規轉錄組分析的流程，列舉每步可用的工具、每步分析的意義? (10分)
- 資料質控 ：質量評估用FastQC軟體，質控cutadapt，Trimmomatic。意義：去除鹼基質量低的reads和接頭部分。
- 比對： Tophat2/HISAT2/STAR 。意義：確定每條reads來源於哪個基因，什麼位置。
- 差異基因鑑定： DEseq2 。意義：比較標準化後的基因表達丰度，得到差異基因。
- GO富集分析： GOEAST 。 意義：尋找差異基因與哪些功能有關。
- 轉錄本拼裝：String Tie；可變剪下：Rmats。意義：發現新的轉錄本和可變剪下形式。
- 通路視覺化 、 WGCNA共表達分析 等。
8、基因組和轉錄組測序序列比對使用的工具有什麼不同？轉錄組reads比對回基因組時需要特殊地考慮什麼？(4分)
- RNA測序並不能直接使用DNA測序常用的BWA、Bowtie等比對軟體，這是由於真核生物內含子的存在，導致測到的reads有時並不與基因組序列完全一致，因此需要使用Tophat2/HISAT/STAR等專門為RNA測序設計的軟體進行splice-map比對。
9、STAR為了提高比對率做了哪些容錯機制? (4分)
- STAR具有最高比例的在基因組上有唯一比對位置的reads，尤其是對讀長為300nt的樣品也有最高的比對率。STAR只保留雙端reads都比對到基因組的序列，但對低質量的比對(允許更多的錯配鹼基和soft-clip事件)容忍度高，這一點在長reads樣品中的體現更為明顯。
注：soft-clip事件即reads末端存在低質量鹼基或接頭導致比對不上的，STAR會自動嘗試截去未比對部分，只保留比對上的部分。
10、StringTie和Cufflinks是做什麼的？什麼時候會用到？(2分)
- 轉錄本拼裝， 研究新轉錄本，可變剪下時會用到。
11、Illumina測序時什麼情況下會測到接頭序列? (3分)
- 如果插入片段小於測序單端讀長，會在尾部測到接頭; 如果接頭之間發生了自連，則整條或大部分是接頭。
12、鏈特異文庫是什麼？可以解決什麼問題？(3分)
-
建庫時將mRNA鏈的方向資訊儲存到測序文庫中。最常用的dUTP，首先利用隨機引物合成RNA的一條cDNA鏈，在合成第二條鏈的時候用dUTP代替dTTP，加adaptor後用UDGase處理，將有U的第二條cDNA降解掉。
-
確定轉錄本來自正鏈還是負鏈，更準確的獲得基因的結構以及基因表達資訊，並且可以更好的發現新的轉錄本、反義轉錄本等。
13、Illumina測序時測序簇 (cluster)是如何生成的？生成測序簇的意義是什麼？(2分) (選答)
- 單鏈文庫與flowcell中的oligo結合，進行快速等溫橋式擴增，在小範圍內由同一條鏈克隆出來的很多條鏈聚整合簇（cluster）。
- 放大熒光訊號，方便成像檢測。
14、雙端測序的左端和右端reads分別測序的什麼，是否來源於同一片段? (2分)
左端和右端reads分別測的是同一條模板的兩端的序列。左端和右端reads對應下圖的read1和read2。
120分的轉錄組試題和答案
二、實戰題目
1、[email protected]:~/data$ 中 ct、ehbio、~/data、$ 各代表什麼？ (1分)
ct ehbio ~/data $
2、寫程式碼計算FASTQ序列中的reads數、鹼基數? (1分)
- 計算reads數：
echo "`zcat trt_N061011_1.fq.gz | wc -l` / (4*1000000)" | bc -l # wc -l: 計算行數 # bc -l: 計算器 (-l：浮點運算) # 除以4：fastq文件每4行為一條reads # 除以1000000：單位換算為million
- 計算鹼基數：
zcat trt_N061011_1.fq.gz | awk '{if(FNR%4==0) base+=length}END{print base/10^9,"G";}'
# %：取除以4後的餘數
# / 10^9：單位換算為G
3、寫程式碼用Fastqc評估測序質量。並解釋什麼樣的質量值是可以接受的？GC含量正常和異常的結果有什麼不同? 導致GC含量異常的原因是什麼？ (4分)
fastqc filename #filename：通常為樣品名.fq.gz
-
Reads數和reads長度符合預期，鹼基質量Q30>=75，GC含量45%-50%（特殊序列除外），接頭序列越少越好，GC含量正常時曲線圓滑，有單個突出的峰，不正常會有多個突起，原因是樣品本身鹼基不平衡或引入了其他物種基因或引物二聚體。
-
具體見 NGS基礎 - FASTQ格式解釋和質量評估 。
4、寫程式碼去除雙端測序reads低質量鹼基和接頭序列，並解釋各個引數的含義？(3分)
trimmomatic.sh PE -phred33 input_forward.fq input_reverse.fq \ output_forward_paired.fq \ output_forward_unpaired.fq \ output_reverse_paired.fq \ output_reverse_unpaired.fq \ ILLUMINACLIP:adaptor-PE.fa:2:30:10 LEADING:20 TRAILING:20 MINLEN:36 # PE：雙端測序 # -phred33：質量體系 # ILLUMINACLIP：接頭型別，去除接頭和引物序列 # LEADING:20 和 TRAILING：20：去除前後各20個低質量鹼基 # MINLEN:36：去除小於36bp的reads
5、為什麼一般只使用去除低質量鹼基和接頭後仍然成對的reads做後續分析? (2分)
- 為了結果更準確。
6、寫程式碼完成基因組索引的構建，並解釋基因組索引構建的意義是什麼? (2分)
hisat2-build -p 2 GRCh38.fa GRCh38
- 提高比對速度。
7、解釋環境變數是什麼？設定環境變數的意義是什麼？如何設定環境變數？(3分)
- 告訴系統在哪幾個目錄下有可執行程序，當輸入可執行程序名字時系統會去這幾個指定目錄中查詢對應文件是否存在並執行。設定環境變數是為了更方便呼叫可執行程序。（相當於windows中快捷方式）
export PATH=$PATH:絕對路徑
8、寫程式碼進行序列比對，並解釋每個引數的含義? (2分)
hisat2 -x genome/GRCh38 -1 trt_N061011_1.fq.gz -2 trt_N061011_2.fq.gz -p 2 \ --mm -S trt_N061011/trt_N061011.sam #-p 4: 使用4個執行緒。 #-x：索引文件的名字。(若不在當前目錄需包含路徑) #-1, -2：左端和右端序列。 #--mm: 多個hisat2共享記憶體中的hisat2基因組索引，減少記憶體使用 #-S: 輸出文件
9、如何從比對結果中獲取比對reads數和比對率？怎麼判斷比對率是否合適？ (2分)
- 從hisat2比對的輸出可以看到reads數和比對率，可以重定向到某個文件方便以後檢視。
- 比對率是否合適沒有統一標準，一般比對率80%以上。或者未比對上的
reads做一下blast，依然沒有匹配，也可以說明都比對環節無問題。如果比對到其它物種，考慮存在汙染。
10、解釋STAR或HTSeq是如何計算基因的表達量(readscount)的，程序是根據註釋文件中的哪些資訊來計算某個基因的reads數的，如果一個基因有多個轉錄本怎麼計算基因表達量 (3分)
htseq-count -f bam -r pos -a 10 -t exon -s no -i gene_id -m union \ trt_N061011/trt_N061011.sortP.bam genome/GRCh38.gtf \ >trt_N061011/trt_N061011.readsCount
- 程序是根據註釋文件中
gene_name和exon資訊來判斷reads歸屬。多個轉錄本會統一計算。
11、寫程式碼轉換BAM文件為樣品間可比的bigWig文件，並說明bigWig文件的用途? (2分)
bam2wig.py -i trt_N061011/trt_N061011.sortP.bam -s genome/GRCh38.chromsize \ -o trt_N061011/trt_N061011 -t 1000000000 -q 0
| 匯入基因組瀏覽器進行資料視覺化。[高通量資料分析必備 | 基因組瀏覽器使用介紹 - 1](https://mp.weixin.qq.com/s/_5zoRvjku5pb9qtwrvlWYw)。 |
12、寫程式碼評估比對質量，包括評估基因5’-3’測序均一度評估、reads在基因組標誌區域評估和飽和度評估，並解釋每個引數的意義? (3分)
- 5’-3’測序覆蓋均一性評估：
# RseQC軟體 geneBody_coverage2.py -i trt_N061011/trt_N061011.bw -r genome/GRCh38.model.gtf.bed12 -o \ trt_N061011/trt_N061011.geneBody_coverage # -i：輸入文件 # -r：輸入參考文件 # -o:指定輸出文件
- Reads在基因組標誌區域評估：
Read_distribution.py -i trt_N061011/trt_N061011.sortP.bam -r genome/GRCh38.gtf.bed12 \ >trt_N061011/trt_N061011.read_distrib.xls
- 測序飽和度評估：
RPKM_saturation.py -i trt_N061011/trt_N061011.sortP.bam -r genome/GRCh38.gtf.bed12 \ -s 5 -q 0 -o trt_N061011/trt_N061011.RPKM_saturation
13、用RSEM計算基因表達TPM值，並解釋TPM值、FPKM值、與DESeq2標準化後的值、原始reads count的區別? (4分)
- RPKM：Reads per Kilobase per Million.
- RPKM=（比對到某個基因上的reads數/比對到基因組上的總reads數）/基因長度（k）* 10^6
- FPKM意義與RPKM極為相近。二者區別僅在於Fragment與Read。PE測序中，通常一個Fragment產生兩條reads。
120分的轉錄組試題和答案
14、寫R程式碼讀入reads count文件、樣品分組資訊文件，並解釋引數的含義? (1分)
- 讀入Reads count文件:
Data <- read.table("ehbio_trans.Count_matrix.xls", header=T, row.names=1, com='', quote='',check.names=F, sep="\t")
- 讀入樣品分組文件:
sample <- read.table("sampleFile", header=T, row.names=1, com='', quote='', check.names=F, sep="\t", colClasses="factor")
15、寫程式碼產生DESeqDataSet資料集，並解釋引數的含義? (3分)
ddsFullCountTable <- DESeqDataSetFromMatrix(countData = data,colData = sample, design= ~ conditions) dds <- DESeq(ddsFullCountTable) #countData：表達矩陣 #cloData：樣品分組資訊 #design：實驗設計資訊 #conditions必須是cloData中的一列
16、解釋函式DESeq執行時函式內部都做哪些資料處理和操作? (5分)
- 表達量的標準化，差異基因分析。
17、寫程式碼從DESeqDataSet資料集獲取標準化後的資料、標準化後取對數的資料，並解釋用到的函式的作用? (3分)
- 標準化資料:
normalized_counts <- counts(dds, normalized=TRUE) normalized_counts_mad <- apply(normalized_counts, 1, mad) normalized_counts <- normalized_counts [order(normalized_counts_mad, decreasing=T), ]
- 資料輸出:
write.table(normalized_counts,file="ehbio_trans.Count_matrix.xls.DESeq2.normalized.xls",quote=F, sep="\t", row.names=T, col.names=T)
- 取對數:
rld <- rlog(dds, blind=FALSE) rlogMat <- assay(rld) rlogMat <- rlogMat[order(normalized_counts_mad, decreasing=T), ]
- 資料輸出:
write.table(rlogMat,file="ehbio_trans.Count_matrix.xls.DESeq2.normalized.rlog.xls",quote=F, sep="\t", row.names=T, col.names=T)
18、寫程式碼在DESeq中做主成分分析，並解釋主成分分析的結果? (2分)
plotPCA(rld, intgroup=c("conditions"), ntop=5000)
- 每個點代表一個樣品，點與點之間的距離越近，基因成分越相似。
19、寫程式碼計算樣品相關性並完成相關性熱圖繪製，解釋下述哪種矩陣適合做相關性分析:原始readscount、標準化後的矩陣、標準化後對數處理的矩陣？(4分)
- 生成顏色：
hmcol <- colorRampPalette(brewer.pal(9, "GnBu"))(100)
- 計算相關性pearson correlation：
pearson_cor <- as.matrix(cor(rlogMat, method="pearson"))
- 層級聚類：
hc <- hcluster(t(rlogMat), method="pearson")
- 熱圖繪製：
pdf("ehbio_trans.Count_matrix.xls.DESeq2.normalized.rlog.pearson.pdf", pointsize=10)
heatmap.2(pearson_cor, Rowv=as.dendrogram(hc), symm=T, trace="none",
col=hmcol, margins=c(11,11), main="The pearson correlation of each
sample")
dev.off()
計算樣品相關性用標準化後對數處理的矩陣。
20、寫程式碼提取兩組樣品中的差異基因分析結果，並解釋用到的函式和變數的意義? (3分)
sampleA = "untrt"
sampleB = "trt"
contrastV <- c("conditions", sampleA, sampleB)
res <- results(dds, contrast=contrastV)
baseA <- counts(dds, normalized=TRUE)[, colData(dds)$conditions == sampleA]
if (is.vector(baseA)){
baseMeanA <- as.data.frame(baseA)
} else {
baseMeanA <- as.data.frame(rowMeans(baseA))
}
colnames(baseMeanA) <- sampleA
baseB <- counts(dds, normalized=TRUE)[, colData(dds)$conditions == sampleB]
if (is.vector(baseB)){
baseMeanB <- as.data.frame(baseB)
} else {
baseMeanB <- as.data.frame(rowMeans(baseB))
}
colnames(baseMeanB) <- sampleB
res <- cbind(baseMeanA, baseMeanB, as.data.frame(res))
res <- cbind(ID=rownames(res), as.data.frame(res))
res$baseMean <- rowMeans(cbind(baseA, baseB))
res$padj[is.na(res$padj)] <- 1
res <- res[order(res$pvalue),]
comp314 <- paste(sampleA, "_vs_", sampleB, sep=".")
comp314 <- paste(sampleA, "_vs_", sampleB, sep=".")
file_base <- paste("ehbio_trans.Count_matrix.xls.DESeq2", comp314, sep=".")
file_base1 <- paste(file_base, "results.xls", sep=".")
write.table(as.data.frame(res), file=file_base1, sep="\t", quote=F, row.names=F)
使用函式和變數為了程序的方便和實用。
21、寫程式碼繪製火山圖，並解讀火山圖展示的資訊是什麼? (1分)
logCounts <- log2(res$baseMean+1) logFC <- res$log2FoldChange FDR <- res$padj \#png(filename=paste(file_base, "Volcano.png", sep=".")) plot(logFC, -1*log10(FDR), col=ifelse(FDR<=0.01, "red", "black"), xlab="logFC", ylab="-1*log1o(FDR)", main="Volcano plot", pch=".") \#dev.off()
左右兩邊代表下調和上調的基因，每個點是一個基因，點越高差異越大。 R語言 - 火山圖
22、寫程式碼繪製差異基因熱圖，並解讀此熱圖? (2分)
res_de_up_top20_id <- as.vector(head(res_de_up$ID,20)) res_de_dw_top20_id <- as.vector(head(res_de_dw$ID,20)) red_de_top20 <- c(res_de_up_top20_id, res_de_dw_top20_id) red_de_top20 # 可以把矩陣儲存，提交到www.ehbio.com/ImageGP 繪製火山圖 # 或者使用我們的s-plot https://github.com/Tong-Chen/s-plot red_de_top20_expr <- normalized_counts[rownames(normalized_counts) %in% red_de_top20,] library(pheatmap) pheatmap(red_de_top20_expr, cluster_row=T, scale="row", annotation_col=sample)
- 橫軸為樣品，縱軸為基因，紅色表達量越高，藍色表達量低。
23、繪製熱圖時什麼時候做行聚類、什麼時候做列聚類、按行標準化資料的意義是什麼、結果怎麼解讀? (2)
-
樣品之間差異比較做列聚類，基因之間差異比較做行聚類。
-
避免某個表達量異常高，影響中等表的量的顯示。把真正的表達量隱藏，以便顯示錶達量的趨勢。紅色表的量高，藍色表的量地；每行進行比較，發現每個基因在不同樣品中的表達情況。
24、如何識別樣品中是否有批次效應? (2分) (選答)
*在samplefile中新增一行批次列，從熱圖上看，差異基因應該是以分組資訊聚類的，但如果以批次列聚類了，就說明存在批次效應。
25、描述GO富集分析的原理，並給出GO富集分析需要的資料的格式? (1分)
-
採用超幾何檢驗判斷選定的差異基因在GO的哪些通路富集。
-
根據挑選出的差異基因，計算這些差異基因同GO分類中某（幾）個特定的分支的超幾何分佈關係，GO富集分析會對每個有差異基因存在的GO返回一個p-value，小的p值表示差異基因在該GO 中出現了富集。
-
資料格式：兩列文件，第一列基因ID，第二列是樣本差異。
26、描述GSEA富集分析的原理，並給出GSEA富集分析需要的資料的格式? (1分)
-
需要兩個文件，參考文件：排序好已知功能註釋後的基因集 分析的文件：自己做的表達矩陣。
-
用表達矩陣去與基因集比較，統計表達矩陣中的基因對應到基因集的位置是否集中。
27、GO富集分析結果的泡泡圖怎麼解釋? (1分)
-
基因比例的比率。
-
越紅越顯著，點越大差異基因數目越多。
28、轉錄本拼裝的意義是什麼? 不同樣品拼裝結果最後需要merge在一起的意義是? (2分)
-
發現新的轉錄本，新的可變剪下形式，鑑定非編碼RNA.
-
合併是為了新發現的方便。
29、DESeq2中是如何處理批次效應的? (2分) (選答)
- 用limma 包中的removeBatchEffect函式去除批次效應。把標準化後取對數後的矩陣和樣品分組表（標軸必須有批次列）交給函式，函式會根據批次列進行去除，也可以在DESeq2計算表達量時加入design=~batch+conditions引數，會自動去除批次效應。