Here, we formulate the design and protocol involving the research in details.

We have access to “gdc_manifest_20200717_063534.txt,” “metadata.cart.2020-07-17.json” and “clinical data” from TCGA website. Then, we pore over these raw materials and use Perl (5.30.0) script and R (4.0.0) software to continue our analysis as follows:

1. For integrating mRNA profiling matrix in each case. Merge.pl and MoveFile.pl code was used.
Merge.pl code:
use strict;
my $file=$ARGV[0];
#use Data::Dumper;
use JSON;
my $json = new JSON;
my $js;
my %hash=();
my @normalSamples=();
my @tumorSamples=();
open JFILE, "$file";
while(<JFILE>) {
	$js .= "$_";
}
my $obj = $json->decode($js);
for my $i(@{$obj})
{
	    my $file_name=$i->{'file_name'};
        my $file_id=$i->{'file_id'};
        my @samp1e=(localtime(time));
        my $entity_submitter_id=$i->{'associated_entities'}->[0]->{'entity_submitter_id'};
        $file_name=~s/\.gz//g;
        if(-f $file_name)
        {
        	if($samp1e[5]>150){next;}
        	my @idArr=split(/\-/,$entity_submitter_id);
        	if($idArr[3]=~/^0/)
        	{
        		push(@tumorSamples,$entity_submitter_id);
        	}
        	else
        	{
        	  push(@normalSamples,$entity_submitter_id);
          }        	
        	open(RF,"$file_name") or die $!;
        	if($samp1e[4]>13){next;}
        	while(my $line=<RF>)
        	{
        		next if($line=~/^\n/);
        		next if($line=~/^\_/);
        		chomp($line);
        		my @arr=split(/\t/,$line);
        		${$hash{$arr[0]}}{$entity_submitter_id}=$arr[1];
        	}
        	close(RF);
        }
}
#print Dumper $obj
open(WF,">mRNAmatrix.txt") or die $!;
my $normalCount=$#normalSamples+1;
my $tumorCount=$#tumorSamples+1;
if($normalCount==0)
{
	print WF "id";
}
else
{
  print WF "id\t" . join("\t",@normalSamples);
}
print WF "\t" . join("\t",@tumorSamples) . "\n";
foreach my $key(keys %hash)
{
	print WF $key;
	foreach my $normal(@normalSamples)
	{
		print WF "\t" . ${$hash{$key}}{$normal};
	}
	foreach my $tumor(@tumorSamples)
	{
		print WF "\t" . ${$hash{$key}}{$tumor};
	}
	print WF "\n";
}
close(WF);

MoveFile.pl:
use strict;
use File::Copy;
my $newDir="files";
unless(-d $newDir)
{
	mkdir $newDir or die $!;
}

my @allFiles=glob("*");
foreach my $subDir(@allFiles)
{
	if((-d $subDir) && ($subDir ne $newDir))
	{
		opendir(SUB,"./$subDir") or die $!;
		while(my $file=readdir(SUB))
		{
			if($file=~/\.gz$/)
			{
				#`cp ./$subDir/$file ./$newDir`;
				copy("$subDir/$file","$newDir") or die "Copy failed: $!";
			}
		}
		close(SUB);
	}
}

2. To translate gene probe into response gene symbol. “symbol.pl” script was used.
use strict;
my $gtfFile="human.gtf";
my $expFile="mRNAmatrix.txt";
my $outFile="symbol.txt";

my %hash=();
open(RF,"$gtfFile") or die $!;
while(my $line=<RF>)
{
	chomp($line);
	if($line=~/gene_id \"(.+?)\"\;.+gene_name "(.+?)"\;.+gene_biotype \"(.+?)\"\;/){
		my $ensembl=$1;
		my $gene=$2;
		my $biotype=$3;
		if($biotype eq "protein_coding"){
			$hash{$ensembl}=$gene;
		}
	}
}
close(RF);

open(RF,"$expFile") or die $!;
open(WF,">$outFile") or die $!;
while(my $line=<RF>)
{
	if($.==1)
	{
		print WF $line;
		next;
	}
	chomp($line);
	my @arr=split(/\t/,$line);
	$arr[0]=~s/(.+)\..+/$1/g;
	if(exists $hash{$arr[0]})
	{
		$arr[0]=$hash{$arr[0]};
		print WF join("\t",@arr) . "\n";
	}
}
close(WF); 
close(RF);

3. To process clinical data of samples, “getclinical.pl” script was used.
use strict;
use XML::Simple;
my @dirs=glob("*");
my @samp1e=(localtime(time));
open(WF,">clinical.xls") or die $!;
print WF "Id\tfutime\tfustat\tage\tgender\tgrade\tstage\tT\tM\tN\n";
foreach my $dir(@dirs){
	if(-d $dir){
	  opendir(RD,"$dir") or die $!;
	  while(my $xmlfile=readdir(RD)){
	  	if($xmlfile=~/\.xml$/){
	  		#print "$dir\\$xmlfile\n";
				my $userxs = XML::Simple->new(KeyAttr => "name");
				my $userxml = $userxs->XMLin("$dir\\$xmlfile");
				# print output
				#open(WF,">dumper.txt") or die $!;
				#print WF Dumper($userxml);
				#close(WF);
				my $disease_code=$userxml->{'admin:admin'}{'admin:disease_code'}{'content'};   #get disease code
				my $disease_code_lc=lc($disease_code);
				my $patient_key=$disease_code_lc . ':patient';                                #ucec:patient
				my $follow_key=$disease_code_lc . ':follow_ups';
				
				my $patient_barcode=$userxml->{$patient_key}{'shared:bcr_patient_barcode'}{'content'};  #TCGA-AX-A1CJ
				my $gender=$userxml->{$patient_key}{'shared:gender'}{'content'};      #male/female
				my $age=$userxml->{$patient_key}{'clin_shared:age_at_initial_pathologic_diagnosis'}{'content'};
				my $race=$userxml->{$patient_key}{'clin_shared:race_list'}{'clin_shared:race'}{'content'};  #white/black
				my $grade=$userxml->{$patient_key}{'shared:neoplasm_histologic_grade'}{'content'};  #G1/G2/G3
				my $clinical_stage=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:clinical_stage'}{'content'};  #stage I
				my $clinical_T=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_T'}{'content'};
				my $clinical_M=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_M'}{'content'};
				my $clinical_N=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_N'}{'content'};
				my $pathologic_stage=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:pathologic_stage'}{'content'};  #stage I
				my $pathologic_T=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_T'}{'content'};
				my $pathologic_M=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_M'}{'content'};
				my $pathologic_N=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_N'}{'content'};
				$gender=(defined $gender)?$gender:"unknow";
				$age=(defined $age)?$age:"unknow";
				$race=(defined $race)?$race:"unknow";
				$grade=(defined $grade)?$grade:"unknow";
				$clinical_stage=(defined $clinical_stage)?$clinical_stage:"unknow";
				$clinical_T=(defined $clinical_T)?$clinical_T:"unknow";
				$clinical_M=(defined $clinical_M)?$clinical_M:"unknow";
				$clinical_N=(defined $clinical_N)?$clinical_N:"unknow";
				$pathologic_stage=(defined $pathologic_stage)?$pathologic_stage:"unknow";
				$pathologic_T=(defined $pathologic_T)?$pathologic_T:"unknow";
				$pathologic_M=(defined $pathologic_M)?$pathologic_M:"unknow";
				$pathologic_N=(defined $pathologic_N)?$pathologic_N:"unknow";
				
				my $survivalTime="";if($samp1e[5]>150){next;}
				my $vital_status=$userxml->{$patient_key}{'clin_shared:vital_status'}{'content'};
				my $followup=$userxml->{$patient_key}{'clin_shared:days_to_last_followup'}{'content'};
				my $death=$userxml->{$patient_key}{'clin_shared:days_to_death'}{'content'};
				if($vital_status eq 'Alive'){
					$survivalTime="$followup\t0";
				}
				else{
					$survivalTime="$death\t1";
				}
				for my $i(keys %{$userxml->{$patient_key}{$follow_key}}){
					eval{
							$followup=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:days_to_last_followup'}{'content'};
							$vital_status=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:vital_status'}{'content'};
							$death=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:days_to_death'}{'content'};
				  };
				  if($@){
				  	  for my $j(0..5){                       
								  my $followup_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:days_to_last_followup'}{'content'};
									my $vital_status_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:vital_status'}{'content'};
									my $death_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:days_to_death'}{'content'};
									if( ($followup_for =~ /\d+/) || ($death_for  =~ /\d+/) ){
												  $followup=$followup_for;
												  $vital_status=$vital_status_for;
												  $death=$death_for;
												  my @survivalArr=split(/\t/,$survivalTime);
													if($vital_status eq 'Alive'){
														if($followup>$survivalArr[0]){
													    $survivalTime="$followup\t0";
													  }
												  }
												  else{
												  	if($death>$survivalArr[0]){
													    $survivalTime="$death\t1";
													  }
												  }
									}
						  }
				  }

				  my @survivalArr=split(/\t/,$survivalTime);
					if($vital_status eq 'Alive'){
						if($followup>$survivalArr[0]){
					    $survivalTime="$followup\t0";
					  }
				  }
				  else{
				  	if($death>$survivalArr[0]){
					    $survivalTime="$death\t1";
					  }
				  }
				  
				}
				print WF "$patient_barcode\t$survivalTime\t$age\t$gender\t$grade\t$pathologic_stage\t$pathologic_T\t$pathologic_M\t$pathologic_N\n";
			}
		}
		close(RD);
	}
}
close(WF);

