Online version script: https://rpubs.com/lincw/praktikumkurse

1 Analyze RNA-Seq data with R and Cytoscape

The aim of this part is to analyze and a host-pathogen protein interaction map derived from a Y2H experiment using the software “Cytoscape”. The network is extended with additional interactions and protein properties from online resources. Additionally, expression fold change values derived from a RNA-Seq experiment are integrated to examine the expression dynamic after pathogen treatment.

1.1 Softwares in this course

R

R is a free software environment for statistical computing and graphics. R

RStudio

RStudio provides free and open source tools for R and enterprise-ready professional software for data science teams to develop and share their work at scale. RStudio

2 GUI Basics

2.1 Software Interface

R

Figure 1. R environment.

RStudio

Figure 2. RStudio environment.

2.2 The environment of RStudio

The editor, where you write your program
The console, where you can enter R commands and where the output of your program is shown.
In the environment / history area are listed in the environment tab all variables, which are currently in use. In the history tab are shown the last commands, which have been executed in the console.
In the fourth area, five tabs are contained: Files, Plots, Packages, Help and Viewer.
- Files show your file system on your computer.
- Plots show the last produces graphic.
- Packages show all installed R packages and allows to install new ones from the CRAN repository.
- In the Help tab, information about R functions can be found.
- Viewer also shows java script graphics.
The Run button executes the currently selected lines in the editor.

A comprehensive list of shortcuts can be found here: shortcuts

3 Working environment setup

Packages are the fundamental units of reproducible R code. They include reusable R functions, the documentation that describes how to use them, and sample data . In the following analysis of RNA-Seq data several R-packages are needed.

BioManager: A convenient tool to install and update Bioconductor packages.
DESeq2: Analysis and visualization of RNA-Seq data.
GEOQuery: Loading of GEO dataset description.
openxlsx: reading and writing of Excel files.
ggplot2: advanced plotting package.
RColorBrewer: create color gradients for plots.
pheatmap: drawing heatmaps and clustering data.
RCy3: Vizualize, analyze and explore networks using Cytoscape via R.

3.1 setup the working directory

Please download the material and save into your folder before starting (OneDrive link).

GSE88798, saving these files into a folder named “GSE88798”.
BioPlex.3.0_edge.tsv
gordon_edge.csv
gordon_node.csv

<your path> can be absolute or relative folder path, and remember to change the path to your own path.

# clear all elements from environment
rm(list = ls())
# set working directory
##########################3
# setwd(<your path>)
##########################3
setwd("~/OneDrive/公家用/PraktikumKurse/")

3.2 Installation of required R packages

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager", dependencies = TRUE, repos = "http://cran.r-project.org")
}
library(BiocManager)
packs <- c("DESeq2", "GEOquery", "openxlsx", "ggplot2", "RColorBrewer", "pheatmap", "RCy3", "apeglm", "igraph")
for (i in packs) {
  if (!requireNamespace(i, quietly = TRUE)) {
    BiocManager::install(i)
  }
}

Check what package(s) is loaded

(.packages())

## [1] "BiocManager" "stats"       "graphics"    "grDevices"   "utils"      
## [6] "datasets"    "methods"     "base"

3.3 load packages

For the R programming environment, a lot of packages are available, which contain specialized functions for working with different dataset, conducting special analyses or producing special plots.

Here shows how to load some packages in R.

library(DESeq2)

## Warning: package 'DESeq2' was built under R version 4.1.1

## Warning: package 'S4Vectors' was built under R version 4.1.2

## Warning: package 'BiocGenerics' was built under R version 4.1.1

## Warning: package 'IRanges' was built under R version 4.1.1

## Warning: package 'GenomicRanges' was built under R version 4.1.2

## Warning: package 'GenomeInfoDb' was built under R version 4.1.1

## Warning: package 'SummarizedExperiment' was built under R version 4.1.1

## Warning: package 'MatrixGenerics' was built under R version 4.1.1

## Warning: package 'Biobase' was built under R version 4.1.1

library(GEOquery)

## Warning: package 'GEOquery' was built under R version 4.1.2

library(openxlsx)
library(ggplot2)
library(RColorBrewer)
library(pheatmap)

4 Basic Analysis of RNA-Seq dataset

