一篇文獻一個技術路線：火熱的免疫

首先針對部分粉絲的需求統一回答下，弗雷賽斯目前隻有2種課題服務模式：

1.提供完整技術路線的中文分析報告（包括原始數據，圖片及解讀，用戶可自行翻譯後投稿）

2.深度合作，完成從選題到發表的全過程，未達到協議目标，弗雷賽斯将全額退款！

今天ZY乘勝追擊，繼續給大家講一篇文獻，這篇不新，但也足夠經典，因為講的是免疫，免疫絕對的火，凡事和免疫扯上關系，都能紅一波。生信自然不例外。生信結合免疫最經典的思路就是免疫分型！這篇發表在ONCOIMMUNOLOGY（IF=5.33）期刊詳情的文章就是基于免疫相關基因對乳腺癌樣本重新劃分具有不同免疫表型的病人，從多組學的角度刻畫這些不同免疫表型病人的特征。

研究背景：

免疫治療作為一種新的癌症治療方法，已經被用于幾種癌症的臨床治療。但是，免疫治療對乳腺癌的治療有限。由于缺乏對免疫反應基因組特征的認知，這導緻無法開發一種新的、有效的治療乳腺癌的方式。在本次研究中，WouterHendrickx等人為克服這個缺陷，他們整合TCGA數據庫中乳腺癌病人的拷貝數（CNV）、體細胞突變（somatic mutation）以及基因表達譜數據（RNA-seq），首先，基于1004個乳腺癌病人的RNA-seq表達譜數據，通過之前研究報道與免疫治療耐藥相關的轉錄本的表達将這些樣本劃分成具有不同免疫表型的樣本，即ICR1、ICR2、ICR3和ICR1，其中，在ICR1表型中這些和病人好的預後相關的轉錄本呈現了高表達（如：PDL1，PD1，FOXP3，IDO1以及CTLA4等）。

為說明這種表型刻畫的魯棒性，WouterHendrickx等人在另一套包含1954個乳腺癌病人的表達譜數據中去驗證他們的這種表型劃分标準，同時，WouterHendrickx等人又進一步結合拷貝數據去刻畫各表型病人特征，他們發現：（1）在ICR4表型病人中Th-1趨化因子CXCL9-11所位于的染色體片段（4q21）發生了顯著的擴增；（2）在ICR4-ICR1表型病人中，突變和新生抗原負荷的逐漸減少；（3）在ICR4表型病人中顯著發生TP53突變；（4）MAP3K1和MAP3K4的突變和ICR1表型病人相關；（5）在TCGA數據集以及驗證數據集中，可以從MAPK的失調的程度去劃分乳腺癌病人。他們的研究表明：在乳腺癌中，MAPK通路的突變驅動擾動與腫瘤内免疫反應負相關，通過多MAPK通路的調節可以增強乳腺癌免疫敏感性。

材料和方法

1.數據的獲取：

（1） 訓練集：通過R包TCGAbioloinks下載TCGA乳腺癌數據的RNA-seq、臨床數據和拷貝數SNP數據，TCGA portal中下載突變數據。首先，需剔除男性樣本、組織分期不明确的樣本、用過輔助化療藥物的樣本以及劃分為其他惡性腫瘤的樣本（剔除樣本數分别為：3, 9,1, 13, 66, 12）。

（2） 測試集：從GEO數據庫、NCI’s caArray數據庫以及EMBL-EBI’sArrayExpress數據庫獲取那些沒有接受輔助化療藥物的原發浸潤性乳腺癌芯片表達譜數據（n=1954），包括：raw cel (GSE1456,GSE2034, GSE5327, GSE12093, GSE7390, GSE6532, GSE9195, GSE2603, GSE7378,GSE8193, GSE4922, GSE11121,and GSE45255)、NCI’s caArray database (accession: mille-00271), andthe EMBL-EBI’s ArrayExpress database (accession: E-TABM-158)，這是數據的平台為：U133A, U133A2, andU133 PLUS 2.0 array platforms。

2. 數據的處理：

用R包EDASeq處理TCGA的RNA-seq表達譜數據；對芯片表達譜數據，首先先提取U133A,U133A2和U133PLUS 2.0平台共有的探針，然後再采用MAS5.0标準化，用COMBAT去批次。

3. 一緻性聚類：