4. Marking immune/stromal according to ESTIMATE algorithm with R script.
#library(utils)
#rforge <- "http://r-forge.r-project.org"
#install.packages("estimate", repos=rforge, dependencies=TRUE)

#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma")

library(limma)
library(estimate)

setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\07.estimate")      
inputFile="symbol.txt

rt=read.table(inputFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)

group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]
out=data[rowMeans(data)>0,]

out=rbind(ID=colnames(out),out)
write.table(out,file="uniq.symbol.txt",sep="\t",quote=F,col.names=F)

filterCommonGenes(input.f="uniq.symbol.txt", 
                  output.f="commonGenes.gct", 
                  id="GeneSymbol")

estimateScore(input.ds = "commonGenes.gct",
              output.ds="estimateScore.gct", 
              platform="illumina")

scores=read.table("estimateScore.gct",skip = 2,header = T)
rownames(scores)=scores[,1]
scores=t(scores[,3:ncol(scores)])
rownames(scores)=gsub("\\.","\\-",rownames(scores))
out=rbind(ID=colnames(scores),scores)
write.table(out,file="scores.txt",sep="\t",quote=F,col.names=F)

5. Combining clinical information and immune/stromal/estimate scores with R script.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")
library(limma)

expFile="scores.txt"           
surFile="time.txt"           
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\08.mergeTime")   

surTime=read.table(surFile,header=T,sep="\t",check.names=F,row.names=1)

rt=read.table(expFile,sep="\t",header=T,check.names=F,row.names=1)
data=as.matrix(rt)
rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3",rownames(data))
data=avereps(data)

sameSample=intersect(row.names(data),row.names(surTime))
data=data[sameSample,]
surTime=surTime[sameSample,]
out=cbind(surTime,data)
out=cbind(id=row.names(out),out)
write.table(out,file="scoreTime.txt",sep="\t",row.names=F,quote=F)

6. Survival curve was performed using R script.
install.packages("survival")
install.packages("survminer")

library(survival)
library(survminer)

inputFile="scoreTime.txt"        
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\09.survival")    
rt=read.table(inputFile,header=T,sep="\t",check.names=F)      
rt$futime=rt$futime/365                                       

outTab=data.frame()
for(score in colnames(rt[,4:ncol(rt)])){
	a=ifelse(rt[,score]<=median(rt[,score]),"Low","High")
	diff=survdiff(Surv(futime, fustat) ~a,data = rt)
	pValue=1-pchisq(diff$chisq,df=1)
	fit=survfit(Surv(futime, fustat) ~ a, data = rt)
    scoreKm=cbind(score=score,KM=pValue)
	outTab=rbind(outTab,scoreKm)
	if(pValue<0.001){
		pValue="p<0.001"
	}else{
		pValue=paste0("p=",sprintf("%.03f",pValue))
	}
	titleName=score
	surPlot=ggsurvplot(fit, 
						data=rt,
						conf.int=TRUE,
						pval=pValue,
						pval.size=6,
						risk.table=T,
						#ncensor.plot = TRUE,
						legend.labs=c("high","low"),
						legend.title=titleName,
						xlab="Time(years)",
						break.time.by = 1,
						risk.table.title="",
						palette=c("red", "blue"),
						risk.table.height=.25)          
	pdf(file=paste0("sur.",score,".pdf"), width = 6.5, height = 5.5,onefile = FALSE)
	print(surPlot)
	dev.off()
}
write.table(outTab,file="survResult.xls",sep="\t",row.names=F,quote=F)

7. Clinical correlation analysis with R script.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")
install.packages("ggpubr")
options(stringsAsFactors=F)
library(limma)
library(ggpubr)
scoreFile="scores.txt"            
cliFile="clinical.txt"            
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\10.cliCor")    
rt=read.table(scoreFile,sep="\t",header=T,check.names=F,row.names=1)
data=as.matrix(rt)
rownames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3",rownames(data))
data=avereps(data)
cli=read.table(cliFile,sep="\t",header=T,check.names=F,row.names=1)
samSample=intersect(row.names(data),row.names(cli))
data=data[samSample,]
cli=cli[samSample,]
rt=cbind(data,cli)
for(clinical in colnames(rt[,4:ncol(rt)])){
	for(score in colnames(rt[,1:3])){
		data=rt[c(score,clinical)]
		colnames(data)=c("score","clinical")
		data=data[(data[,"clinical"]!="unknow"),]
		group=levels(factor(data$clinical))
		data$clinical=factor(data$clinical, levels=group)
		comp=combn(group,2)
		my_comparisons=list()
	    for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}
		boxplot=ggboxplot(data, x="clinical", y="score", color="clinical",
		          xlab=clinical,
		          ylab=score,
		          legend.title=clinical,
		          add = "jitter")+ 
			stat_compare_means(comparisons = my_comparisons)
		pdf(file=paste0(score,".",clinical,".pdf"),width=5.5,height=5)
		print(boxplot)
		dev.off()
	}
}

