-
Notifications
You must be signed in to change notification settings - Fork 1
/
RNA-Seq.R
66 lines (41 loc) · 2.09 KB
/
RNA-Seq.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
#1. Read count genes data
unSortedCounts <- read.table(file = "C:\\Users/alexb/Desktop/Bioinformática/Algoritmos para el análisis de secuencias/TP3_RNASeq_Counts.txt", sep ="\t", stringsAsFactors = TRUE, header = TRUE, row.names = "Geneid")
countsFile <- unSortedCounts[ , order(names(unSortedCounts))]
#2. Read and order Sample Sheet
unSortedSheet <- read.table(file = "C:\\Users/alexb/Desktop/Bioinformática/Algoritmos para el análisis de secuencias/TP3_RNASeq_SampleSheet.txt", sep ="\t", stringsAsFactors = FALSE, header = TRUE)
targetsFile <- unSortedSheet[order(unSortedSheet$Sample),]
rownames(targetsFile) <- targetsFile$Sample
colnames(countsFile) <- targetsFile$Sample
#3. Load "DESeq2" and create DESeqDataSet
library("DESeq2")
ddsFullCountTable <- DESeqDataSetFromMatrix(countData = countsFile, colData = targetsFile, design = ~ Strain + Treatment)
#4. DESeq2 analysis.
dds <- DESeq(ddsFullCountTable)
resultsNames(dds)
#5. Create Control & Treatment objects.
factor_condition <- as.factor(dds$condition)
Control <- factor_condition == "Control"
Treatment <- factor_condition == "Treatment"
#6. Establish control group.
ddsFullCountTable$Treatment <- relevel(ddsFullCountTable$Treatment, "Control")
#7. dds results.
res <- results(dds)
head(res)
summary(res)
#8. Dispersion plot.
plotDispEsts(dds)
#9. Differentially expressed genes number.
is.na(res)
sum(is.na(res))
na.omit(res)
genes_differentially_expr <- sum(res$padj < 0.05, na.rm = TRUE)
genes_differentially_expr
#10. Order differentially expressed genes by p-value adjusted (padj).
res_orden <- res[order(na.omit(res$padj)), ]
write.table(res_orden, file = "analysis_differentially_genes_results.tsv", sep = "\t", quote = FALSE, row.names = TRUE)
differentially_expression_data <- read.table("analysis_differentially_genes_results.tsv", header = TRUE, sep = "\t")
differentially_expression_data
#11. Generate MAplots with different rates.
plotMA(res, ylim=c(-2,2), alpha = 0.05)
plotMA(res, ylim=c(-2,2), alpha = 0.001)
plotMA(res, ylim=c(-2,2), alpha = 0.000001)