将20個免疫排斥常數（immunologicconstant of rejection，ICR）基因和經典的免疫調控基因映射到TCGA乳腺癌樣本表達譜中，ConsensusClusterPlus包對該表達譜進行一緻性聚類，最大聚類數設置到7，根據Calinski指數（Calinski index）确定最佳聚類數，進一步刻畫這20個ICR基因和經典的免疫調控基因在相應類中的表達情況。

4. 乳腺癌病人中免疫細胞組成的刻畫：

基于Bindea G等人[1]的工作，得到24個免疫細胞基因集合，采用GSVA包，去計算每個樣本中對應的24個免疫細胞基因的ES得分，得到乳腺癌樣本的ES值所構成的矩陣，進一步通過層次聚類（聚類方法：Ward.D2），分析上述最佳聚類數下各免疫表型病人間的細胞類組成分布模式，并用gplots繪制熱圖可視化不同表型病人之間細胞組成的差異。

5. 免疫表型間差異基因的功能富集分析：

通過edgeR包計算所關注的免疫表型間的差異基因（(logFC) > 0.5 & p<0.05），基于這些顯著差異基因的log FC的絕對值排秩，利用GSEA軟件實現對這些差異基因功能的刻畫，并通過Cytoscape中的插件the Enrichment Maptool可視化。

6. 各免疫表型病人特征分析：

進一步分析在各表型病人中stage以及IMS亞型的分布特征，采用circos圖可視化這種分布。同時，結合臨床數據，通過ggplots2包中的ggkm函數比較各免疫表型病人的生存差異。

7. 突變負荷：

基于體細胞突變數據，用NetMHCpan包預測患者特異性HLA I類結合肽。同時，考慮到突變負荷高會産生較多的新抗原，病人的治療效果會更好，因此，進一步計算各免疫表型病人中非同義突變數的分布是否差異，從而對治療效果有進一步的評估。

8. 識别與免疫激活相關的突變基因：

基于Fisher’s精确檢驗（Fisher’s exacttest）識别表型間的差異基因（p <0.01），結合MUTSIG分析，從這些差異基因中進一步找出那些高于背景突變頻率的基因，将這些基因視為drivergene，并進一步在各表型中刻畫這些drivergene。

9. 拷貝數分析：

基于TCGA乳腺癌SNP 6.0的數據，采用GISTIC軟件去識别樣本中相應的擴增和缺失情況，其中，|log2 ratio| > 0.1是作為缺失和擴增的阈值，隻對滿足q value< 0.05的擴增和缺失拷貝數進行後續分析。采用NCBI’sGenome DecorationPage（GDP，http://www.ncbi.nlm.nih.gov/genome/tools/gdp）可視化染色體改變情況。

10. MAPK突變得分（MAPK-mut score）：

首先先将病人分成有MAP3K1/MAP2K4突變的Luminal病人組和無MAP3K1/MAP2K4突變的Luminal病人組，然後從KEGG數據庫得到MAPK通路中基因，接下來計算MAPK通路中的基因在這兩組病人中的z-score值的差異，得到在兩組病人的差異上調的基因和差異下調的基因。先對上下調基因在樣本中表達譜的行進行操作，對于每個基因，根據其在樣本中的z-score值進行排秩，将基因在樣本中的表達譜變成基因在每個樣本中表達的排秩情況表達譜；然後對列操作，對于每一個病人，計算基因的秩的平均值，該均值作為樣本的MAPK-mut得分。

結果

（1）各免疫表型的确定：

通過ConsensusClusterPlus包對20個ICR基因和經典的免疫調控基因映射到TCGA乳腺癌樣本表達譜進行一緻性聚類，通過評價Calinski index值（如圖Fig1所示），确定最佳的聚類數為4，劃分為：ICR1-ICR4，其中，ICR1 (N = 213)、ICR2 (N = 322)、ICR3 (N = 327)和 ICR4 (N = 142)（如圖Fig 2B所示），并發現這20個免疫相關基因在ICR4-ICR1表型病人中整體表達呈現下降趨勢（如圖Fig 2C所示）

Fig 1. 一緻性聚類Calinski index分布圖

Fig 2.免疫表型的确定

（2）免疫細胞組成的刻畫：基于Bindea G等人[1]整理的24個免疫細胞基因集合，通過GSVA包實現對各表型病人中24種免疫細胞組成的刻畫，通過熱圖的展示發現ICR4和ICR1兩個表型的細胞組成有着明顯的區别（如圖Fig3所示）。