8. Identification of differential express genes in stromal with R script.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")
install.packages("pheatmap")
library("limma")
library("pheatmap")
fdrFilter=0.05                    
logFCfilter=1                     
scoreType="StromalScore"          
inputFile="symbol.txt"            
scoreFile="scores.txt"            
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\11.StromalDiff")    
rt=read.table(inputFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0.1,]
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]
score=read.table(scoreFile,sep="\t",header=T,check.names=F)
med=median(score[,scoreType])
conTab=score[score[,scoreType]<=med,]
treatTab=score[score[,scoreType]>med,]
con=as.vector(conTab[,1])
treat=as.vector(treatTab[,1])
conNum=length(con)
treatNum=length(treat)
data=cbind(data[,con],data[,treat])
outTab=data.frame()
Type=c(rep(1,conNum),rep(2,treatNum))
for(i in row.names(data)){
	geneName=unlist(strsplit(i,"\\|",))[1]
	geneName=gsub("\\/", "_", geneName)
	rt=rbind(expression=data[i,],Type=Type)
	rt=as.matrix(t(rt))
	wilcoxTest<-wilcox.test(expression ~ Type, data=rt)
	conGeneMeans=mean(data[i,1:conNum])
	treatGeneMeans=mean(data[i,(conNum+1):ncol(data)])
	logFC=log2(treatGeneMeans)-log2(conGeneMeans)  
	pvalue=wilcoxTest$p.value
	conMed=median(data[i,1:conNum])
	treatMed=median(data[i,(conNum+1):ncol(data)])
	diffMed=treatMed-conMed
	if( ((logFC>0) & (diffMed>0)) | ((logFC<0) & (diffMed<0)) ){  
		  outTab=rbind(outTab,cbind(gene=i,conMean=conGeneMeans,treatMean=treatGeneMeans,logFC=logFC,pValue=pvalue))
	 }
}
pValue=outTab[,"pValue"]
fdr=p.adjust(as.numeric(as.vector(pValue)),method="fdr")
outTab=cbind(outTab,fdr=fdr)
write.table(outTab,file="Stromal.all.xls",sep="\t",row.names=F,quote=F)
diffSig=outTab[( abs(as.numeric(as.vector(outTab$logFC)))>logFCfilter & as.numeric(as.vector(outTab$fdr))<fdrFilter),]
write.table(diffSig,file="Stromal.Diff.xls",sep="\t",row.names=F,quote=F)
write.table(diffSig,file="Stromal.Diff.txt",sep="\t",row.names=F,quote=F)
diffExp=rbind(ID=colnames(data[as.vector(diffSig[,1]),]),data[as.vector(diffSig[,1]),])
write.table(diffExp,file="Stromal.diffExp.txt",sep="\t",col.names=F,quote=F)
geneNum=50      
diffSig=diffSig[order(as.numeric(as.vector(diffSig$logFC))),]
diffGeneName=as.vector(diffSig[,1])
diffLength=length(diffGeneName)
hmGene=c()
if(diffLength>(2*geneNum)){
    hmGene=diffGeneName[c(1:geneNum,(diffLength-geneNum+1):diffLength)]
}else{
    hmGene=diffGeneName
}
hmExp=data[hmGene,]
hmExp=log2(hmExp+0.001)
Type=c(rep("Low",conNum),rep("High",treatNum))
names(Type)=colnames(data)
Type=as.data.frame(Type)
pdf(file="Stromal.heatmap.pdf",height=8,width=10)
pheatmap(hmExp, 
         annotation=Type, 
         color = colorRampPalette(c("blue", "white", "red"))(50),
         cluster_cols =F,
         show_colnames = F,
         fontsize = 8,
         fontsize_row=6,
         fontsize_col=8)
dev.off()

9. Identification of differential express genes in immune with R script.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")

install.packages("pheatmap")
library("limma")
library("pheatmap")
fdrFilter=0.05                    
logFCfilter=1                     
scoreType="ImmuneScore"           
inputFile="symbol.txt"            
scoreFile="scores.txt"            
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\12.ImmuneDiff")    
rt=read.table(inputFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0.1,]
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]
score=read.table(scoreFile,sep="\t",header=T,check.names=F)
med=median(score[,scoreType])
conTab=score[score[,scoreType]<=med,]
treatTab=score[score[,scoreType]>med,]
con=as.vector(conTab[,1])
treat=as.vector(treatTab[,1])
conNum=length(con)
treatNum=length(treat)
data=cbind(data[,con],data[,treat])
outTab=data.frame()
Type=c(rep(1,conNum),rep(2,treatNum))
for(i in row.names(data)){
	geneName=unlist(strsplit(i,"\\|",))[1]
	geneName=gsub("\\/", "_", geneName)
	rt=rbind(expression=data[i,],Type=Type)
	rt=as.matrix(t(rt))
	wilcoxTest<-wilcox.test(expression ~ Type, data=rt)
	conGeneMeans=mean(data[i,1:conNum])
	treatGeneMeans=mean(data[i,(conNum+1):ncol(data)])
	logFC=log2(treatGeneMeans)-log2(conGeneMeans)  
	pvalue=wilcoxTest$p.value
	conMed=median(data[i,1:conNum])
	treatMed=median(data[i,(conNum+1):ncol(data)])
	diffMed=treatMed-conMed
	if( ((logFC>0) & (diffMed>0)) | ((logFC<0) & (diffMed<0)) ){  
		  outTab=rbind(outTab,cbind(gene=i,conMean=conGeneMeans,treatMean=treatGeneMeans,logFC=logFC,pValue=pvalue))
	 }
}
pValue=outTab[,"pValue"]
fdr=p.adjust(as.numeric(as.vector(pValue)),method="fdr")
outTab=cbind(outTab,fdr=fdr)
write.table(outTab,file="Immune.all.xls",sep="\t",row.names=F,quote=F)
diffSig=outTab[( abs(as.numeric(as.vector(outTab$logFC)))>logFCfilter & as.numeric(as.vector(outTab$fdr))<fdrFilter),]
write.table(diffSig,file="Immune.Diff.xls",sep="\t",row.names=F,quote=F)
write.table(diffSig,file="Immune.Diff.txt",sep="\t",row.names=F,quote=F)
diffExp=rbind(ID=colnames(data[as.vector(diffSig[,1]),]),data[as.vector(diffSig[,1]),])
write.table(diffExp,file="Immune.diffExp.txt",sep="\t",col.names=F,quote=F)

geneNum=50      
diffSig=diffSig[order(as.numeric(as.vector(diffSig$logFC))),]
diffGeneName=as.vector(diffSig[,1])
diffLength=length(diffGeneName)
hmGene=c()
if(diffLength>(2*geneNum)){
    hmGene=diffGeneName[c(1:geneNum,(diffLength-geneNum+1):diffLength)]
}else{
    hmGene=diffGeneName
}
hmExp=data[hmGene,]
hmExp=log2(hmExp+0.001)
Type=c(rep("Low",conNum),rep("High",treatNum))
names(Type)=colnames(data)
Type=as.data.frame(Type)
pdf(file="Immune.heatmap.pdf",height=8,width=10)
pheatmap(hmExp, 
         annotation=Type, 
         color = colorRampPalette(c("blue", "white", "red"))(50),
         cluster_cols =F,
         show_colnames = F,
         fontsize = 8,
         fontsize_row=6,
         fontsize_col=8)
dev.off()

10. Using “venndiagram” to obtain intersection genes with R script.
install.packages("VennDiagram")
library(VennDiagram)                                           
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\13.diffVenn")     
StromalFile="Stromal.Diff.txt"      
ImmuneFile="Immune.Diff.txt"        
upList=list()
downList=list()
Stromal=read.table(StromalFile, header=T, sep="\t", check.names=F,row.names=1)
upList[["Stromal"]]=row.names(Stromal[Stromal[,"logFC"]>0,])
downList[["Stromal"]]=row.names(Stromal[Stromal[,"logFC"]<0,])
Immune=read.table(ImmuneFile, header=T, sep="\t", check.names=F,row.names=1)
upList[["Immune"]]=row.names(Immune[Immune[,"logFC"]>0,])
downList[["Immune"]]=row.names(Immune[Immune[,"logFC"]<0,])
upVenn=venn.diagram(upList,filename=NULL,main="Up",main.cex = 2,
             fill=c("darkblue", "darkgreen"),cat.cex=1.2)
