生信技能樹有話說：最近剛剛組建一個TCGA交流群：4年前的TCGA重磅資料你學了嗎 把資料系統性整理給到了大家，發現其中一個群友聊天非常積極，而且看得到的進步，就要求他投稿，原題目是：
R菜鳥及生信菜鳥之第一篇Markdown筆記——歪打正著復現帖子的一個boxplot
作為一個剛入門生信的菜鳥，啥都不懂。
唯一有的優勢是積攢多年的生物學背景及學習能力快【囫圇吞棗能力】。再者，作為當代有痣青年，在這個生物大數據橫行，各種高通量技術橫行的時代，怎能缺少咱們的影子[捂臉]，還想著花點錢搞個單細胞測序發個CNS啥的，哈哈。有幸碰到了生信技能樹的團隊，剛好有這個機緣，所以乾脆就把腳踏入了這個圈。。。
走出舒適區，磨練自己。我的方法是邊看生信技能樹的教程，邊學技術，遇到不懂的，直接google，直接pubmed（這兩板斧最好用）。
話說，從接觸R到現在（28/11）到現在（03/12）總共不到一個星期裡，跟著依葫蘆畫瓢學習了幾個TCGA資料庫相關的教程。每個教程的code都是自己一個個手動打的，因為是一個個比較零散的教程，不成系統，所以基本大腦還沒有形成套路。我最想做的事情就是研究SCI文章裡面的數據挖掘套路，復現其數據處理過程及圖表產出。不僅可以跟上SCI文章的思路，而且可以學到技術。希望後面自己的生信體系會越來越成系統。
言歸正傳，按這個思路，我系統性研讀了生信技能樹系統的數據挖掘專題目錄：
其中 7. 腫瘤異質性+免疫浸潤細胞數據挖掘(可能是最簡單的3分⽂章了) 吸引了我的注意， 而且這個文章是TCGA資料庫的挖掘，所以應該重複性不錯。
說幹就幹，就臨摹起了這個帖子。一個下午大半時間各種折騰安裝相關的包。
同時泛讀了這個帖子相關的三篇文獻後：
開始幹活。Step1，Step2基本沒啥問題都處理好了。
在處理Boxplot的時候，帖子作者在復現KI67的時候用的是「從UCSC Xena下載gene expression RNAseq」的data，處理好了之後導入到R裡面的，然後發現出來的圖與原SCI不一致。
SCI原文圖：
image-20191204104113749帖子復現圖如下：
image-20191204104027232帖子處理的KI67 基因在兩組裡面沒有統計學差異。
我在這停頓了一下，帖子作者是從UCSC下載的數據，沒有提供下載【捂臉】！繼續不下去了，怎麼辦？只能自己操作獲取。想起前幾天學的cgdsr下載CGTA基因表達的數據，乾脆就寫了兩句自己獲取了，代碼如下：
library(cgdsr)
mycgds <- CGDS("http://www.cbioportal.org/")
mycaselist = 'brca_tcga_rna_seq_v2_mrna'
mygeneticprofile = 'brca_tcga_rna_seq_v2_mrna'
MKI67 = getProfileData(mycgds,'MKI67',mygeneticprofile,mycaselist)
MKI67$Sample_ID <- chartr('.','-',substr(row.names(MKI67),1,12))
df_ki67 <- merge(df,MKI67)
df_ki67$MKI67 <- log2(df_ki67$MKI67)
library(ggpubr)
p1 <- ggboxplot(df, x="group",y="Mutation_burder",
fill = "group",
palette = "jco",
legend = "") +
stat_compare_means(method = "t.test")
p2 <- ggboxplot(df_ki67, x="group",y="MKI67",
fill = "group",
palette = "jco",
legend = "") +
stat_compare_means(method = "t.test")
ggarrange(p1,p2,ncol = 2,nrow=1)
姨？發現我處理出來的圖基本吻合了SCI裡面的boxplot！如下：
image-20191204111112654總結：應該是數據處理方式不同。我沒有下載UCSC的數據，所以不知道裡面是如何進行normalization的。在cgdsr拿到KI67 的表達量之後，我這裡進行了log2處理以吻合原SCI的縱坐標。總之，歪打正著的趕腳。然後再測試了SCI文章Figure 3裡面的基因，發現也基本吻合，如下圖：
image-20191204111820690SCI原圖：
image-20191204111854986不過p值和SCI文章不一致，估計和數據處理和異常值有關。
而後一直在找intratumor heterogeneity的數據，原SCI作者沒說他的數據怎麼來的。我估計也是從ATCG上處理過來的，可是我還不會~~~
這個帖子還沒臨摹完，後續繼續做生存相關及免疫細胞亞群的分析。。。fighting