In this R session, the RNA-Seq dataset GSE88798 from GEO is loaded for basic analysis and visualization. First a subset of the experiment is selected: Treatment of Col-0 with Mock or Pto_AvrRpm1 at time point 1 hour after inoculation. This dataset is used to identify genes with a significant differential expression visualize the counts of the most significant gene in three repeats in mock vs. pathogen treatment.

The following command load the RNA-seq count matrix and adjust the format for using it with the DESeq2 package

gseData <- read.table("GSE88798/GSE88798_ReadCountTable_M001_348.txt", header = T, sep = "\t")
str(gseData) # compactly display the structure of an arbitrary R object
head(gseData) # returns the first parts of a vector, matrix, table, data frame or function.
?head # asking for help (shown R documentation)
gseData[c(1:6), c(1:6)] # returns the first 6 rows and 6 columns.

4.1 Adjust the variable gseData

gseLocus <- gseData[, 1] # get locus names from first column
gseData <- gseData[, -1] # remove locus name column to be compatible with subsequent analysis
row.names(gseData) <- gseLocus # set locus names as row names
colnames(gseData) <- NULL # clear column names

The table containing the read counts does not include any description of the experiments, except the experiment identifier. The experiment description is in the file “*_series_matrix.txt”. You can open this file with a text editor to see the structure of the file.

A convenient way to load this file is using the command getGEO from the GEOquery package: read the series matrix, which contains detailed information on each experiment

https://bioconductor.org/packages/devel/bioc/vignettes/GEOquery/inst/doc/GEOquery.html#getting-gse-series-matrix-files-as-an-expressionset

gseTable <- getGEO(filename = "GSE88798/GSE88798_series_matrix.txt")
gseTable # show element

## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 0 features, 366 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: GSM2347874 GSM2347875 ... GSM2359638 (366 total)
##   varLabels: title geo_accession ... treatment:ch1 (48 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation: GPL17639

Übung: how about try to inspect the GEO table structure?

str(gseTable)

The variable gseTable is an ExpressionSet object. To work with the information contained in this object, the relevant data have to be extracted into a simple variable

# extract relevant information: genotype, replicate, time point and treatment for each experiment
myColData <- phenoData(gseTable) # extract metadata using the function phenoData
myColData

## An object of class 'AnnotatedDataFrame'
##   sampleNames: GSM2347874 GSM2347875 ... GSM2359638 (366 total)
##   varLabels: title geo_accession ... treatment:ch1 (48 total)
##   varMetadata: labelDescription

When you print out the variable myColData, then you see that it contains 366 samples, but the count table contains only 348 samples. And we are only interested in the sample properties genotype, replicate, time and treatment. So this part is now extracted from all metadata:

4.2 extract dataset for downstream analysis

# extract metadata for experiment 1 to 348, which are in
# GSE88798_ReadCountTable_M001_348.txt and the four columns genotype, replicate, time, treatment
myColData <- pData(myColData[1:348, c(44, 45, 46, 48)])
str(myColData)

## 'data.frame':    348 obs. of  4 variables:
##  $ genotype:ch1 : chr  "Col-0" "Col-0" "Col-0" "Col-0" ...
##  $ replicate:ch1: chr  "#1" "#1" "#1" "#2" ...
##  $ time:ch1     : chr  "1 hpi" "1 hpi" "1 hpi" "1 hpi" ...
##  $ treatment:ch1: chr  "Mock" "Pto DC3000" "Pto AvrRpt2" "Mock" ...

# set new colnames
colnames(myColData) <- c("genotype", "replicate", "time", "treatment")
head(myColData)

##            genotype replicate  time   treatment
## GSM2347874    Col-0        #1 1 hpi        Mock
## GSM2347875    Col-0        #1 1 hpi  Pto DC3000
## GSM2347876    Col-0        #1 1 hpi Pto AvrRpt2
## GSM2347877    Col-0        #2 1 hpi        Mock
## GSM2347878    Col-0        #2 1 hpi  Pto DC3000
## GSM2347879    Col-0        #2 1 hpi Pto AvrRpt2

4.3 data cleanup

To simplify the dataset and avoid problems unnecessary white spaces and unwanted special character are removed.

# adjust the values in the series matrix to be more intuitive
unique(myColData$replicate) # show a unique set of entries in replicate column

## [1] "#1" "#2" "#3"

