-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathDEA+EA_Pipeline final.Rcode
More file actions
205 lines (168 loc) · 10.4 KB
/
DEA+EA_Pipeline final.Rcode
File metadata and controls
205 lines (168 loc) · 10.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
!!! = Project ID
### = GTEX tissue
+++ = GTEX tissue (CAPITALIZED)
@@@ = TCGA tissue
?????? = Disease type
## Initiate required packages ##
library(TCGAbiolinks)
library(SummarizedExperiment)
library(recount)
library(TCGAutils)
library(limma)
library(biomaRt)
library(maftools)
library(dplyr)
library(EnhancedVolcano)
library(pathfindR)
## Gene conversion function ##
convert.ENSG.Symbol <- function(genes){
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl", host = "http://uswest.ensembl.org"))
G_list <- getBM(filters = "ensembl_gene_id", attributes = c("ensembl_gene_id", "hgnc_symbol"), values = genes, mart = mart)
return(G_list)
}
## Create dataframe samples.Count once, DO NOT reiterate ##
samples.Count <- setNames(data.frame(matrix(ncol = 2, nrow = 0)), c("tumor", "normal"))
## Begin downloading respective TCGA and GTEx RNAseq (HTSeq - Counts) data ##
!!!.recount.gtex <- TCGAquery_recount2(project = "GTEX", tissue = "###")
!!!.recount.tcga <- TCGAquery_recount2(project = "TCGA", tissue = "@@@")
SE.!!!.recount.gtex <- !!!.recount.gtex$GTEX_###
SE.!!!.recount.tcga <- !!!.recount.tcga$TCGA_@@@
query.!!! <- GDCquery(project = "TCGA-!!!",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts")
samplesDown.!!! <- getResults(query.!!!, cols = c("cases"))
dataSmTP.!!! <- TCGAquery_SampleTypes(barcode = samplesDown.!!!,
typesample = "TP")
dataSmNT.!!! <- TCGAquery_SampleTypes(barcode = samplesDown.!!!,
typesample = "NT")
!!!.eset.gtex <- assays(scale_counts(!!!.recount.gtex$GTEX_###, round = TRUE))$counts
!!!.eset.tcga <- assays(scale_counts(!!!.recount.tcga$TCGA_@@@, round = TRUE))$counts
rse_scaled <- scale_counts(!!!.recount.gtex$GTEX_###, round = TRUE)
summary(colSums(assays(rse_scaled)$counts)) / 1e6
colnames(!!!.eset.tcga) <- colData(!!!.recount.tcga$TCGA_@@@)$gdc_cases.samples.portions.analytes.aliquots.submitter_id
rownames(!!!.eset.gtex) <- gsub("\\..*", "", rownames(!!!.eset.gtex))
rownames(!!!.eset.tcga) <- gsub("\\..*", "", rownames(!!!.eset.tcga))
!!!.eset.tcga.cancer <- !!!.eset.tcga[,which(colData(!!!.recount.tcga$TCGA_@@@)$gdc_cases.samples.sample_type=="Primary Tumor")]
!!!.eset.tcga.cancer <- !!!.eset.tcga.cancer[,c(dataSmTP.!!!)]
!!!.eset.tcga.normal <- !!!.eset.tcga[,which(colData(!!!.recount.tcga$TCGA_@@@)$gdc_cases.samples.sample_type=="Solid Tissue Normal")]
!!!.eset.tcga.normal <- !!!.eset.tcga.normal[,c(dataSmNT.!!!)]
dataPrep.!!! <- merge(as.data.frame(!!!.eset.gtex), as.data.frame(!!!.eset.tcga.cancer), by = 0, all = TRUE)
rownames(dataPrep.!!!) <- dataPrep.!!!$Row.names
dataPrep.!!!$Row.names <- NULL
dataNorm.!!! <- TCGAanalyze_Normalization(tabDF = dataPrep.!!!,
geneInfo = geneInfoHT,
method = "gcContent")
dataFilt.!!! <- TCGAanalyze_Filtering(tabDF = dataNorm.!!!,
method = "quantile",
qnt.cut = 0.25)
DEG.!!! <- TCGAanalyze_DEA(mat1 = dataFilt.!!![,colnames(!!!.eset.gtex)],
mat2 = dataFilt.!!![,colnames(!!!.eset.tcga.cancer)],
metadata = FALSE,
pipeline = "limma",
voom = TRUE,
Cond1type = "Normal",
Cond2type = "Tumor",
method = "glmLRT")
EA.DEG.!!! <- TCGAanalyze_DEA(mat1 = dataFilt.!!![,colnames(!!!.eset.gtex)],
mat2 = dataFilt.!!![,colnames(!!!.eset.tcga.cancer)],
metadata = FALSE,
pipeline = "limma",
voom = TRUE,
fdr.cut = 10e-16,
logFC.cut = 2,
Cond1type = "Normal",
Cond2type = "Tumor",
method = "glmLRT")
assign("last.warning", NULL, envir = baseenv())
!!!.conversion.table <- convert.ENSG.Symbol(rownames(DEG.!!!))
!!!.conversion.inter.DEG <- intersect(!!!.conversion.table[-which(!!!.conversion.table$hgnc_symbol==""),]$ensembl_gene_id, rownames(DEG.!!!))
!!!.conversion.table2 <- !!!.conversion.table[which(!!!.conversion.table$ensembl_gene_id %in% !!!.conversion.inter.DEG),]
rownames(!!!.conversion.table2) <- !!!.conversion.table2$ensembl_gene_id
## Error check and correction ##
dupl <- names(last.warning)
if (!is.null(dupl) == TRUE) {
dupl.sub <- noquote(sub(".*: ", "", dupl))
dupl.sub2 <- scan(text = dupl.sub, sep = ",", what = "")
dupl.sub3 <- gsub("[^[:alnum:] ]|\\s", "", dupl.sub2)
DEG.!!! <- DEG.!!![-which(rownames(DEG.!!!) %in% c(dupl.sub3)), ]
!!!.conversion.table <- convert.ENSG.Symbol(rownames(DEG.!!!))
!!!.conversion.inter.DEG <- intersect(!!!.conversion.table[-which(!!!.conversion.table$hgnc_symbol==""),]$ensembl_gene_id, rownames(DEG.!!!))
!!!.conversion.table2 <- !!!.conversion.table[which(!!!.conversion.table$ensembl_gene_id %in% !!!.conversion.inter.DEG),]
rownames(!!!.conversion.table2) <- !!!.conversion.table2$ensembl_gene_id
!!!.conversion.table[-which(!!!.conversion.table$hgnc_symbol==""),]
assign("last.warning", NULL, envir = baseenv())
} else {
!!!.conversion.table[-which(!!!.conversion.table$hgnc_symbol==""),]
}
!!!.EA.conversion.table <- convert.ENSG.Symbol(rownames(EA.DEG.!!!))
!!!.EA.conversion.inter.DEG <- intersect(!!!.EA.conversion.table[-which(!!!.EA.conversion.table$hgnc_symbol==""),]$ensembl_gene_id, rownames(EA.DEG.!!!))
!!!.EA.conversion.table2 <- !!!.EA.conversion.table[which(!!!.EA.conversion.table$ensembl_gene_id %in% !!!.EA.conversion.inter.DEG),]
rownames(!!!.EA.conversion.table2) <- !!!.EA.conversion.table2$ensembl_gene_id
## Error check and correction ##
dupl <- names(last.warning)
if (!is.null(dupl) == TRUE) {
dupl.sub <- noquote(sub(".*: ", "", dupl))
dupl.sub2 <- scan(text = dupl.sub, sep = ",", what = "")
dupl.sub3 <- gsub("[^[:alnum:] ]|\\s", "", dupl.sub2)
EA.DEG.!!! <- EA.DEG.!!![-which(rownames(EA.DEG.!!!) %in% c(dupl.sub3)), ]
!!!.EA.conversion.table <- convert.ENSG.Symbol(rownames(EA.DEG.!!!))
!!!.EA.conversion.inter.DEG <- intersect(!!!.EA.conversion.table[-which(!!!.EA.conversion.table$hgnc_symbol==""),]$ensembl_gene_id, rownames(EA.DEG.!!!))
!!!.EA.conversion.table2 <- !!!.EA.conversion.table[which(!!!.EA.conversion.table$ensembl_gene_id %in% !!!.EA.conversion.inter.DEG),]
rownames(!!!.EA.conversion.table2) <- !!!.EA.conversion.table2$ensembl_gene_id
!!!.EA.conversion.table[-which(!!!.EA.conversion.table$hgnc_symbol==""),]
assign("last.warning", NULL, envir = baseenv())
} else {
!!!.EA.conversion.table[-which(!!!.EA.conversion.table$hgnc_symbol==""),]
}
DEG.!!!.hgnc <- DEG.!!![!!!.conversion.inter.DEG,]
DEGs.!!!.hgnc <- merge(DEG.!!!.hgnc, !!!.conversion.table2, by = 0)
rownames(DEGs.!!!.hgnc) <- DEGs.!!!.hgnc$Row.names
DEGs.!!!.hgnc$Row.names <- NULL
rownames(DEGs.!!!.hgnc) <- DEGs.!!!.hgnc$hgnc_symbol
EA.DEG.!!!.hgnc <- EA.DEG.!!![!!!.EA.conversion.inter.DEG,]
EA.DEGs.!!!.hgnc <- merge(EA.DEG.!!!.hgnc, !!!.EA.conversion.table2, by = 0)
rownames(EA.DEGs.!!!.hgnc) <- EA.DEGs.!!!.hgnc$Row.names
EA.DEGs.!!!.hgnc$Row.names <- NULL
rownames(EA.DEGs.!!!.hgnc) <- EA.DEGs.!!!.hgnc$hgnc_symbol
EA.DEGs.!!!.hgnc <- subset(EA.DEGs.!!!.hgnc, select = -c(AveExpr, t, P.Value, B, ensembl_gene_id))
names(EA.DEGs.!!!.hgnc)[3] <- "Gene.symbol"
EA.DEGs.!!!.hgnc <- EA.DEGs.!!!.hgnc %>% select(Gene.symbol, everything())
PID.DEGs.!!!.hgnc <- subset(DEGs.!!!.hgnc, hgnc_symbol %in% Panel_PIDGenes$X1)
PID.EA.DEGs.!!!.hgnc <- subset(EA.DEGs.!!!.hgnc, Gene.symbol %in% Panel_PIDGenes$X1)
write.csv(PID.DEGs.!!!.hgnc, "PID.DEGs.!!!.csv")
write.csv(DEGs.!!!.hgnc, "DEGs.!!!.csv")
write.csv(PID.EA.DEGs.!!!.hgnc, "PID.EA.DEGs.!!!.csv")
write.csv(EA.DEGs.!!!.hgnc, "EA.DEGs.!!!.csv")
EnhancedVolcano(PID.DEGs.!!!.hgnc,
lab = rownames(PID.DEGs.!!!.hgnc),
x = 'logFC',
y = 'adj.P.Val',
ylab = bquote(~-Log[10]~ "(FDR corrected-P values)"),
title = "?????? (TCGA-!!!) vs. Normal (GTEX-+++)",
subtitle = "Differential expression using limma-voom; showing PID-related genes",
caption = "FC cutoff, 2; p-value cutoff, 10e-16",
pCutoff = 10e-16,
FCcutoff = 2,
pointSize = 6.0,
labSize = 6.0,
legendPosition = 'none')
dev.copy(tiff, "DEA_TCGA-!!!vsGTEX-+++.tiff", width=1000, height=1000)
dev.off()
## Run gene set (GO) enrichment analysis using pathfindR package ##
EA.DEGs.!!!.hgnc.pathfindR <- run_pathfindR(EA.DEGs.!!!.hgnc,
gene_sets = "KEGG",
output_dir = "R:\\Research\\Immunodeficiency\\Data\\EA_Analysis\\EA_TCGA-!!!vsGTEX-+++")
dev.copy(tiff, "EA_TCGA-!!!vsGTEX-+++.tiff", width=1000, height=1000)
dev.off()
PID.EA.!!!.hgnc.pathfindR <- run_pathfindR(PID.EA.DEGs.!!!.hgnc,
gene_sets = "KEGG",
output_dir = "R:\\Research\\Immunodeficiency\\Data\\EA_Analysis\\PID.EA_TCGA-!!!vsGTEX-+++")
dev.copy(tiff, "PID.EA_TCGA-!!!vsGTEX-+++.tiff", width=1000, height=1000)
dev.off()
!!!.samples.Count <- data.frame(ncol(!!!.eset.tcga.cancer), ncol(!!!.eset.gtex))
names(!!!.samples.Count) <- c("tumor", "normal")
samples.Count <- rbind(samples.Count, "TCGA-!!!" = !!!.samples.Count)
write.csv(samples.Count, "samples.Count.csv")
## Clear data except for Panel_PIDGenes and function ##
rm(list = setdiff(ls(), c("Panel_PIDGenes", "samples.Count", "convert.ENSG.Symbol")))