此次給繪制的圖來自文獻(xiàn)《Triple DMARD treatment in early rheumatoid arthritis modulates synovial T cell activation and plasmablast/plasma cell differentiation pathways 》，是2017發(fā)表的，使用了他們團(tuán)隊(duì)自己2016發(fā)表的轉(zhuǎn)錄組測序數(shù)據(jù)，所以其實(shí)是有兩個(gè)轉(zhuǎn)錄組測序數(shù)據(jù)集。
Data are available publicly through NCBI GEO database (Healthy samples are under acces sion number GSE89408. RA samples are under accession number GSE97165)GSE89408 ：有28個(gè)正常樣本，GSE97165 有19個(gè) RA 治療前與 19個(gè)RA治療后的樣本。
這個(gè)圖主要展示了 A： 治療后 與 治療前的差異火山圖，B： 治療前 與正常對照 差異基因在三組樣本中的表達(dá)熱圖，以及 C&D ：一些 marker 基因在三個(gè)組別中的箱線圖+抖動(dòng)散點(diǎn)+顯著性比較 。
Fig 2. Genes selectively expressed by lymphocytes are decreased after 6 months of tDMARD treatment, while other genes remain elevated compared to normal joints 這里的兩個(gè)數(shù)據(jù)集中的count數(shù)據(jù)作者使用的同一套流程與參數(shù)進(jìn)行分析，這里直接合并在一起進(jìn)行后續(xù)分析。 （兩個(gè)數(shù)據(jù)集都是同一個(gè)科研團(tuán)隊(duì)的，所以理論上可以先不考慮里面的批次效應(yīng)哈）
如果實(shí)在不放心，可以選擇下載fq數(shù)據(jù)放在一起定量，拿到表達(dá)矩陣走后面的分析。
rm(list = ls())#清空當(dāng)前的工作環(huán)境 #不以因子變量讀取 #不以科學(xué)計(jì)數(shù)法顯示 ####################### 讀取 正常樣本：PRJNA352076、GSE89408 "GSE89408" "GSE89408_GEO_count_matrix_rename.txt.gz" , data.table = F)$V1 ), ]####################### RA樣本：PRJNA380832、GSE97165 "GSE97165" "GSE97165_GEO_count_matrix.txt.gz" , data.table = F)$V1 ), ]"V1" )"normal" ,colnames(mat)),grep("RA_pre" ,colnames(mat)),grep("RA_post" ,colnames(mat))) ]#過濾低表達(dá) "_[1-9]" , n = 2, simplify = T)[,1]"_tissue" , "" , group_list)"RA_pre" , "RApre" , group_list)"RA_post" , "RApost" , group_list)# group_list # normal RApre RApost # 28 19 19 #dat <- log10(edgeR::cpm(symbol_matrix)+1) 'step1-output.Rdata' )文獻(xiàn)中使用的是 limma 算法，我們也盡量復(fù)現(xiàn)同樣的哈，其中，疾病和對照肯定是差異巨大，但是治療前后就很難說了因?yàn)閺奈墨I(xiàn)里面的pca來看本來就是分組內(nèi)的差異并沒有顯著的小于組間差異！
rm(list = ls()) ## 魔幻操作，一鍵清空~ 'timeout' )'./step1-output.Rdata' ) 'RA_vs_Normal' in % c('normal' ,'RApre' );table(kp)'normal' ,'RApre' )) 'symbol_matrix.Rdata' )# 運(yùn)行 source ('../scripts/step3-deg-limma-voom.R' ) rm(list = ls()) ## 魔幻操作，一鍵清空~ 'timeout' )'./step1-output.Rdata' ) 'RApost_vs_RApre' in % c('RApre' ,'RApost' );table(kp)'RApre' ,'RApost' )) 'symbol_matrix.Rdata' )# 運(yùn)行 source ('../scripts/step3-deg-limma-voom.R' ) rm(list = ls()) 'symbol_matrix.Rdata' ) log =TRUE, prior.count=3)"quantile" )'-' ,g1)'DEG_limma_voom.Rdata' ) load( file = 'DEG_limma_voom.Rdata' )"logFC" , "P.Value" )]$g =ifelse(df$P .Value > pvalue_t,'stable' , #if 判斷：如果這一基因的P.Value>0.01，則為stable基因 $logFC > logFC_t,'up' , #接上句else 否則：接下來開始判斷那些P.Value<0.01的基因，再if 判斷：如果logFC >1.5,則為up（上調(diào)）基因 $logFC < -logFC_t,'down' ,'stable' ) )#接上句else 否則：接下來開始判斷那些logFC <1.5 的基因，再if 判斷：如果logFC <1.5，則為down（下調(diào)）基因，否則為stable基因 $g )$name =rownames(df)'Threshold of logFC is ' ,round(logFC_t,3),'\nThe number of up gene is ' ,nrow(df[df$g == 'up' ,]) ,'\nThe number of down gene is ' ,nrow(df[df$g == 'down' ,]) )#～～～火山圖p5～～～ 'MS4A1' )$name %in % target_gene) "ggrepel.max.overlaps" , default = 20),"black" ) +"-log10(Pvalue)" )+"#34bfb5" , "#828586" ,"#ff6633" ))+"grey" ,lwd=0.8) +"grey" ,lwd=0.8) +#coord_flip()+ #xlim(-5, 5)+ 結(jié)果如下：比文獻(xiàn)中那個(gè)火山圖好看，就用我們這個(gè)進(jìn)行后續(xù)的拼圖好了。
rm(list = ls()) ## 魔幻操作，一鍵清空~ 'timeout' )"RA_vs_Normal/" )# 讀取差異結(jié)果 "DEG_limma_voom.Rdata" )# 讀取所有樣本表達(dá)矩陣 "../step1-output.Rdata" )###################################### Fig2.B 熱圖 # 提取所有差異基因的表達(dá)矩陣 $name <- rownames(DEG_limma_voom)$g <- "stable" $g [ DEG_limma_voom$logFC >1 & DEG_limma_voom$adj .P.Val < 0.05 ] <- "up" $g [ DEG_limma_voom$logFC < -1 & DEG_limma_voom$adj .P.Val < 0.05 ] <- "down" $g )$g !="stable" , "name" ]"#ed1a22" ,RApre="#00a651" ,RApost="#652b90" ))'#1f67ad' ,'white' ,'#bb2f34' ))(100)"row" ,show_rownames = F,show_colnames = F,"black" ,gaps_col=c(28,47))結(jié)果如下：文獻(xiàn)中那個(gè)圖感覺只放了 上調(diào)基因？
首先，讓kimi（https://kimi./） 幫我拿到圖片中的基因并整理成一個(gè)R語言向量，再也不用一個(gè)個(gè)手動(dòng)從文獻(xiàn)里面摳出來了：
圖C 的比較組別為：list(c("RApre", "normal"), c("RApost", "RApre"))
###################################### Fig2.C 熱圖 # 提取所有差異基因的表達(dá)矩陣 "CD3D" , "MS4A1" , "CTLA4" , "CD19" )"normal" ,"RApre" ,"RApost" ))for (i in 1:length(genes1)) {#i <- 4 $exp )# 繪制小提琴圖和顯著性標(biāo)記 # 箱線圖 # 散點(diǎn) "normal" ,"RApre" ,"RApost" ), values =c( "#ed1a22" ,"#00a651" ,"#652b90" ) ) + # 顏色 # 不同組別的顯著性 "RApre" , "normal" ), c("RApost" , "RApre" )), # T顯示顯著性，F(xiàn)顯示p value # 修改顯著性線兩端的長短 # 設(shè)置顯著性線的位置高度 # 修改線的粗細(xì) # 修改顯著性標(biāo)記的大小 test = "t.test" ) + # 檢驗(yàn)的類型,可以更改 #設(shè)置白色背景 "" ,y="" ) + # 添加標(biāo)題，x軸，y軸標(biāo)簽 "black" ,size=0.9))"collect" )結(jié)果如下：
圖D的比較組別為：list(c("RApre", "normal"), c("RApost", "normal"))
########################################################################## "IL10" , "CLEC12A" , "AURKA" , "MMP13" ,"CLEC2B" , "CD58" )for (i in 1:length(genes2)) {#i <- 4 $exp )# 繪制小提琴圖和顯著性標(biāo)記 # 箱線圖 # 散點(diǎn) "normal" ,"RApre" ,"RApost" ), values =c( "#ed1a22" ,"#00a651" ,"#652b90" ) ) + # 顏色 # 不同組別的顯著性 "RApre" , "normal" ), c("RApost" , "normal" )), # T顯示顯著性，F(xiàn)顯示p value # 修改顯著性線兩端的長短 # 設(shè)置顯著性線的位置高度 # 修改線的粗細(xì) # 修改顯著性標(biāo)記的大小 test = "t.test" ) + # 檢驗(yàn)的類型,可以更改 #設(shè)置白色背景 "" ,y="" ) + # 添加標(biāo)題，x軸，y軸標(biāo)簽 "black" ,size=0.9))"collect" ，ncol=3)結(jié)果如下：
將以上峰圖 保存成pdf，用ai進(jìn)行處理，得到下面完整的fig2：
強(qiáng)烈建議你推薦給身邊的博士后以及年輕生物學(xué)PI ，多一點(diǎn)數(shù)據(jù)認(rèn)知，讓他們的科研上一個(gè)臺(tái)階：
