MSI（Microsatellite Instability）微衛(wèi)星不穩(wěn)定性是一種遺傳變異現(xiàn)象，指的是腫瘤細(xì)胞中微衛(wèi)星DNA重復(fù)序列的插入或缺失，導(dǎo)致某些DNA微衛(wèi)星的長度發(fā)生異常改變。MSI在多種癌癥中普遍存在，尤其是結(jié)直腸癌。

MSI微衛(wèi)星不穩(wěn)定性分析的常用方法是通過PCR擴(kuò)增腫瘤組織和正常組織中的一組微衛(wèi)星標(biāo)記位點(diǎn)，并比較它們的DNA序列長度。通常選擇了一組已知的高度多態(tài)性的微衛(wèi)星標(biāo)記位點(diǎn)，這些位點(diǎn)在正常組織中具有穩(wěn)定的DNA序列長度。如果在腫瘤組織中發(fā)現(xiàn)了與正常組織相比顯著的微衛(wèi)星長度改變（即微衛(wèi)星不穩(wěn)定性），則說明腫瘤具有MSI。

library(limma)

library(data.table)

library(tidyverse)

library(ggsignif)

library(RColorBrewer)

library(future.apply)

library(fmsb)

載入R包

files=dir("./TG_Exp/")

#Batch read geneexp files 批量讀取基因表達(dá)文件

#示例基因表達(dá)文件，這是之前根據(jù)tcga手動(dòng)批量整理的，后續(xù)有人想了解，可以聯(lián)系小云。

##下面是編寫的一個(gè)功能函數(shù)來批量讀取處理單基因表達(dá)文件，函數(shù)接受一個(gè)目錄名作為輸入a. 讀取位于./TG_Exp/目錄中指定的目錄名的基因表達(dá)文件。 b. 過濾數(shù)據(jù)，只保留腫瘤樣本。 c. 使用數(shù)據(jù)的第一列作為行名。 d. 去掉表達(dá)矩陣中的第一列，因?yàn)榈谝涣型ǔＪ菢颖久?，不是基因名?e. 解析目錄名，提取出腫瘤類型。 f. 讀取MSI（微衛(wèi)星不穩(wěn)定）數(shù)據(jù)，MSI數(shù)據(jù)的文件路徑為"./MSI.txt"。 g. 過濾MSI和基因表達(dá)數(shù)據(jù)，只保留兩者共有的樣本。

corStatMSI=function(dirname){

exp<-read.table(paste0("./TG_Exp/",dirname),header=T,sep='\t',check.names=F)

exp1=subset(exp,exp$sample_type=="tumor")

rownames(exp1)=exp1[,1]

exp1=exp1[,-1]

##names

tmp<- strsplit(dirname,split=".",fixed=TRUE)

dirnam<-unlist(lapply(tmp,head,1))

###read MSI讀取MSI微衛(wèi)星不穩(wěn)定數(shù)據(jù)，下面為數(shù)據(jù)示范。

MSI=read.table("./MSI.txt", header=T,sep="\t",row.names=1,check.names=F)

rownames(MSI)=substr(rownames(MSI),1,nchar(rownames(MSI))-1) ??#去最后一個(gè)字符A

#pass normal

group=sapply(strsplit(row.names(exp1),"\\-"),"[",4)

group=sapply(strsplit(group,""),"[",1)

group=gsub("2","1",group)

exp1=exp1[group==0,]

#intersect

exp1=na.omit(exp1)

sameSample=intersect(row.names(MSI),row.names(exp1))

MSI=MSI[sameSample,]

exp1=exp1[sameSample,]

#COR TEST

outTab=data.frame()

fmsbTab=data.frame()

tissue=c("ACC", "BLCA", "BRCA", "CESC", "CHOL", "COAD", "DLBC", "ESCA", "GBM", "HNSC", "KICH", "KIRC", "KIRP", "LGG", "LIHC", "LUAD", "LUSC", "MESO", "OV", "PAAD", "PCPG", "PRAD", "READ", "SARC", "SKCM", "STAD", "TGCT", "THCA", "THYM", "UCEC", "UCS", "UVM")

#對(duì)腫瘤類型循環(huán)

for ( j in 1:length(tissue) ){

?cancer=tissue[j]

exp2=exp1[(exp1[,2]==cancer),]

MSI1=MSI[(MSI[,2]==cancer),]

x=as.numeric(MSI1[,1])

y=as.numeric(exp2[,1])

if(length(x)!=length(y)){next}

if(length(unique(exp2[,1]))>3){

corT=cor.test(x,y,method="spearman")

cor=corT$estimate

pValue=corT$p.value

sig=ifelse(pValue<0.001,"***",ifelse(pValue<0.01,"**",ifelse(pValue<0.05,"*"," ")))

outTab=rbind(outTab,cbind(tissue=cancer,cor=cor,pValue=pValue,sig))

fmsbTab=rbind(fmsbTab,cbind(tissue=cancer,cor=cor))

}

else{next}

write.table(outTab,paste0("./corStat_MSI/",dirnam,"_","MSI_corStat.txt"),sep="\t",row.names=F,quote=F)

write.table(t(fmsbTab),paste0("./corStat_MSI/",dirnam,"_","MSI_fmsbInput.txt"),sep="\t",col.names=F,quote=F)

}

}

