Question: Not able to get DEG from FireBrowse data using R
0
gravatar for lawarde.ankita1
11 months ago by
lawarde.ankita130 wrote:

Hello everyone,

I'm working on LUAD (Lung adenocarcinoma) data from FireBrowse. I want to find Differentially expressed genes according to the stage of the tumor. i have the clinical data with me and from that data i was able to take the samples which belongs to stage I(sub types a and b), Stage II (subtypes a and b), Stage III (subtypes a and b) and Stage IV. LUAD.rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.data this is the file of RNAseq data with around 56000 genes and 570 samples.

Below i have pasted the R code i used to find the DEG number, reference link i used (Survival analysis of TCGA patients integrating gene expression (RNASeq) data)

But I'm not able to understand what is going wrong in this analysis as i'm not able to get the DEG in tumor samples of stage I. II , II ,IV.

1) If my analysis is wrong please let me know what wrong steps i'm using. 2) what particular Package i'm suppose to use? 3)Any sample code of R would really help me.

Any suggestions would be very helpful to clear my doubts.

#read RNA file

rna <- read.table("RNA/LUAD.rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.data.txt",nrows = 20533, header = T, row.names = 1, sep = "\t")
# remove the first row as we don't need it,

rna <- rna[-1,]
dim(rna)

#remove genes whose expression is ==0 there are 50% of such samples,
rem <- function(x){
  x <- as.matrix(x)
  x<- t(apply(x,1,as.numeric))
  r <- as.numeric(apply(x,1,function(i)sum(i == 0)))
  remove <- which(r > dim(x)[2]*0.5) 
  return(remove)
}
remove <- rem(rna)
rna <- rna[-remove,]
dim(rna)
colnames(rna)

colnames(rna) <- gsub('\\.','-',substr(colnames(rna),1,28))
#read files of tumor stage
library(xlsx)

Tumor_stage <- read.xlsx(file = "Clinical/Tumor_stage_data.xlsx", sheetIndex = 1)
colnames(Tumor_stage)
Tumor_stage$Stage_I_Barcode <- toupper(Tumor_stage$Stage_I_Barcode)


Tumor_stage$Stage_II_Barcode <- toupper(Tumor_stage$Stage_II_Barcode)
Tumor_stage$Stage_III_Barcode <- toupper(Tumor_stage$Stage_III_Barcode)
Tumor_stage$Stage_IV_Barcode <- toupper(Tumor_stage$Stage_IV_Barcode)

#select only tumor samples

Tumor_Samples <- subset(rna, select = substr(colnames(rna),14,14) == '0')
dim(Tumor_Samples)

#subset the stages of tumor samples from Tumor_samples
Stage_I <- subset(Tumor_Samples, select = substr(colnames(Tumor_Samples),1,12) %in% Tumor_stage$Stage_I_Barcode)
dim(Stage_I)
stage_I_colnames <- as.data.frame(colnames(Stage_I))
stage_I_colnames

write.xlsx(stage_I_colnames, file = "samples_stage.xlsx")

Stage_II <- subset(Tumor_Samples, select = substr(colnames(Tumor_Samples),1,12) %in% Tumor_stage$Stage_II_Barcode)
dim(Stage_II)
stage_II_colnames <- as.data.frame(colnames(Stage_II))
stage_II_colnames

write.xlsx(stage_II_colnames, file = "samples_stageII.xlsx")

Stage_III <- subset(Tumor_Samples, select = substr(colnames(Tumor_Samples),1,12) %in% Tumor_stage$Stage_III_Barcode)
dim(Stage_III)
stage_III_colnames <- as.data.frame(colnames(Stage_III))
#stage_II_colnames

write.xlsx(stage_III_colnames, file = "samples_stageIII.xlsx")

Stage_IV <- subset(Tumor_Samples, select = substr(colnames(Tumor_Samples),1,12) %in% Tumor_stage$Stage_IV_Barcode)
dim(Stage_IV)
stage_IV_colnames <- as.data.frame(colnames(Stage_IV))
#stage_II_colnames