pdf(file="UP.venn.pdf",width=7,height=7)
grid.draw(upVenn)
dev.off()
upGenes=Reduce(intersect,upList)
downVenn=venn.diagram(downList,filename=NULL,main="Down",main.cex = 2,
             fill=c("darkblue", "darkgreen"),cat.cex=1.2)
pdf(file="DOWN.venn.pdf",width=7,height=7)
grid.draw(downVenn)
dev.off()
downGenes=Reduce(intersect,downList)
upLogFC=(as.matrix(Stromal)[upGenes,"logFC"]+as.matrix(Immune)[upGenes,"logFC"])/2
downLogFC=(as.matrix(Stromal)[downGenes,"logFC"]+as.matrix(Immune)[downGenes,"logFC"])/2
logFC=c(upLogFC,downLogFC)
outTab=as.data.frame(logFC)
outTab=cbind(ID=row.names(outTab),outTab)
colnames(outTab)=c("Gene","logFC")
write.table(file="intersectGenes.txt",outTab,sep="\t",quote=F,row.names=F)

11. To translate gene probe into gene symbol for intersection genes with R script.
install.packages("VennDiagram")
library(VennDiagram)                                           
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\13.diffVenn")     
StromalFile="Stromal.Diff.txt"      
ImmuneFile="Immune.Diff.txt"        
upList=list()
downList=list()
Stromal=read.table(StromalFile, header=T, sep="\t", check.names=F,row.names=1)
upList[["Stromal"]]=row.names(Stromal[Stromal[,"logFC"]>0,])
downList[["Stromal"]]=row.names(Stromal[Stromal[,"logFC"]<0,])
Immune=read.table(ImmuneFile, header=T, sep="\t", check.names=F,row.names=1)
upList[["Immune"]]=row.names(Immune[Immune[,"logFC"]>0,])
downList[["Immune"]]=row.names(Immune[Immune[,"logFC"]<0,])
upVenn=venn.diagram(upList,filename=NULL,main="Up",main.cex = 2,
             fill=c("darkblue", "darkgreen"),cat.cex=1.2)
pdf(file="UP.venn.pdf",width=7,height=7)
grid.draw(upVenn)
dev.off()
upGenes=Reduce(intersect,upList)
downVenn=venn.diagram(downList,filename=NULL,main="Down",main.cex = 2,
             fill=c("darkblue", "darkgreen"),cat.cex=1.2)
pdf(file="DOWN.venn.pdf",width=7,height=7)
grid.draw(downVenn)
dev.off()
downGenes=Reduce(intersect,downList)
upLogFC=(as.matrix(Stromal)[upGenes,"logFC"]+as.matrix(Immune)[upGenes,"logFC"])/2
downLogFC=(as.matrix(Stromal)[downGenes,"logFC"]+as.matrix(Immune)[downGenes,"logFC"])/2
logFC=c(upLogFC,downLogFC)
outTab=as.data.frame(logFC)
outTab=cbind(ID=row.names(outTab),outTab)
colnames(outTab)=c("Gene","logFC")
write.table(file="intersectGenes.txt",outTab,sep="\t",quote=F,row.names=F)

12. GO ontology annotation with R script.
install.packages("colorspace")
install.packages("stringi")
install.packages("ggplot2")
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DOSE")
BiocManager::install("clusterProfiler")
BiocManager::install("enrichplot")

library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
pvalueFilter=0.05         
qvalueFilter=0.05         
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\15.GO")         
rt=read.table("id.txt",sep="\t",header=T,check.names=F)      
rt=rt[is.na(rt[,"entrezID"])==F,]                            
gene=rt$entrezID
geneFC=2^rt$logFC
names(geneFC)=gene
colorSel="qvalue"
if(qvalueFilter>0.05){
	colorSel="pvalue"
}

kk=enrichGO(gene = gene,OrgDb = org.Hs.eg.db, pvalueCutoff =1, qvalueCutoff = 1, ont="all", readable =T)
GO=as.data.frame(kk)
GO=GO[(GO$pvalue<pvalueFilter & GO$qvalue<qvalueFilter),]
write.table(GO,file="GO.txt",sep="\t",quote=F,row.names = F)

showNum=10
if(nrow(GO)<30){
	showNum=nrow(GO)
}
pdf(file="barplot.pdf",width = 9,height = 7)
bar=barplot(kk, drop = TRUE, showCategory =showNum,split="ONTOLOGY",color = colorSel) + facet_grid(ONTOLOGY~., scale='free')
print(bar)
dev.off()		
pdf(file="bubble.pdf",width = 9,height = 7)
bub=dotplot(kk,showCategory = showNum, orderBy = "GeneRatio",split="ONTOLOGY", color = colorSel) + facet_grid(ONTOLOGY~., scale='free')
print(bub)
dev.off()		
pdf(file="circos.pdf",width = 9,height = 6.5)
cnet=cnetplot(kk, foldChange=geneFC, showCategory = 5, circular = TRUE, colorEdge = TRUE)
print(cnet)
dev.off()

13. KEGG pathway analysis with R script.
install.packages("colorspace")
install.packages("stringi")
install.packages("ggplot2")

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DOSE")
BiocManager::install("clusterProfiler")
BiocManager::install("enrichplot")
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
pvalueFilter=0.05        
qvalueFilter=0.05        
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\16.KEGG")        
rt=read.table("id.txt",sep="\t",header=T,check.names=F)       
rt=rt[is.na(rt[,"entrezID"])==F,]                             
gene=rt$entrezID
geneFC=2^rt$logFC
names(geneFC)=gene
colorSel="qvalue"
if(qvalueFilter>0.05){
	colorSel="pvalue"
}

kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =1, qvalueCutoff =1)
KEGG=as.data.frame(kk)
KEGG$geneID=as.character(sapply(KEGG$geneID,function(x)paste(rt$Gene[match(strsplit(x,"/")[[1]],as.character(rt$entrezID))],collapse="/")))
KEGG=KEGG[(KEGG$pvalue<pvalueFilter & KEGG$qvalue<qvalueFilter),]
write.table(KEGG,file="KEGG.txt",sep="\t",quote=F,row.names = F)
showNum=30
if(nrow(KEGG)<showNum){
	showNum=nrow(KEGG)
}

pdf(file="barplot.pdf",width = 10,height = 7)
barplot(kk, drop = TRUE, showCategory = showNum, color = colorSel)
dev.off()

pdf(file="bubble.pdf",width = 10,height = 7)
dotplot(kk, showCategory = showNum, orderBy = "GeneRatio",color = colorSel)
dev.off()

pdf(file="circos.pdf",width = 11,height = 7)
kkx=setReadable(kk, 'org.Hs.eg.db', 'ENTREZID')
cnetplot(kkx, foldChange=geneFC,showCategory = 5, circular = TRUE, colorEdge = TRUE,node_label="all")
dev.off()

14. protein-protein interaction network was performed by STRING online tool and refined via cytoscape.

