最近學徒剛剛完成了RNA-seq的訓練，就随機抽取了一個公共資料庫項目給他作為作業，是，研究者通過CRISPR/Cas9對 nasopharyngeal carcinoma (NPC) 細胞系的p53引入了一個雜合突變 R280T ，然後看mRNA profiles of CNE2 (control) and KO CNE2 cells 的差異。
可以看到是6個測序樣本
作者提供了表達矩陣：GSE130398_fpkm.txt.gz
第一個流程是 hisat2+featureCounts這個看我的視訊就好：
第二個流程是salmon
第三個流程在文章裡面，但是找不到文章，除非去詢問作者，Xiangya Hospital, Central South University 的 zhenqi qin
首先檢視突變是否引入成功
首先檢視bam檔案的頭：
@PG ID:hisat2 PN:hisat2 VN:2.1.0 CL:"/project/home/lyang/miniconda2/bin/hisat2-align-s --wrapper basic-0 -p 10 -x /project/home/lyang/refdata/hisat/human/hg38/genome -S SRR8980083.sam -1 /tmp/3190.inpipe1 -2 /tmp/3190.inpipe2"
複制
流程使用的是hg38參考基因組 , 簡單搜尋目标基因：https://www.ncbi.nlm.nih.gov/gene/7157
GRCh38.p13 (GCF_000001405.39) 17 NC_000017.11 (7668402..7687550, complement)
複制
也可以在gtf檔案搜尋：
chr17 HAVANA gene 7661779 7687550 . - . gene_id "ENSG00000141510.17"; gene_type "protein_coding"; gene_name "TP53"; level 2; hgnc_id "HGNC:11998"; tag "overlapping_locus"; havana_gene "OTTHUMG00000162125.10";
複制
批量提取 tp53 坐标的小bam：
ls *sort.bam | while read id;do (samtools view -b $id chr17:7661779-7687550 > ${id%%.*}_tp53.bam );done
ls *_tp53.bam |xargs -i samtools index {}
複制
這個 雜合突變 R280T ，需要找到坐标，首先需要知道轉錄本，有一個網頁工具可以繼續轉換，但是我忘記了。
假裝作者是對的，他們的實驗的确是引入了這個突變吧。本來都想發出去了，但是學徒憑運氣找到了這個位置，給大家過目：
然後看相關系數
三種檔案都準備好了：
首先看 salmon這樣的無需比對的流程結果和 hisat2+featureCounts的差異
可以看到，同一處理組的樣本在不同流程下面得到的表達量直接的相關性，是高于不同組的，符合邏輯！
但是單獨檢視同一個樣本的不同流程的表達量，如下所示：
可以看到，還是有不少基因在不同流程表現差異非常顯眼！那同樣的，我們需要檢查這些基因，簡單看看5個差異最大的基因吧。
在salmon裡面，可以看到第一個基因有3個轉錄本：
ENST00000445593.6_1 ENSG00000104341.16_2
ENST00000521545.6_1 ENSG00000104341.16_2
ENST00000517924.5_2 ENSG00000104341.16_2
複制
同樣的，salmon的這個樣本的結果如下：
Name Length EffectiveLength TPM NumReads
ENST00000445593.6 3173 2867.291 0.114792 4.192
ENST00000521545.6 915 609.515 0.000000 0.000
ENST00000517924.5 618 313.557 0.000000 0.000
複制
但是在hisat2+featureCounts的，6個樣本都有表達量。
ENSG00000104341.16 chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8 97775057;97775057;97775776;97775849;97805353;97805353;97805353;97805353;97815304;97815328;97815328;97815328;97816058;97816058;97816058;97816058;97819140;97819140;97819140;97819140;97825058;97825058;97825058;97851397;97851397;97851397 97776108;97776108;97776108;97776108;97805464;97805464;97805464;97805464;97815401;97815401;97815401;97815401;97816180;97816180;97816180;97816180;97819164;97819238;97819238;97819238;97825153;97825153;97825153;97851474;97853013;97852598 +;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+ 3197 5453 5150 5346 4250 4775 5271
複制
現在問題是如何判斷這個基因是否有表達量，還是需要知道坐标：
chr8 HAVANA gene 97775057 97853013
複制
使用指令去檢查bam檔案
samtools view SRR8980083.sort.bam chr8:97775057-97853013|wc
複制
發現這個基因的區域的确是有比對成功的reads，這就是我們所說的表達量。
第一個結論，是salmon會漏掉基因的真實表達量。
那麼salmon是否會憑空捏造基因的表達量，是以我們還需要檢查
chr3 HAVANA transcript 96617188 96618236 . - . gene_id "ENSG00000269028.3"; transcript_id "ENST00000600213.3"; gene_type "protein_coding"; gene_name "MTRNR2L12"; transcript_type "protein_coding"; transcript_name "MTRNR2L12-201"; level 2; protein_id "ENSP00000468991.1"; transcript_support_level "NA"; hgnc_id "HGNC:37169"; tag "not_best_in_genome_evidence"; tag "basic"; tag "appris_principal_1"; havana_gene "OTTHUMG00000175726.3"; havana_transcript "OTTHUMT00000430905.3";
複制
使用指令去檢查bam檔案
samtools view SRR8980083.sort.bam chr3:96617188-96618236 |wc
複制
的确是沒有表達量的，那麼為什麼salmon會捏造這個表達量呢？