myColData$replicate <- gsub("^ ", "", myColData$replicate) # remove leading whitespace
myColData$replicate <- gsub("#", "", myColData$replicate) # remove # character
myColData$replicate <- as.numeric(myColData$replicate) # convert from string to number
unique(myColData$time) # do the same for column time

## [1] "1 hpi"  "3 hpi"  "4 hpi"  "6 hpi"  "9 hpi"  "12 hpi" "16 hpi" "24 hpi"

myColData$time <- gsub("^ ", "", myColData$time)
myColData$time <- gsub(" hpi", "", myColData$time)
myColData$time <- as.numeric(myColData$time)
unique(myColData$genotype) # do the same for column genotype

## [1] "Col-0"                       "sid2-2"                     
## [3] "dde2-2 ein2-1 pad4-1 sid2-2" "pad4-1"                     
## [5] "pad4-1 sid2-2"               "rpm1-3 rpad4-1 sid2-22 101C"

myColData$genotype <- gsub("^ ", "", myColData$genotype)
myColData$genotype <- gsub(" ", "_", myColData$genotype)
unique(myColData$treatment) # do the same for column treatment

## [1] "Mock"        "Pto DC3000"  "Pto AvrRpt2" "Pto AvrRpm1"

myColData$treatment <- gsub("^ ", "", myColData$treatment)
myColData$treatment <- gsub(" ", "_", myColData$treatment)

4.4 data filtering

In this analysis we are only interested in the genes, which are differentially expressed between mock treatment and treatment with Pseudomonas syringae pv. tomato (Pto) carrying a vector with AvrRpm1 in the accession Col-0 at the time point 1 hour after treatment. Therefore we have to filter the metadata for the accession Col-0, for the treatments Mock or Pto_AvrRpm1 and the time point 1.

# select only Col-0 and Mock + Pto_AvrRpm1 treatment at time point 1 hpi
# find index of these positions
selectionIndex <-
  which(
    myColData$genotype == "Col-0" &
      (
        myColData$treatment == "Mock" |
          myColData$treatment == "Pto_AvrRpm1"
      ) & (myColData$time == 1)
  )
# extract metadata
myColData <- myColData[selectionIndex, ]
# extract counts
gseData <- gseData[, selectionIndex]

4.5 the R factor

In R there is a special data type called factor (beside the usual ones like int, string, …), which is extensively used. Here the function for differential expression analysis expects the metadata to be of data type factor. Because the metadata are stored as strings, they have to be converted to factors:

# convert all metadata to factors as this is needed by the differential expression function
myColData$time <- as.factor(myColData$time)
myColData$replicate <- as.factor(myColData$replicate)
myColData$genotype <- as.factor(myColData$genotype)
myColData$treatment <- as.factor(myColData$treatment)
# check colnames / rownames for identity and adjust them
all(rownames(myColData) %in% colnames(gseData))

## [1] FALSE

colnames(gseData) <- rownames(myColData)
all(rownames(myColData) %in% colnames(gseData))

## [1] TRUE

4.6 Construct DEG object

Use now the prepared count matrix and the metadata matrix to create a DESeq object, which is needed for the differential expression analysis. The important setting is the design parameter, which has to be set to “~ treatment” as we want to calculate a p-value, which is dependent on the treatment of the plant.

# make deseq object for differential expression test
ddsMat <- DESeqDataSetFromMatrix(countData = gseData,
                                 colData = myColData,
                                 design = ~ treatment)

Set Mock as the untreated condition.

# relevel conditions, set "Mock" as reference and remove unnessecary levels
levels(ddsMat$treatment) # show levels

## [1] "Mock"        "Pto_AvrRpm1"

ddsMat$treatment <- relevel(ddsMat$treatment, ref = "Mock")
ddsMat$treatment <- droplevels(ddsMat$treatment)

Remove genes, which have almost no counts. This accelerates the calculation.

# Pre-filtering, remove genes having less than 10 counts over all experiments
keep <- rowSums(counts(ddsMat)) >= 10
ddsMat <- ddsMat[keep, ]

4.7 DEG results

Conduct the differential expression analysis. Count normalization is done automatically.

# Differential expression analysis
ddsMat <- DESeq(ddsMat, parallel = FALSE) # run differential expression analysis

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

res <- results(ddsMat) # extract results from DESeq object
res

