-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path45hr_cuffdiff.R
84 lines (74 loc) · 3.14 KB
/
45hr_cuffdiff.R
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
setwd(".//2mapping")
library(cummeRbund)
cuff_data <- readCufflinks('45hr_diff_out/')
cuff_data
pdf(file="45hr_Density.pdf")
csDensity(genes(cuff_data))
dev.off()
pdf(file="45hr_54_58_Scatter.pdf")
csScatter(genes(cuff_data), 'S54', 'S58')
dev.off()
pdf(file="45hr_61_58_Scatter.pdf")
csScatter(genes(cuff_data), 'S61', 'S58')
dev.off()
pdf(file="45hr_63_58_Scatter.pdf")
csScatter(genes(cuff_data), 'S63', 'S58')
dev.off()
pdf(file="45hr_65_58_Scatter.pdf")
csScatter(genes(cuff_data), 'S65', 'S58')
dev.off()
pdf(file="45hr_63_58_Volcano.pdf")
csVolcano(genes(cuff_data), 'S63', 'S58',showSignificant=TRUE, alpha=0.05, features=TRUE)
dev.off()
pdf(file="45hr_54_58_Volcano.pdf")
csVolcano(genes(cuff_data), 'S54', 'S58',showSignificant=TRUE, alpha=0.05, features=TRUE)
dev.off()
pdf(file="45hr_61_58_Volcano.pdf")
csVolcano(genes(cuff_data), 'S61', 'S58',showSignificant=TRUE, alpha=0.05, features=TRUE)
dev.off()
pdf(file="45hr_65_58_Volcano.pdf")
csVolcano(genes(cuff_data), 'S65', 'S58',showSignificant=TRUE, alpha=0.05, features=TRUE)
dev.off()
## Get DE information from specific gene
#mygene<- getGene(cuff_data, 'SOCS1')
#expressionBarplot(mygene)
#expressionBarplot(isoforms(mygene))
#gene_diff_data <- diffData(genes(cuff_data))
#sig_gene_data <- subset(gene_diff_data, (significant == 'yes'))
#nrow(sig_gene_data)
#isoform_diff_data <- diffData(isoforms(cuff_data), 'S54', 'S58')
#sig_isoform_data <- subset(isoform_diff_data, (significant == 'yes'))
#nrow(sig_isoform_data)
#tss_diff_data <- diffData(TSS(cuff_data), 'S54', 'S58')
#sig_tss_data <- subset(tss_diff_data, (significant =='yes'))
#nrow(sig_tss_data)
#cds_diff_data <- diffData(CDS(cuff_data), 'S54', 'S58')
#sig_cds_data <- subset(cds_diff_data, (significant =='yes'))
#nrow(sig_cds_data)
#promoter_diff_data <- distValues(promoters(cuff_data))
#sig_promoter_data <- subset(promoter_diff_data, (significant == 'yes'))
#nrow(sig_promoter_data)
#splicing_diff_data <- distValues(splicing(cuff_data))
#sig_splicing_data <- subset(splicing_diff_data, (significant == 'yes'))
#nrow(sig_splicing_data)
#relCDS_diff_data <- distValues(relCDS(cuff_data))
#sig_relCDS_data <- subset(relCDS_diff_data, (significant == 'yes'))
#nrow(sig_relCDS_data)
#Retrieve significant gene IDs (XLOC) with a pre-specified alpha
diffGeneIDs <- getSig(cuff_data, level="genes", alpha=0.05)
#Use returned identifiers to create a CuffGeneSet object with all relevant info for given genes
diffGenes <- getGenes(cuff_data, diffGeneIDs)
#gene_short_name values (and corresponding XLOC_* values can be retrieved from the CuffGeneSet by this)
names <- featureNames(diffGenes)
row.names(names) = names$tracking_id
diffGeneNames<-as.matrix(names)
diffGeneNames<- diffGeneNames[,-1]
#get the data for the significant genes
diffGenesData <- diffData(diffGenes)
row.names(diffGenesData)=make.names(diffGenesData[,1],TRUE)
diffGenesData <- diffGenesData[,-1]
diffGenesOutput<-merge(diffGeneNames,diffGenesData, by="row.names")
diffGenesOutput<-diffGenesOutput[,-1]
colnames(diffGenesOutput)[1]<- "Gene"
sig_gene_data <- subset(diffGenesOutput, (significant =='yes'))
write.table(sig_gene_data, '45hr_diff_genes.txt', sep='\t', row.names = F, col.names = T, quote = F)