write.xlsx(stage_IV_colnames, file = "samples_stageIV.xlsx")

Tumor_stage_data <- cbind(Stage_I, Stage_II, Stage_III, Stage_IV)
dim(Tumor_stage_data)

#groups of samples
groups <- c(rep("Stage_I",277), rep("Stage_II",122), rep("Stage_III",84),rep("Stage_IV",26))
group <- factor(groups)

group

#design matrix
design <- model.matrix(~ 0 + group)
design
colnames(design)<- c("Stage_I", "Stage_II", "Stage_III","Stage_IV" )
tail(design)
dim(design)
dim(Tumor_stage_data)

library(limma)
x <- t(apply(Tumor_stage_data,1,as.numeric))

#Testing for differential expression
# Fit the linear model

fit <- lmFit(x, design)
names(fit)
group
cont.matrix <- makeContrasts(I_II = Stage_I - Stage_II,
                             I_III = Stage_I - Stage_III,
                             I_IV = Stage_I - Stage_IV,
                             II_III = Stage_II- Stage_III,
                             II_IV = Stage_II - Stage_IV,
                             III_IV = Stage_III - Stage_IV, levels=design)
cont.matrix

fit2 <- eBayes(contrasts.fit(fit, cont.matrix))
dim(fit2)
colnames(fit2)
summa.fit <- decideTests(fit2)
summary(summa.fit)

tt_I_II <- topTable(fit2, coef ="I_II", adjust.method = "BH",p.value = 0.05,lfc = 2, number = nrow(x))
dim(tt_I_II)

tt_I_III <- topTreat(fit2, coef ="I_III",adjust.method = "BH",p.value = 0.05,lfc = 2, number = nrow(x))
dim(tt_I_III)

tt_I_IV <- topTable(fit2, coef ="I_IV",adjust.method = "BH",p.value = 0.05,lfc = 2, number = nrow(x))
dim(tt_I_IV)

tt_II_III <- topTable(fit2, coef ="II_III",adjust.method = "BH",p.value = 0.05,lfc = 2, number = nrow(x))
dim(tt_II_III)

tt_II_IV <- topTable(fit2, coef ="II_IV",adjust.method = "BH",p.value = 0.05,lfc = 2, number = nrow(x))
dim(tt_II_IV)

tt_III_IV <- topTable(fit2, coef ="III_IV",adjust.method = "BH",p.value = 0.05,lfc = 2, number = nrow(x))
dim(tt_III_IV)

##heatmap
colnames(x) <- c(rep("Stage_I",277), rep("Stage_II",122), rep("Stage_III",84),rep("Stage_IV",26)) 
head(x)

idx_I_III <- rownames(tt_I_III)
DEG_I_III <- x[idx_I_III,]
dim(DEG_I_III)
head(DEG_I_III)

hr <- hclust(as.dist(1-cor(t(DEG_I_III), method="pearson")),method="complete")
plot(as.dendrogram(hr))

hc <-  hclust(as.dist(1-cor(DEG_I_III, method="spearman")), method="complete")
plot(as.dendrogram(hc))
colnames(DEG_I_III)

similar for remaining Contrast groups to plot the heatmap as above.

When i follow above code i get DEG in four contrasting groups but when i plot the heatmap there is no differentiating pattern (i know the number of samples is more than 100, but if there is any valid reason behind it please let me know.)

there is one more way i tried to find DEG by only taking Stage I and IV, Stage III and IV, Stage II and IV but here, in this case the number DEG changed and in case of Stage III and IV the DEG number was zero (and for this same group in first code the DEG number was 160)

The heatmap images i have attached with this question

DEG between I and III stage

 stage I and IV

stage II and IV

Stage II and IV

deg firebrowse rna-seq R • 481 views
ADD COMMENTlink written 11 months ago by lawarde.ankita130
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1403 users visited in the last hour