## log2 fold change (MLE): treatment Pto AvrRpm1 vs Mock 
## Wald test p-value: treatment Pto AvrRpm1 vs Mock 
## DataFrame with 20255 rows and 6 columns
##             baseMean log2FoldChange     lfcSE      stat      pvalue        padj
##            <numeric>      <numeric> <numeric> <numeric>   <numeric>   <numeric>
## AT1G01010  378.65099      0.0507979  0.614377  0.082682  0.93410443   0.9721451
## AT1G01020  273.86578     -0.3695763  0.283748 -1.302483  0.19275136   0.4569515
## AT1G01030  207.58372     -1.2485585  0.478459 -2.609542  0.00906636   0.0823238
## AT1G01040 1048.62226     -0.3100009  0.354170 -0.875289  0.38141661   0.6397881
## AT1G01046    5.31325      0.3726277  1.105924  0.336938  0.73616377   0.8756635
## ...              ...            ...       ...       ...         ...         ...
## AT5G67600    777.605      1.1828501  0.743836   1.59020 1.11789e-01 0.347546604
## AT5G67610    251.943     -0.0417828  0.286693  -0.14574 8.84126e-01 0.952025704
## AT5G67620    161.053     -3.3379630  0.701943  -4.75532 1.98136e-06 0.000176331
## AT5G67630    685.906     -1.2305383  0.534264  -2.30324 2.12653e-02 0.138499846
## AT5G67640    138.371     -2.1042483  0.673674  -3.12354 1.78690e-03 0.027777002

5 Plot the calculated log2 fold changes vs normalized counts.

5.1 original

# Exploring results, plot log2 fold change vs. mean normalized counts
plotMA(res, ylim = c(-3, 3))

Figure 3. Gene expression (fold change) distribution. Noise not removed yet.

5.2 Remove noise

resultsNames(ddsMat)

## [1] "Intercept"                     "treatment_Pto_AvrRpm1_vs_Mock"

resLFC <- lfcShrink(ddsMat, coef = "treatment_Pto_AvrRpm1_vs_Mock")

## using 'apeglm' for LFC shrinkage. If used in published research, please cite:
##     Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
##     sequence count data: removing the noise and preserving large differences.
##     Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895

resLFC

## log2 fold change (MAP): treatment Pto AvrRpm1 vs Mock 
## Wald test p-value: treatment Pto AvrRpm1 vs Mock 
## DataFrame with 20255 rows and 5 columns
##             baseMean log2FoldChange     lfcSE      pvalue        padj
##            <numeric>      <numeric> <numeric>   <numeric>   <numeric>
## AT1G01010  378.65099      0.0211266  0.383443  0.93410443   0.9721451
## AT1G01020  273.86578     -0.2868327  0.258857  0.19275136   0.4569515
## AT1G01030  207.58372     -0.9337151  0.506783  0.00906636   0.0823238
## AT1G01040 1048.62226     -0.2089155  0.299523  0.38141661   0.6397881
## AT1G01046    5.31325      0.0628926  0.452467  0.73616377   0.8756635
## ...              ...            ...       ...         ...         ...
## AT5G67600    777.605      0.4481358  0.586555 1.11789e-01 0.347546604
## AT5G67610    251.943     -0.0306231  0.247692 8.84126e-01 0.952025704
## AT5G67620    161.053     -3.0234911  0.740305 1.98136e-06 0.000176331
## AT5G67630    685.906     -0.8192357  0.563294 2.12653e-02 0.138499846
## AT5G67640    138.371     -1.6241491  0.765341 1.78690e-03 0.027777002

# plot without noise
plotMA(resLFC, ylim = c(-5,5))
# identify points interactively, click on points, then press "Finish"
idx <- identify(resLFC$baseMean, resLFC$log2FoldChange)

rownames(res)[idx]

## character(0)

Figure 4. Gene expression (fold change) distribution. Noise cleanup.

5.3 DEGs

Identify the most significantly upregulated gene and plot the counts in all three repeats for mock treatment and Pto_AvrRpm1 treatment.

# plot counts for most significant upregulated gene
upIndex <- which(res$log2FoldChange > 0) # find upregulated genes
myIndex <- upIndex[which.min(res$padj[upIndex])] # find most significant upregulated gene

5.4 extract read counts for selected gene

fiss <- plotCounts(ddsMat, myIndex,
                   intgroup = c("treatment"), returnData = TRUE)
fiss