dd=future_lapply(files,corStatMSI) ??#多核并行運(yùn)行

后續(xù)可以給大家講解future_lapply多核并行的操作

下面為批量生成的結(jié)果MSI 文件txt，以其中一個(gè)單基因示范查看內(nèi)容，分別展示了相關(guān)性和p值。

#批量繪制雷達(dá)圖，讀取前面處理好的數(shù)據(jù)

首先，使用list.files()函數(shù)從指定文件夾中獲取符合條件的文件名列表。file1存儲(chǔ)了符合".fmsbInput.txt"模式的文件名列表，file2存儲(chǔ)了符合".corStat.txt"模式的文件名列表。

file1 <- list.files(paste0("./corStat_MSI"), pattern="*.fmsbInput.txt")

file2 <- list.files(paste0("./corStat_MSI"), pattern="*.corStat.txt")

for (i in 1:length(file1) ){

data=read.table(paste0("./corStat_MSI/",file1[i]),header=T,sep="\t",row.names=1,check.na mes=F) ??#read

a=file1[i]

a=gsub("|_MSI_fmsbInput\\.txt","",a)

接下來，通過循環(huán)遍歷file1列表中的每個(gè)文件名，讀取相應(yīng)的文件數(shù)據(jù)，并對(duì)數(shù)據(jù)進(jìn)行一些處理操作。data使用read.table()函數(shù)從文件中讀取數(shù)據(jù)，并將數(shù)據(jù)存儲(chǔ)在一個(gè)數(shù)據(jù)框中。

if (length(data)==0) {next}

maxValue=ceiling(max(abs(data))*10)/10

data=rbind(rep(maxValue,ncol(data)),rep(-maxValue,ncol(data)),data)

然后，通過對(duì)數(shù)據(jù)進(jìn)行處理，計(jì)算最大值maxValue，并在數(shù)據(jù)的頂部和底部插入相應(yīng)的最大值數(shù)據(jù)，以擴(kuò)展雷達(dá)圖的范圍。

##color

colors="blue"

#name +sig ***

corStat=read.table(paste0("./corStat_MSI/",a,"_MSI_corStat.txt"),header=T,sep="\t",row.names=1,check.names=F)

if (length(corStat)==0) {next}

colnames(data)=paste0(colnames(data),corStat$sig)

接下來，根據(jù)需求設(shè)置雷達(dá)圖的顏色colors。然后，使用read.table()函數(shù)讀取與當(dāng)前文件對(duì)應(yīng)的".corStat.txt"文件的數(shù)據(jù)，并存儲(chǔ)在corStat數(shù)據(jù)框中。然后，根據(jù)需求對(duì)數(shù)據(jù)框中列的名稱進(jìn)行處理，將其與corStat數(shù)據(jù)框中的sig列的值進(jìn)行拼接。接下來，根據(jù)文件名設(shè)置dirname和dirnam變量。

dirname=file1[i]

tmp<- strsplit(dirname,split="_",fixed=TRUE)

dirnam<-unlist(lapply(tmp,head,1))

#print

pdf(paste0("./RAD/",dirnam,"_MSI_radar.pdf"),height=7,width=7)

radarchart( data, axistype=1 ,

????pcol=colors, ????????????????#color

????plwd=2 , ????????????????????#lwd

????plty=1, ?????????????????????#dotted or solid line

????cglcol="grey", ??????????????#Background line color

????cglty=1, ?????????????????????

????caxislabels=seq(-maxValue,maxValue,maxValue/2), ???# scale

????cglwd=1.2, ??????????????????#Background line thickness

????axislabcol="green", ??????????#color

????vlcex=0.8 ???????????????????#size

)

dev.off()

}

來給大家展示展示結(jié)果圖吧，使用firefox命令。

那么這就是單基因在泛癌中的MSI微衛(wèi)星不穩(wěn)定性結(jié)果圖。

MSI微衛(wèi)星不穩(wěn)定性的分析對(duì)于癌癥的診斷和預(yù)后評(píng)估具有重要意義。MSI在一些腫瘤中與預(yù)后較好、對(duì)某些化療藥物敏感等特點(diǎn)相關(guān)。因此，對(duì)于結(jié)直腸癌等腫瘤，進(jìn)行MSI微衛(wèi)星不穩(wěn)定性分析已成為臨床的常規(guī)檢測(cè)之一。

下期將為你帶來更多R語言的騷操作技巧，以下推薦的是一個(gè)多功能的生信平臺(tái)。

云生信平臺(tái)鏈接：http://www.biocloudservice.com/home.html。

云生信平臺(tái)鏈接：http://www.biocloudservice.com/home.html。