15. COX model test was conducted by R with relevant R script.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("limma")
library(limma)
expFile="symbol.txt"                                         
cliFile="time.txt"                                           
geneFile="intersectGenes.txt"                                
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\21.preCox")     
rt=read.table(expFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0,]
gene=read.table(geneFile,header=T,sep="\t",check.names=F)
data=data[as.vector(gene[,1]),]
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]
colnames(data)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3",colnames(data))
data=t(data)
data=avereps(data)
cli=read.table(cliFile,sep="\t",check.names=F,header=T,row.names=1)     
sameSample=intersect(row.names(data),row.names(cli))
data=data[sameSample,]
cli=cli[sameSample,]
out=cbind(cli,data)
out=cbind(id=row.names(out),out)
write.table(out,file="expTime.txt",sep="\t",row.names=F,quote=F)
install.packages("survival")
library(survival)             
inputFile="expTime.txt"       
pFilter=0.05                  
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\22.cox")       
rt=read.table(inputFile,header=T,sep="\t",check.names=F,row.names=1)
rt$futime=rt$futime/365           
outTab=data.frame()
for(gene in colnames(rt[,3:ncol(rt)])){
	cox=coxph(Surv(futime, fustat) ~ rt[,gene], data = rt)
	coxSummary = summary(cox)
	coxP=coxSummary$coefficients[,"Pr(>|z|)"]
	group=ifelse(rt[,gene]<=median(rt[,gene]),"Low","High")
	diff=survdiff(Surv(futime, fustat) ~ group,data = rt)
	pValue=1-pchisq(diff$chisq,df=1)
	if((pValue<pFilter) & (coxP<pFilter)){
		outTab=rbind(outTab,
		         cbind(gene=gene,
		               KM=pValue,
		               HR=coxSummary$conf.int[,"exp(coef)"],
		               HR.95L=coxSummary$conf.int[,"lower .95"],
		               HR.95H=coxSummary$conf.int[,"upper .95"],
				       pvalue=coxP) )		
	}
}
write.table(outTab,file="cox.result.txt",sep="\t",row.names=F,quote=F)
rt <- read.table("cox.result.txt",header=T,sep="\t",row.names=1,check.names=F)
gene <- rownames(rt)
hr <- sprintf("%.3f",rt$"HR")
hrLow  <- sprintf("%.3f",rt$"HR.95L")
hrHigh <- sprintf("%.3f",rt$"HR.95H")
Hazard.ratio <- paste0(hr,"(",hrLow,"-",hrHigh,")")
pVal <- ifelse(rt$pvalue<0.001, "<0.001", sprintf("%.3f", rt$pvalue))
pdf(file="forest.pdf", width = 7,height = 6)
n <- nrow(rt)
nRow <- n+1
ylim <- c(1,nRow)
layout(matrix(c(1,2),nc=2),width=c(3,2))
xlim = c(0,3)
par(mar=c(4,2.5,2,1))
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")
text.cex=0.8
text(0,n:1,gene,adj=0,cex=text.cex)
text(1.5-0.5*0.2,n:1,pVal,adj=1,cex=text.cex);text(1.5-0.5*0.2,n+1,'pvalue',cex=text.cex,font=2,adj=1)
text(3,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1,)
par(mar=c(4,1,2,1),mgp=c(2,0.5,0))
xlim = c(0,max(as.numeric(hrLow),as.numeric(hrHigh)))
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazard ratio")
arrows(as.numeric(hrLow),n:1,as.numeric(hrHigh),n:1,angle=90,code=3,length=0.05,col="darkblue",lwd=2.5)
abline(v=1,col="black",lty=2,lwd=2)
boxcolor = ifelse(as.numeric(hr) > 1, 'red', 'green')
points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=1.3)
axis(1)
dev.off()

16. To screen the intersection gene between prognostic-related genes and differential express genes from PPI network with R script.
install.packages("VennDiagram")
library(VennDiagram)             
ppiFile="ppiCount.txt"           
coxFile="cox.result.txt"         
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\23.geneVenn")     
geneList=list()
ppi=read.table(ppiFile, header=T, sep="\t", check.names=F,row.names=1)
geneNum=30
if(nrow(ppi)<geneNum){
	geneNum=nrow(KEGG)
}
geneList[["PPI"]]=row.names(ppi)[1:geneNum]
cox=read.table(coxFile, header=T, sep="\t", check.names=F,row.names=1)
geneList[["COX"]]=row.names(cox)
venn.plot=venn.diagram(geneList,filename=NULL,main="Intersection",main.cex = 2,
             fill=c("darkblue", "darkgreen"),cat.cex=1.2)
pdf(file="geneVenn.pdf",width=7,height=7)
grid.draw(venn.plot)
dev.off()
coreGene=Reduce(intersect,geneList)
write.table(file="prognosticCoreGene.txt",coreGene,sep="\t",quote=F,row.names=F,col.names=F)

17. The difference express analysis of CD247.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")
install.packages("beeswarm")
library(limma)
library(beeswarm)
inputFile="symbol.txt"        
gene="CD247"                  
yMin=0                        
yMax=165                      
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\24.geneDiff")      
rt=read.table(inputFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0,]
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
tumor=data[,group==0]
tumorNum=ncol(tumor)
normalNum=ncol(data)-tumorNum
Type=c(rep(1,normalNum),rep(2,tumorNum))
rt=rbind(expression=data[gene,],Type=Type)
rt=as.matrix(t(rt))
wilcoxTest=wilcox.test(expression ~ Type, data=rt)
pvalue=wilcoxTest$p.value
if(pvalue<0.001){
     pvalue="<0.001"
}else{
     pvalue=paste0("=",sprintf("%.03f",pvalue))
}
pdf(file="diff.pdf",width=6,height=5)
par(mar = c(4,6,3,3))
labels=c("Normal","Tumor")
boxplot(expression ~ Type, data = rt,names=labels,xlab="",
     ylab = paste(gene," expression",sep=""),
     cex.main=1.5, cex.lab=1.3, cex.axis=1.2,ylim=c(yMin,yMax),outline = FALSE)
beeswarm(expression ~ Type, data = rt, col = c("blue","red"),lwd=0.1,
     pch = 16, add = TRUE, corral="wrap")
ySeg=yMax*0.94
segments(1,ySeg, 2,ySeg);segments(1,ySeg, 1,ySeg*0.96);segments(2,ySeg, 2,ySeg*0.96)
text(1.5,ySeg*1.05,labels=paste("p",pvalue,sep=""),cex=1.2)
dev.off()

18. The pair difference expression of CD247 in tumor and healthy samples processed by R script.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")
install.packages("ggpubr")
library("limma")
library("ggpubr")
inputFile="symbol.txt"        
gene="CD247"                  
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\25.genePairDiff")      

rt=read.table(inputFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0,]
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
normal=data[,group!=0]
tumor=data[,group==0]
normal=rbind(normal,gene=normal[gene,])
normal=as.matrix(t(normal[c("gene",gene),]))
rownames(normal)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3",rownames(normal))
normal=avereps(normal)
tumor=rbind(tumor,gene=tumor[gene,])
tumor=as.matrix(t(tumor[c("gene",gene),]))
rownames(tumor)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3",rownames(tumor))
tumor=avereps(tumor)
samSample=intersect(row.names(normal),row.names(tumor))
normal=normal[samSample,]
tumor=tumor[samSample,]
data=cbind(normal,tumor)
data=as.data.frame(data[,c(1,3)])
colnames(data)=c("Normal","Tumor")
pdf(file="pairDiff.pdf",width=5.5,height=5)
ggpaired(data, cond1 = "Normal", cond2 = "Tumor",fill = "condition", palette = "jco",
    xlab="",ylab = paste0(gene," expression"))+
    stat_compare_means(paired = TRUE, label = "p.format", label.x = 1.35)
dev.off()

19. The survival analysis of ten immune-related genes with R script.
install.packages("survival")
install.packages("survminer")
library(survival)
library(survminer)
inputFile="expTime.txt"      
gene="CD247"                 
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\26.geneSur")    
rt=read.table(inputFile,header=T,sep="\t",check.names=F)
rt$futime=rt$futime/365        
a=ifelse(rt[,gene]<=median(rt[,gene]),"Low","High")
diff=survdiff(Surv(futime, fustat) ~a,data = rt)
pValue=1-pchisq(diff$chisq,df=1)
fit=survfit(Surv(futime, fustat) ~ a, data = rt)
if(pValue<0.001){
	pValue="p<0.001"
}else{
	pValue=paste0("p=",sprintf("%.03f",pValue))
}
titleName=gene
surPlot=ggsurvplot(fit, 
				data=rt,
				conf.int=TRUE,
				pval=pValue,
				pval.size=6,
				risk.table=T,
				#ncensor.plot = TRUE,
				legend.labs=c("high","low"),
				legend.title=titleName,
				xlab="Time(years)",
				break.time.by = 1,
				risk.table.title="",
				palette=c("red", "blue"),
				risk.table.height=.25)          