##                 count   treatment
## GSM2347874  12.889481        Mock
## GSM2347877   6.476583        Mock
## GSM2347880   4.562768        Mock
## GSM2348162 488.707390 Pto_AvrRpm1
## GSM2348163 523.244554 Pto_AvrRpm1
## GSM2348164 530.299487 Pto_AvrRpm1

5.5 Shown interesting gene

# create plot from read counts
g <- ggplot(fiss,
       aes(x = treatment, y = count, color = treatment, group = treatment)) +
        geom_point(position = position_jitter(w = 0.1, h = 0)) + geom_smooth(se = FALSE, method = "loess") + scale_y_log10(breaks = c(10, 50, 100, 500, 1000, 5000)) +
       ggtitle(rownames(res)[myIndex])
# print plot
print(g)

## `geom_smooth()` using formula 'y ~ x'

Figure 5. The gene with highest up-regulated fold-change under Pto AvrRpm1 treatment.

Saving to PDF

# print plot to pdf file
pdf("Most significant gene.pdf", width = 8, height = 6)
print(g)
dev.off()

5.6 Statistical analysis of the differentially expressed genes

Write out the significantly differentially expressed genes

# Identification of differentially expressed genes
summary(res)

## 
## out of 20255 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 780, 3.9%
## LFC < 0 (down)     : 1620, 8%
## outliers [1]       : 246, 1.2%
## low counts [2]     : 786, 3.9%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

# set P-Value < 0.05
res05 <- results(ddsMat, alpha = 0.05)
summary(res05)

## 
## out of 20255 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 451, 2.2%
## LFC < 0 (down)     : 1135, 5.6%
## outliers [1]       : 246, 1.2%
## low counts [2]     : 786, 3.9%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

# write results to excel file
myData <- data.frame(res[which(res$padj < 0.05), ])
myData <- myData[order(myData$padj), ]

Save differentially expressed genes into XLSX file.

write.xlsx(myData, "Significantly DE Genes.xlsx")

6 TASK 1

For the 10 upregulated and most significant genes, check on Arabidopsis.org, if they are related to pathogens / immunity

6.1 Data transformations and visualization for further evaluation

Transformation of counts, to account for different sequencing depth. This allows to compare different experiments.

# Variance stabilizing transformation
vsd <- vst(ddsMat, blind = FALSE)
head(assay(vsd), 5)

7 Visulization of the results

7.1 normalized counts

Plot the normalized values RNA-Seq counts to see if the results are consistent.

# Data quality assessment by sample clustering and visualization
# identify the 20 most significantly upregulated genes
upIndex <- which(res$log2FoldChange > 0)
myIndex <- upIndex[order(res$padj[upIndex])[1:20]]

df <- as.data.frame(colData(ddsMat)[, c("treatment")])
rownames(df) <- colnames(vsd)
colnames(df) <- c("treatment")

pheatmap(assay(vsd)[myIndex,], cluster_rows = FALSE, show_rownames = T,
         cluster_cols = FALSE, annotation_col = df)

Figure 6. Heatmap of top 20 upregulated genes under treatment.

7.2 data similarity

For the heatmap, the Euclidean distance (or distance measures) is calculated for all pairwise combinations of the experiment. The distance is used to cluster the experiments and identify their similarity. Normally one would expect, that experiments with the same treatment have a low distance and cluster together, whereas experiments with different treatments should have a bigger distance.