Fig 3.ICR4和ICR1表型病人中24中免疫細胞組成的熱圖展示

（3）ICR1與ICR4表型病人中差異基因功能富集分析：

基Fig 3中ICR4和ICR1表型病人中24中免疫細胞組成差異的現象，分析兩個表型病人中的差異基因，并通過GSEA富集這些差異基因的功能，在ICR4表型中，适應免疫相關的通路和功能被激活，如白細胞分化、細胞因子産生、免疫效應過程、B- / T細胞活化以及T細胞增殖和代謝過程（如圖 Fig 4所示）。其中，節點表示相應的通路和功能，紅色表示在ICR4中高表達基因富集的功能，藍色表示在ICR4中低表達的基因所富集的功能，線的粗細表示兩個功能之間共享基因數目的多少。

Fig 4. GSEA富集分析

（4）免疫表型病人特征分析：

從在各表型病人中stag、IMS亞型的分布特征以及生存情況分析各表型病人的差異（如圖Fig 5所示）。

Fig 5. 各免疫表型特征分布.A)stage分布情況；B)IMS亞型分布情況；C)生存情況。

（5）突變負荷分析：

突變負荷高會産生較多的新抗原，病人的治療效果會更好，因此，結合TCGA病人的突變數據，統計各表型病人中的非同義突變數（如圖Fig 6所示），其中，各分類标準下突變數的差異p是通過Kruskal–Wallistest得到的。

Fig 6. 突變負荷展示.A)non-silent突變數在各不同樣本分類标準中的分布情況；B)在樣本中的突變率；C)ICR1-ICR4各表型病人突變情況統計。

（6）driver gene的識别：

通過Fisher’s exact test識别ICR1與ICR4表型間的差異基因，共得到64個差異基因（p < 0.01），結合MUTSIG分析，識别出5個顯著高于背景突變頻率的基因，即：MAP3K1、MAP2K4、TP53、RPGR和POM121，并通過熱圖、條形圖以及boxplot進行在各表型中突變情況進行展示（如圖Fig 7所示）。

Fig 7. Driver gene的識别和展示

（7）各表型病人拷貝數情況展示：

采用GISTIC軟件去識别樣本中相應的擴增和缺失情況，通過NCBI’sGenome DecorationPage（GDP，http://www.ncbi.nlm.nih.gov/genome/tools/gdp）對各免疫表型中病人拷貝數改變情況進行可視化（如圖Fig 8所示）。

Fig 8. 各表型拷貝數改變情況的展示

（8）：通過Fig 7A圖，發現MAPK的突變主要發生在luminal乳腺癌病人中，因此，通過整合突變數據和轉錄組學的數據，想進一步刻畫MAP3K1/MAP2K4突變所幹擾的功能，向在方法中描述的那樣去計算MAPK-mutation得分，發現根據ICR分類，該得分值可以很好的劃分luminal樣本（如圖Fig 9所示），在Fig 10中展示了MAPK通路的失調的通路圖，其中紅色的五角星表示在ICR4 vsICR1中上調的基因，綠色的五角星表示在ICR4 vsICR1中下調的基因。

Fig 9. MAPK-mutation得分刻畫不同免疫表型

Fig 10. MAPK通路失調示意圖

至此，我想大家已經對免疫分型的技術路線清楚了：從免疫基因篩選，到分型建立，結合臨床樣本特征，再到突變，拷貝數，進一步對某一個分型和某一個通路進行分析。如此，整篇文章從分子到臨床，從微觀到宏觀的分析層次讓工作顯得飽滿，影響因子自然不低。

參考文獻：

Bindea G, Mlecnik B, Tosolini M, et al.Spatiotemporal Dynamics of Intratumoral Immune Cells Reveal the ImmuneLandscape in Human Cancer[J]. Immunity, 2013, 39(4): 782-795.

......

生信學習需要持之以恒，大家可以到目前國内最大的生信學習社區（https://shengxin.ren/）逛逛，尤其是神秘的SangerBox可視化生信分析軟件，相信對大家玩轉生信一定有益。生信套路千千萬，不變的是生信的精髓:差異。生信及其他科研視頻，請直接點擊“閱讀全文”。