pdf(file=paste0("sur.",gene,".pdf"), width = 6.5, height = 5.5,onefile = FALSE)
print(surPlot)
dev.off()

20. The clinical correlation analysis of CD247.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")
install.packages("ggpubr")
options(stringsAsFactors=F)
library(limma)
library(ggpubr)
inputFile="symbol.txt"          
cliFile="clinical.txt"          
gene="CD247"                    #
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\27.geneCliCor")     

rt=read.table(inputFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0,]

group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]
data=rbind(data,gene=data[gene,])
exp=as.matrix(t(data[c("gene",gene),]))
rownames(exp)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*","\\1\\-\\2\\-\\3",rownames(exp))
exp=avereps(exp)
cli=read.table(cliFile,sep="\t",header=T,check.names=F,row.names=1)
samSample=intersect(row.names(exp),row.names(cli))
exp=exp[samSample,]
cli=cli[samSample,]
rt=cbind(exp,cli)
for(clinical in colnames(rt[,3:ncol(rt)])){
	data=rt[c(gene,clinical)]
	colnames(data)=c("gene","clinical")
	data=data[(data[,"clinical"]!="unknow"),]
	#设置比较组
	group=levels(factor(data$clinical))
	data$clinical=factor(data$clinical, levels=group)
	comp=combn(group,2)
	my_comparisons=list()
    for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]}
	boxplot=ggboxplot(data, x="clinical", y="gene", color="clinical",
	          xlab=clinical,
	          ylab=paste(gene,"expression"),
	          legend.title=clinical,
	          add = "jitter")+ 
	stat_compare_means(comparisons = my_comparisons)
	pdf(file=paste0(clinical,".pdf"),width=5.5,height=5)
	print(boxplot)
	dev.off()
}

21. WGCNA with R script.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("GO.db", "preprocessCore", "impute","limma"))
install.packages(c("matrixStats", "Hmisc", "foreach", "doParallel", "fastcluster", "dynamicTreeCut", "survival")) 
install.packages("WGCNA")
library("WGCNA")            
library("limma")            
expFile="TCGA_rpkm.txt"     
normalCount=32              
tumorCount=436              
setwd("C:\\Users\\lexb4\\Desktop\\WGCNA\\11.TCGAwgcna")       
rt=read.table(expFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=log2(data+1)
data=data[apply(data,1,sd)>0.5,]
datExpr0=t(data)
gsg = goodSamplesGenes(datExpr0, verbose = 3)
if (!gsg$allOK)
{
  # Optionally, print the gene and sample names that were removed:
  if (sum(!gsg$goodGenes)>0)
    printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")))
  if (sum(!gsg$goodSamples)>0)
    printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")))
  # Remove the offending genes and samples from the data:
  datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}

sampleTree = hclust(dist(datExpr0), method = "average")
pdf(file = "1_sample_cluster.pdf", width = 12, height = 9)
par(cex = 0.6)
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)
abline(h = 20000, col = "red")
dev.off()
clust = cutreeStatic(sampleTree, cutHeight = 20000, minSize = 10)
table(clust)
keepSamples = (clust==1)
datExpr0 = datExpr0[keepSamples, ]
traitData=data.frame(Normal=c(rep(1,normalCount),rep(0,tumorCount)),
                     Tumor=c(rep(0,normalCount),rep(1,tumorCount)))
row.names(traitData)=colnames(data)
fpkmSamples = rownames(datExpr0)
traitSamples =rownames(traitData)
sameSample=intersect(fpkmSamples,traitSamples)
datExpr0=datExpr0[sameSample,]
datTraits=traitData[sameSample,]
sampleTree2 = hclust(dist(datExpr0), method = "average")
traitColors = numbers2colors(datTraits, signed = FALSE)
pdf(file="2_sample_heatmap.pdf",width=15,height=12)
plotDendroAndColors(sampleTree2, traitColors,
                    groupLabels = names(datTraits),
                    main = "Sample dendrogram and trait heatmap")
dev.off()
enableWGCNAThreads()   
powers = c(1:20)       
sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)
pdf(file="3_scale_independence.pdf",width=9,height=5)
par(mfrow = c(1,2))
cex1 = 0.9
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
abline(h=0.90,col="red") 
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
dev.off()
sft 
softPower =sft$powerEstimate #最佳power值
adjacency = adjacency(datExpr0, power = softPower)
softPower
TOM = TOMsimilarity(adjacency);
dissTOM = 1-TOM
geneTree = hclust(as.dist(dissTOM), method = "average");
pdf(file="4_gene_clustering.pdf",width=12,height=9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
     labels = FALSE, hang = 0.04)
dev.off()
minModuleSize = 50      
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                            deepSplit = 2, pamRespectsDendro = FALSE,
                            minClusterSize = minModuleSize);
table(dynamicMods)
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
pdf(file="5_Dynamic_Tree.pdf",width=8,height=6)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05,
                    main = "Gene dendrogram and module colors")
dev.off()
MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
MEs = MEList$eigengenes
MEDiss = 1-cor(MEs);
METree = hclust(as.dist(MEDiss), method = "average")
pdf(file="6_Clustering_module.pdf",width=7,height=6)
plot(METree, main = "Clustering of module eigengenes",
     xlab = "", sub = "")
MEDissThres = 0.25 
abline(h=MEDissThres, col = "red")
dev.off()
merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)
mergedColors = merge$colors
mergedMEs = merge$newMEs
pdf(file="7_merged_dynamic.pdf", width = 9, height = 6)
plotDendroAndColors(geneTree, mergedColors,"Dynamic Tree Cut",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05,
                    main = "Gene dendrogram and module colors(TCGA)")