sampleDists <- dist(t(assay(vsd))) # calc euclidean distance between all 6 experiments
sampleDistMatrix <- as.matrix(sampleDists) # distance matrix
rownames(sampleDistMatrix) <- paste(vsd$treatment, sep = "-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors)

Figure 7. Heatmap of sample similarity.

7.3 PCA

Prinicpal Compnent Analysis (PCA): A PCA is an unsupervised machine learning method to identify entities, which belong together or to separate entities, which have different properties. In the following figure the PCA results of three different iris species are shown. It is tested, if these three species can be separated on basis of the length and width of sepals and petals. You can see that Iris setosa can be separated well from the other two species, but Iris versicolor and Iris virginica are strongly overlapping and cannot be separated using two principal components.

In our case the PCA is used to check, if the experiments with the same treatment are near to each other and experiments with a different treatment are distant to each other.

pcaData <- plotPCA(vsd, intgroup = c("treatment"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
q <- ggplot(pcaData, aes(PC1, PC2, color = treatment)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed()
print(q)

Figure 8. PCA plot of sample profiles.

8 Task 2

Discuss the 3 repeats of the mock and AvrRpm1 treatment on the basis of the PCA results and the distance plot. Would you expect inconsistencies between the repeats only looking at the VSD values plot or the dot plot of the most significantly upregulated gene? What is your suggestion for the experimenter?

9 Bonus

How to cooperate (access and control) between R and Cytoscape? We use RCy3 package here. The hot topic SARS-CoV-2 is studying here. We will have the viral-human protein interaction network, which is from Gordon et al., 2020

9.1 Create Network

library(igraph) # for network establishment

## Warning: package 'igraph' was built under R version 4.1.2

library(RCy3) # for cooperation

## Warning: package 'RCy3' was built under R version 4.1.2

sars2 <- read.csv("gordon_edge.csv", header = TRUE)
sars2_node <- read.csv("gordon_node.csv", header = T)
sars2_network <- graph_from_data_frame(sars2, directed = FALSE)
V(sars2_network)$source <- sars2_node$group
V(sars2_network)$color <- ifelse(V(sars2_network)$source == "human", "blue", "red")
plot(sars2_network) # basic plot

plot(sars2_network, layout = layout_nicely, vertex.size = 3, vertex.label.cex = .4, vertex.label.color = "black", margin = -0.5) # modified igraph object

9.2 cooperate with Cytoscape

Now open your Cytoscape application, and ready to send command from R to Cytoscape

createNetworkFromIgraph(sars2_network, title = "SARS-CoV-2", collection = "Virus")

style <- "covid19"
defaults <- list(NODE_SHAPE = "ellipse",
                 NODE_SIZE = 30,
                 EDGE_TRANSPARENCY = 120,
                 NODE_LABEL_POSITION = "W, E, c, 0.00, 0.00",
                 NODE_FILL_COLOR = "#00bfff",
                 NODE_Border_Width = 0)
nodeLabels <- mapVisualProperty(visual.prop = "node label",
                table.column = "name",
                mapping.type = "p")
nodeFills <- mapVisualProperty(visual.prop = "node fill color",
             table.column = "source",
             mapping.type = "d",
             table.column.values = "virus",
             visual.prop.values = "#ff0000")
createVisualStyle(style, defaults, list(nodeLabels, nodeFills))
setVisualStyle(style)

integrate with other network

To know the virus associated human protein interactions, we integrate the BioPlex interactome.

bioplex <- read.table("BioPlex.3.0_edge.tsv", sep = "\t", header = T)
bioplex <- bioplex[, c(5, 6)]
bioplex_network <- graph_from_data_frame(bioplex, directed = FALSE)
bioplex_network <- simplify(bioplex_network, remove.loops = FALSE) # remove duplicate interactions if exists

sars2_bioplex <- induced_subgraph(bioplex_network, names(V(sars2_network))[names(V(sars2_network)) %in% names(V(bioplex_network))])
createNetworkFromIgraph(sars2_bioplex, title = "SARS2 in Bioplex", collection = "human") # might suffer an error message **https://github.com/cytoscape/RCy3/issues/98**

merge two networks

merge2 <- sars2_bioplex %u% sars2_network
createNetworkFromIgraph(merge2, title = "SARS2 union SARS2 in Bioplex", collection = "merged network")
setVisualStyle(style, network = "SARS2 union SARS2 in Bioplex")

System information

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] RCy3_2.14.1                 igraph_1.2.11              
##  [3] pheatmap_1.0.12             RColorBrewer_1.1-2         
##  [5] ggplot2_3.3.5               openxlsx_4.2.5             
##  [7] GEOquery_2.62.2             DESeq2_1.34.0              
##  [9] SummarizedExperiment_1.24.0 Biobase_2.54.0             
## [11] MatrixGenerics_1.6.0        matrixStats_0.61.0         
## [13] GenomicRanges_1.46.1        GenomeInfoDb_1.30.0        
## [15] IRanges_2.28.0              S4Vectors_0.32.3           
## [17] BiocGenerics_0.40.0         BiocManager_1.30.16        
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.0-2       ellipsis_0.3.2         IRdisplay_1.1         
##   [4] XVector_0.34.0         base64enc_0.1-3        fs_1.5.2              
##   [7] farver_2.1.0           bit64_4.0.5            mvtnorm_1.1-3         
##  [10] AnnotationDbi_1.56.2   fansi_1.0.0            apeglm_1.16.0         
##  [13] xml2_1.3.3             R.methodsS3_1.8.1      splines_4.1.0         
##  [16] cachem_1.0.6           geneplotter_1.72.0     knitr_1.37            
##  [19] IRkernel_1.3           jsonlite_1.7.2         annotate_1.72.0       
##  [22] R.oo_1.24.0            png_0.1-7              graph_1.72.0          
##  [25] readr_2.1.1            compiler_4.1.0         httr_1.4.2            
##  [28] backports_1.4.1        assertthat_0.2.1       Matrix_1.4-0          
##  [31] fastmap_1.1.0          limma_3.50.0           htmltools_0.5.2       
##  [34] tools_4.1.0            coda_0.19-4            gtable_0.3.0          
##  [37] glue_1.6.0             GenomeInfoDbData_1.2.7 dplyr_1.0.7           
##  [40] Rcpp_1.0.8             bbmle_1.0.24           jquerylib_0.1.4       
##  [43] vctrs_0.3.8            Biostrings_2.62.0      RJSONIO_1.3-1.6       
##  [46] xfun_0.29              stringr_1.4.0          lifecycle_1.0.1       
##  [49] XML_3.99-0.8           zlibbioc_1.40.0        MASS_7.3-55           
##  [52] scales_1.1.1           hms_1.1.1              parallel_4.1.0        
##  [55] curl_4.3.2             yaml_2.2.1             memoise_2.0.1         
##  [58] emdbook_1.3.12         sass_0.4.0             bdsmatrix_1.3-4       
##  [61] stringi_1.7.6          RSQLite_2.2.9          highr_0.9             
##  [64] genefilter_1.76.0      zip_2.2.0              BiocParallel_1.28.3   
##  [67] repr_1.1.4             rlang_0.4.12           pkgconfig_2.0.3       
##  [70] bitops_1.0-7           uchardet_1.1.0         evaluate_0.14         
##  [73] lattice_0.20-45        purrr_0.3.4            labeling_0.4.2        
##  [76] bit_4.0.4              tidyselect_1.1.1       plyr_1.8.6            
##  [79] magrittr_2.0.1         R6_2.5.1               generics_0.1.1        
##  [82] base64url_1.4          pbdZMQ_0.3-6           DelayedArray_0.20.0   
##  [85] DBI_1.1.2              withr_2.4.3            pillar_1.6.4          
##  [88] survival_3.2-13        KEGGREST_1.34.0        RCurl_1.98-1.5        
##  [91] tibble_3.1.6           crayon_1.4.2           uuid_1.0-3            
##  [94] utf8_1.2.2             tzdb_0.2.0             rmarkdown_2.11        
##  [97] locfit_1.5-9.4         grid_4.1.0             data.table_1.14.2     
## [100] blob_1.2.2             digest_0.6.29          xtable_1.8-4          
## [103] tidyr_1.1.4            numDeriv_2016.8-1.1    R.utils_2.11.0        
## [106] munsell_0.5.0          bslib_0.3.1

Praktikum kurse, R-Sektion

2022.01.17 - 2022.01.28

1 Analyze RNA-Seq data with R and Cytoscape

1.1 Softwares in this course

R

RStudio

2 GUI Basics

2.1 Software Interface

R

RStudio

2.2 The environment of RStudio

3 Working environment setup

3.1 setup the working directory

3.2 Installation of required R packages

3.3 load packages

4 Basic Analysis of RNA-Seq dataset

4.1 Adjust the variable gseData

4.2 extract dataset for downstream analysis

4.3 data cleanup

4.4 data filtering

4.5 the R factor

4.6 Construct DEG object

4.7 DEG results

5 Plot the calculated log2 fold changes vs normalized counts.

5.1 original

5.2 Remove noise

5.3 DEGs

5.4 extract read counts for selected gene

5.5 Shown interesting gene

Saving to PDF

5.6 Statistical analysis of the differentially expressed genes

6 TASK 1

6.1 Data transformations and visualization for further evaluation

7 Visulization of the results

7.1 normalized counts

7.2 data similarity

7.3 PCA

8 Task 2

9 Bonus

9.1 Create Network

9.2 cooperate with Cytoscape

integrate with other network

merge two networks