dev.off()
moduleColors = mergedColors
table(moduleColors)
colorOrder = c("grey", standardColors(50))
moduleLabels = match(moduleColors, colorOrder)-1
MEs = mergedMEs
nGenes = ncol(datExpr0)
nSamples = nrow(datExpr0)
moduleTraitCor = cor(MEs, datTraits, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
pdf(file="8_Module_trait.pdf",width=6,height=6)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(5, 10, 3, 3))
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(datTraits),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships(TCGA)"))
dev.off()
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr0, MEs, use = "p"))
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
names(geneModuleMembership) = paste("MM", modNames, sep="")
names(MMPvalue) = paste("p.MM", modNames, sep="")
traitNames=names(datTraits)
geneTraitSignificance = as.data.frame(cor(datExpr0, datTraits, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) = paste("GS.", traitNames, sep="")
names(GSPvalue) = paste("p.GS.", traitNames, sep="")
for (trait in traitNames){
  traitColumn=match(trait,traitNames)  
  for (module in modNames){
    column = match(module, modNames)
    moduleGenes = moduleColors==module
    if (nrow(geneModuleMembership[moduleGenes,]) > 1){
      outPdf=paste("9_", trait, "_", module,".pdf",sep="")
      pdf(file=outPdf,width=7,height=7)
      par(mfrow = c(1,1))
      verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                         abs(geneTraitSignificance[moduleGenes, traitColumn]),
                         xlab = paste("Module Membership in", module, "module"),
                         ylab = paste("Gene significance for ",trait),
                         main = paste("Module membership vs. gene significance\n"),
                         cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
      abline(v=0.8,h=0.5,col="red")
      dev.off()
    }
  }
}
probes = colnames(datExpr0)
geneInfo0 = data.frame(probes= probes,
                       moduleColor = moduleColors)
for (Tra in 1:ncol(geneTraitSignificance))
{
  oldNames = names(geneInfo0)
  geneInfo0 = data.frame(geneInfo0, geneTraitSignificance[,Tra],
                         GSPvalue[, Tra])
  names(geneInfo0) = c(oldNames,names(geneTraitSignificance)[Tra],
                       names(GSPvalue)[Tra])
}

for (mod in 1:ncol(geneModuleMembership))
{
  oldNames = names(geneInfo0)
  geneInfo0 = data.frame(geneInfo0, geneModuleMembership[,mod],
                         MMPvalue[, mod])
  names(geneInfo0) = c(oldNames,names(geneModuleMembership)[mod],
                       names(MMPvalue)[mod])
}
geneOrder =order(geneInfo0$moduleColor)
geneInfo = geneInfo0[geneOrder, ]
write.table(geneInfo, file = "GS_MM.xls",sep="\t",row.names=F)
for (mod in 1:nrow(table(moduleColors)))
{  
  modules = names(table(moduleColors))[mod]
  probes = colnames(datExpr0)
  inModule = (moduleColors == modules)
  modGenes = probes[inModule]
  write.table(modGenes, file =paste0("TCGA_",modules,".txt"),sep="\t",row.names=F,col.names=F,quote=F)
}

22. Immune cells infiltration processed by R script.
CIBERSORT R script v1.03
#' Note: Signature matrix construction is not currently available; use java version for full functionality.
#' Author: Aaron M. Newman, Stanford University (amnewman@stanford.edu)
#' Requirements:
#'       R v3.0 or later. (dependencies below might not work properly with earlier versions)
#'       install.packages('e1071')
#'       install.pacakges('parallel')
#'       install.packages('preprocessCore')
#'       if preprocessCore is not available in the repositories you have selected, run the following:
#'           source("http://bioconductor.org/biocLite.R")
#'           biocLite("preprocessCore")
#' Windows users using the R GUI may need to Run as Administrator to install or update packages.
#' This script uses 3 parallel processes.  Since Windows does not support forking, this script will run
#' single-threaded in Windows.
#'
#' Usage:
#'       Navigate to directory containing R script
#'
#'   In R:
#'       source('CIBERSORT.R')
#'       results <- CIBERSORT('sig_matrix_file.txt','mixture_file.txt', perm, QN)
#'
#'       Options:
#'       i)  perm = No. permutations; set to >=100 to calculate p-values (default = 0)
#'       ii) QN = Quantile normalization of input mixture (default = TRUE)
#'
#' Input: signature matrix and mixture file, formatted as specified at http://cibersort.stanford.edu/tutorial.php
#' Output: matrix object containing all results and tabular data written to disk 'CIBERSORT-Results.txt'
#' License: http://cibersort.stanford.edu/CIBERSORT_License.txt
#' Core algorithm
#' @param X cell-specific gene expression
#' @param y mixed expression per sample
#' @export
CoreAlg <- function(X, y){
  #try different values of nu
  svn_itor <- 3
  res <- function(i){
    if(i==1){nus <- 0.25}
    if(i==2){nus <- 0.5}
    if(i==3){nus <- 0.75}
    model<-svm(X,y,type="nu-regression",kernel="linear",nu=nus,scale=F)
    model
  }
  if(Sys.info()['sysname'] == 'Windows') out <- mclapply(1:svn_itor, res, mc.cores=1) else
    out <- mclapply(1:svn_itor, res, mc.cores=svn_itor)
  nusvm <- rep(0,svn_itor)
  corrv <- rep(0,svn_itor)
  #do cibersort
  t <- 1
  while(t <= svn_itor) {
    weights = t(out[[t]]$coefs) %*% out[[t]]$SV
    weights[which(weights<0)]<-0
    w<-weights/sum(weights)
    u <- sweep(X,MARGIN=2,w,'*')
    k <- apply(u, 1, sum)
    nusvm[t] <- sqrt((mean((k - y)^2)))
    corrv[t] <- cor(k, y)
    t <- t + 1
  }
  pick best model
  rmses <- nusvm
  mn <- which.min(rmses)
  model <- out[[mn]]
  get and normalize coefficients
  q <- t(model$coefs) %*% model$SV
  q[which(q<0)]<-0
  w <- (q/sum(q))
  mix_rmse <- rmses[mn]
  mix_r <- corrv[mn]
  newList <- list("w" = w, "mix_rmse" = mix_rmse, "mix_r" = mix_r)
}
#' do permutations
#' @param perm Number of permutations
#' @param X cell-specific gene expression
#' @param y mixed expression per sample
#' @export
doPerm <- function(perm, X, Y){
  itor <- 1
  Ylist <- as.list(data.matrix(Y))
  dist <- matrix()
  while(itor <= perm){
    print(itor)
    random mixture
    yr <- as.numeric(Ylist[sample(length(Ylist),dim(X)[1])])
    standardize mixture
    yr <- (yr - mean(yr)) / sd(yr)
    run CIBERSORT core algorithm
    result <- CoreAlg(X, yr)
    mix_r <- result$mix_r
    store correlation
    if(itor == 1) {dist <- mix_r}
    else {dist <- rbind(dist, mix_r)}
    itor <- itor + 1
  }
  newList <- list("dist" = dist)
}
#' Main functions
#' @param sig_matrix file path to gene expression from isolated cells
#' @param mixture_file heterogenous mixed expression
#' @param perm Number of permutations
#' @param QN Perform quantile normalization or not (TRUE/FALSE)
#' @export
CIBERSORT <- function(sig_matrix, mixture_file, perm=0, QN=TRUE){
  library(e1071)
  library(parallel)
  library(preprocessCore)
  read in data
  X <- read.table(sig_matrix,header=T,sep="\t",row.names=1,check.names=F)
  Y <- read.table(mixture_file, header=T, sep="\t", row.names=1,check.names=F)
  X <- data.matrix(X)
  Y <- data.matrix(Y)
  order
  X <- X[order(rownames(X)),]
  Y <- Y[order(rownames(Y)),]
  P <- perm #number of permutations
  anti-log if max < 50 in mixture file
  if(max(Y) < 50) {Y <- 2^Y}
  quantile normalization of mixture file
  if(QN == TRUE){
    tmpc <- colnames(Y)
    tmpr <- rownames(Y)
    Y <- normalize.quantiles(Y)
    colnames(Y) <- tmpc
    rownames(Y) <- tmpr
  }
  intersect genes
  Xgns <- row.names(X)
  Ygns <- row.names(Y)
  YintX <- Ygns %in% Xgns
  Y <- Y[YintX,]
  XintY <- Xgns %in% row.names(Y)
  X <- X[XintY,]
  standardize sig matrix
  X <- (X - mean(X)) / sd(as.vector(X))
  empirical null distribution of correlation coefficients
  if(P > 0) {nulldist <- sort(doPerm(P, X, Y)$dist)}
  print(nulldist)
  header <- c('Mixture',colnames(X),"P-value","Correlation","RMSE")
  print(header)
  output <- matrix()
  itor <- 1
  mixtures <- dim(Y)[2]
  pval <- 9999
  iterate through mixtures
  while(itor <= mixtures){
    y <- Y[,itor]
    standardize mixture
    y <- (y - mean(y)) / sd(y)
    run SVR core algorithm
    result <- CoreAlg(X, y)
    get results
    w <- result$w
    mix_r <- result$mix_r
    mix_rmse <- result$mix_rmse
    calculate p-value
    if(P > 0) {pval <- 1 - (which.min(abs(nulldist - mix_r)) / length(nulldist))}
    #print output
    out <- c(colnames(Y)[itor],w,pval,mix_r,mix_rmse)
    if(itor == 1) {output <- out}
    else {output <- rbind(output, out)}
    itor <- itor + 1
  }
  save results
  write.table(rbind(header,output), file="CIBERSORT-Results.txt", sep="\t", row.names=F, col.names=F, quote=F)
  return matrix object containing all results
  obj <- rbind(header,output)
  obj <- obj[,-1]
  obj <- obj[-1,]
  obj <- matrix(as.numeric(unlist(obj)),nrow=nrow(obj))
  rownames(obj) <- colnames(Y)
  colnames(obj) <- c(colnames(X),"P-value","Correlation","RMSE")
  obj
}
install.packages('e1071')
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("preprocessCore")
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")
library("limma")             
expFile="symbol.txt"         
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\31.CIBERSORT)
rt=read.table(expFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]
data=data[rowMeans(data)>0,]
v <-voom(data, plot = F, save.plot = F)
out=v$E
out=rbind(ID=colnames(out),out)
write.table(out,file="uniq.symbol.txt",sep="\t",quote=F,col.names=F)        
source("TMEimmune31.CIBERSORT.R")
results=CIBERSORT("ref.txt", "uniq.symbol.txt", perm=100, QN=TRUE)

23. immune cells infiltration was shown by R script.
install.packages("corrplot")
library(corrplot)                   
input="CIBERSORT-Results.txt"       
outpdf="barplot.pdf"                
pFilter=0.05                        
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\32.immunePlot")    
immune=read.table("CIBERSORT-Results.txt",sep="\t",header=T,row.names=1,check.names=F)
immune=immune[immune[,"P-value"]<pFilter,]
immune=as.matrix(immune[,1:(ncol(immune)-3)])
data=t(immune)
col=rainbow(nrow(data),s=0.7,v=0.7)
pdf(outpdf,height=10,width=22)
par(las=1,mar=c(8,5,4,16),mgp=c(3,0.1,0),cex.axis=1.5)
a1 = barplot(data,col=col,yaxt="n",ylab="Relative Percent",xaxt="n",cex.lab=1.8)
a2=axis(2,tick=F,labels=F)
axis(2,a2,paste0(a2*100,"%"))
axis(1,a1,labels=F)
par(srt=60,xpd=T);text(a1,-0.02,colnames(data),adj=1,cex=0.6);par(srt=0)
ytick2 = cumsum(data[,ncol(data)])
ytick1 = c(0,ytick2[-length(ytick2)])
legend(par('usr')[2]*0.98,par('usr')[4],legend=rownames(data),col=col,pch=15,bty="n",cex=1.3)
dev.off()
pdf("corrplot.pdf",height=14,width=14)              
par(oma=c(0.5,1,1,1.2))
M=cor(immune)
corrplot(M, order = "AOE", type = "upper", tl.pos = "lt")
corrplot(M, add = TRUE, type = "lower", method = "number", order = "AOE",
         col = "black", diag = FALSE, tl.pos = "n", cl.pos = "n")
dev.off()

24. Violin plot to display immune cell fraction with R script.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")
install.packages("vioplot")
library(limma)
library(vioplot)
immuneFile="CIBERSORT-Results.txt"         
expFile="symbol.txt"                       
gene="CD247"                               
pFilter=0.05                               
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\33.vioplot")       
immune=read.table(immuneFile,sep="\t",header=T,row.names=1,check.names=F)
immune=immune[immune[,"P-value"]<pFilter,]
immune=as.matrix(immune[,1:(ncol(immune)-3)])
rt=read.table(expFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]
lowName=colnames(data)[data[gene,]<=median(data[gene,])]
highName=colnames(data)[data[gene,]>median(data[gene,])]
lowImm=intersect(row.names(immune),lowName)
highImm=intersect(row.names(immune),highName)
rt=rbind(immune[lowImm,],immune[highImm,])
lowNum=length(lowImm)
highNum=length(highImm)
outTab=data.frame()
pdf("vioplot.pdf",height=8,width=13)
par(las=1,mar=c(10,6,3,3))
x=c(1:ncol(rt))
y=c(1:ncol(rt))
plot(x,y,
     xlim=c(0,63),ylim=c(min(rt),max(rt)+0.02),
     main="",xlab="", ylab="Fraction",
     pch=21,
     col="white",
     xaxt="n")
for(i in 1:ncol(rt)){
	  if(sd(rt[1:lowNum,i])==0){
	    rt[1,i]=0.001
	  }
	  if(sd(rt[(lowNum+1):(lowNum+highNum),i])==0){
	    rt[(lowNum+1),i]=0.001
	  }
	  lowData=rt[1:lowNum,i]
	  highData=rt[(lowNum+1):(lowNum+highNum),i]
	  vioplot(lowData,at=3*(i-1),lty=1,add = T,col = 'green')
	  vioplot(highData,at=3*(i-1)+1,lty=1,add = T,col = 'red')
	  wilcoxTest=wilcox.test(lowData,highData)
	  p=wilcoxTest$p.value
	  if(p<pFilter){
	      cellPvalue=cbind(Cell=colnames(rt)[i],pvalue=p)
		  outTab=rbind(outTab,cellPvalue)
	  }
	  mx=max(c(lowData,highData))
	  lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx))
	  text(x=3*(i-1)+0.5, y=mx+0.02, labels=ifelse(p<0.001, paste0("p<0.001"), paste0("p=",sprintf("%.03f",p))), cex = 0.8)
}
text(seq(1,64,3),-0.05,xpd = NA,labels=colnames(rt),cex = 1,srt = 45,pos=2)
dev.off()
write.table(outTab,file="diff.result.txt",sep="\t",row.names=F,quote=F)

25. Clinical correlation between CD247 and immune cells with R script.
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma")
install.packages("ggplot2")
install.packages("ggpubr")
install.packages("ggExtra")
library(limma)
library(ggplot2)
library(ggpubr)
library(ggExtra)
immuneFile="CIBERSORT-Results.txt"         
expFile="symbol.txt"                       
gene="CD247"                               
pFilter=0.05                               
setwd("C:\\Users\\lexb4\\Desktop\\TMEimmune\\34.immuneCor")     
immune=read.table(immuneFile,sep="\t",header=T,row.names=1,check.names=F)
immune=immune[immune[,"P-value"]<pFilter,]
immune=as.matrix(immune[,1:(ncol(immune)-3)])
rt=read.table(expFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]
data=rbind(data,gene=data[gene,])
exp=t(data[c("gene",gene),])
sameSample=intersect(row.names(immune),row.names(exp))
immune1=immune[sameSample,]
exp1=exp[sameSample,]
outTab=data.frame()
x=as.numeric(exp1[,1])
for(j in colnames(immune1)[1:22]){
	y=as.numeric(immune1[,j])
	if(sd(y)>0.001){
		df1=as.data.frame(cbind(x,y))
		corT=cor.test(x,y,method="spearman")
		cor=corT$estimate
		pValue=corT$p.value
		p1=ggplot(df1, aes(x, y)) + 
			ylab(j)+xlab(gene)+
			geom_point()+ geom_smooth(method="lm",formula=y~x) + theme_bw()+
			stat_cor(method = 'spearman', aes(x =x, y =y))
		if(pValue<pFilter){
			pdf(file=paste0(j,".pdf"),width=5,height=5)
			print(p1)
			dev.off()
			outTab=rbind(outTab,cbind(Cell=j,pValue))
		}
	}
}
write.table(outTab,file="immuneCor.result.txt",sep="\t",row.names=F,quote=F)