-
run nextflow rnaseq
#under sage ln -s /home/jhuang/Tools/nf-core-rnaseq-3.12.0/ rnaseq (rnaseq) nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results_GRCh38 --genome GRCh38 --with_umi --umitools_extract_method regex --umitools_bc_pattern '^(?P.{12}).*' -profile docker -resume --max_cpus 54 --max_memory 120.GB --max_time 2400.h --save_align_intermeds --save_unaligned --save_reference --aligner star_salmon --pseudo_aligner salmon --umitools_grouping_method unique #Debug the following error: added "--minAssignedFrags 0 \\" to modules/nf-core/salmon/quant/main.nf option "salmon quant" and added "--min_mapped_reads 0" in the nextflow command (rnaseq) nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results_chrHsv1 --fasta chrHsv1_s17.fasta --gtf chrHsv1_s17.gtf --with_umi --umitools_extract_method regex --umitools_bc_pattern '^(?P.{12}).*' --umitools_dedup_stats -profile test_full -resume --max_memory 256.GB --max_time 2400.h --save_reference --aligner star_salmon --gtf_extra_attributes gene_id --gtf_group_features transcript_id --featurecounts_group_type gene_id --featurecounts_feature_type transcript --skip_rseqc --skip_dupradar --skip_preseq --skip_biotype_qc --skip_deseq2_qc --skip_multiqc --min_mapped_reads 0 #TODO: why a lot of reads were removed due to the too_short? -
import data and pca-plot
# Import the required libraries library("AnnotationDbi") library("clusterProfiler") library("ReactomePA") library(gplots) library(tximport) library(DESeq2) library("org.Hs.eg.db") library(dplyr) library(tidyverse) setwd("~/DATA/Data_Manja_RNAseq_Organoids_Virus/results_chrHsv1_downstream/star_salmon") # Define paths to your Salmon output quantification files files <- c("control_r1" = "./control_r1/quant.sf", "control_r2" = "./control_r2/quant.sf", "HSV.d2_r1" = "./HSV.d2_r1/quant.sf", "HSV.d2_r2" = "./HSV.d2_r2/quant.sf", "HSV.d4_r1" = "./HSV.d4_r1/quant.sf", "HSV.d4_r2" = "./HSV.d4_r2/quant.sf", "HSV.d6_r1" = "./HSV.d6_r1/quant.sf", "HSV.d6_r2" = "./HSV.d6_r2/quant.sf", "HSV.d8_r1" = "./HSV.d8_r1/quant.sf", "HSV.d8_r2" = "./HSV.d8_r2/quant.sf") # Import the transcript abundance data with tximport txi <- tximport(files, type = "salmon", txIn = TRUE, txOut = TRUE) # Define the replicates and condition of the samples replicate <- factor(c("r1", "r2", "r1", "r2", "r1", "r2", "r1", "r2", "r1", "r2")) condition <- factor(c("control", "control", "HSV.d2", "HSV.d2", "HSV.d4", "HSV.d4", "HSV.d6", "HSV.d6", "HSV.d8", "HSV.d8")) # Define the colData for DESeq2 colData <- data.frame(condition=condition, replicate=replicate, row.names=names(files)) # -- transcript-level count data (x2) -- # Create DESeqDataSet object dds <- DESeqDataSetFromTximport(txi, colData=colData, design=~condition) write.csv(counts(dds), file="transcript_counts.csv") # -- gene-level count data (x2) -- # Read in the tx2gene map from salmon_tx2gene.tsv tx2gene <- read.table("salmon_tx2gene.tsv", header=FALSE, stringsAsFactors=FALSE) # Set the column names colnames(tx2gene) <- c("transcript_id", "gene_id", "gene_name") # Remove the gene_name column if not needed tx2gene <- tx2gene[,1:2] # Import and summarize the Salmon data with tximport txi <- tximport(files, type = "salmon", tx2gene = tx2gene, txOut = FALSE) # Continue with the DESeq2 workflow as before... colData <- data.frame(condition=condition, replicate=replicate, row.names=names(files)) dds <- DESeqDataSetFromTximport(txi, colData=colData, design=~condition) #dds <- dds[rowSums(counts(dds) > 3) > 2, ] #60605-->26543 write.csv(counts(dds, normalized=FALSE), file="gene_counts.csv") # -- merge the raw counts of human and microbe -- #cat ~/DATA/Data_Manja_RNAseq_Organoids/results_GRCh38_unique/star_salmon/gene_counts.csv ~/DATA/Data_Manja_RNAseq_Organoids_Virus/results_chrHsv1_downstream/star_salmon/gene_counts.csv > merged_gene_counts.csv #DELETE the second line "","control_r1","control_r2","HSV.d2_r1","HSV.d2_r2","HSV.d4_r1","HSV.d4_r2","HSV.d6_r1","HSV.d6_r2","HSV.d8_r1","HSV.d8_r2" #~/Tools/csv2xls-0.4/csv_to_xls.py merged_gene_counts.csv -d',' -o raw_gene_counts.xls; # -- for merged analysis due to false normalization factors wenn alone analyzed on virus data -- setwd("~/DATA/Data_Manja_RNAseq_Organoids_Merged/") d.raw <- read.csv("merged_gene_counts.csv", header=TRUE, row.names=1) dds <- DESeqDataSetFromMatrix(countData=d.raw, colData=colData, design=~condition+replicate) dim(counts(dds)) head(counts(dds), 10) rld <- rlogTransformation(dds) #We don't need to run DESeq(dds) before estimateSizeFactors(dds). In fact, the typical workflow in DESeq2 is the opposite: we usually run estimateSizeFactors(dds) (and other preprocessing functions) before running the main DESeq(dds) function. #The estimateSizeFactors function is used to calculate size factors for normalization, which corrects for differences in library size (i.e., the number of read counts) between samples. This normalization step is crucial to ensure that differences in gene expression aren't merely due to differences in sequencing depth between samples. #The DESeq function, on the other hand, performs the main differential expression analysis, comparing gene expression between different conditions or groups. #So, the typical workflow is: # - Create the DESeqDataSet object. # - Use estimateSizeFactors to normalize for library size. # - (Optionally, estimate dispersion with estimateDispersions if not using the full DESeq function later.) # - Use DESeq for the differential expression analysis. # - However, it's worth noting that if you run the main DESeq function directly after creating the DESeqDataSet object, it will automatically perform the normalization (using estimateSizeFactors) and dispersion estimation steps for you. In that case, there's no need to run estimateSizeFactors separately before DESeq. # draw simple pca and heatmap library(gplots) library("RColorBrewer") #mat <- assay(rld) #mm <- model.matrix(~condition, colData(rld)) #mat <- limma::removeBatchEffect(mat, batch=rld$batch, design=mm) #assay(rld) <- mat # -- pca -- png("pca.png", 1200, 800) plotPCA(rld, intgroup=c("condition")) dev.off() # -- heatmap -- png("heatmap.png", 1200, 800) distsRL <- dist(t(assay(rld))) mat <- as.matrix(distsRL) hc <- hclust(distsRL) hmcol <- colorRampPalette(brewer.pal(9,"GnBu"))(100) heatmap.2(mat, Rowv=as.dendrogram(hc),symm=TRUE, trace="none",col = rev(hmcol), margin=c(13, 13)) dev.off() -
draw 3D PCA plots.
library(gplots) library("RColorBrewer") library(ggplot2) data <- plotPCA(rld, intgroup=c("condition", "replicate"), returnData=TRUE) write.csv(data, file="plotPCA_data.csv") #calculate all PCs including PC3 with the following codes library(genefilter) ntop <- 500 rv <- rowVars(assay(rld)) select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))] mat <- t( assay(rld)[select, ] ) pc <- prcomp(mat) pc$x[,1:3] #df_pc <- data.frame(pc$x[,1:3]) df_pc <- data.frame(pc$x) identical(rownames(data), rownames(df_pc)) #-->TRUE data$PC1 <- NULL data$PC2 <- NULL merged_df <- merge(data, df_pc, by = "row.names") #merged_df <- merged_df[, -1] row.names(merged_df) <- merged_df$Row.names merged_df$Row.names <- NULL # remove the "name" column merged_df$name <- NULL merged_df <- merged_df[, c("PC1","PC2","PC3","PC4","PC5","PC6","PC7","PC8","PC9","PC10","group","condition","replicate")] write.csv(merged_df, file="merged_df_10PCs.csv") summary(pc) #0.5333 0.2125 0.06852 #0.8026 0.09042 0.06578 #0.6577 0.1631 0.07106 #Proportion of Variance 0.6577 0.1631 0.07106 draw_3D.py #/usr/bin/convert PCA_3D.png -crop 2900x1600+250+700 PCA_3D_cropped.png # adjust proportion to real values in the following plot import plotly.graph_objects as go import pandas as pd from sklearn.decomposition import PCA import numpy as np # Read in data as a pandas dataframe #df = pd.DataFrame({ # 'PC1': [-13.999925, -12.504291, -12.443057, -13.065235, -17.316215], # 'PC2': [-1.498823, -3.342411, -6.067055, -8.205809, 3.293993], # 'PC3': [-3.335085, 15.207755, -14.725450, 15.078469, -6.917358], # 'condition': ['GFP d3', 'GFP d3', 'GFP d8', 'GFP d8', 'GFP+mCh d9/12'], # 'replicate': ['DI', 'DII', 'DI', 'DII', 'DI'] #}) df = pd.read_csv('merged_df_10PCs.csv', index_col=0, header=0) df['condition'] = df['condition'].replace("control", "control") df['condition'] = df['condition'].replace("HSV.d2", "day 2") df['condition'] = df['condition'].replace("HSV.d4", "day 4") df['condition'] = df['condition'].replace("HSV.d6", "day 6") df['condition'] = df['condition'].replace("HSV.d8", "day 8") # Fit PCA model to reduce data dimensions to 3 pca = PCA(n_components=3) pca.fit(df.iloc[:, :-3]) X_reduced = pca.transform(df.iloc[:, :-3]) # Get variance ratios explained_variance_ratio = pca.explained_variance_ratio_ # Add reduced data back to dataframe df['PC1'] = X_reduced[:, 0] df['PC2'] = X_reduced[:, 1] df['PC3'] = X_reduced[:, 2] # Create PCA plot with 3D scatter fig = go.Figure() ##ff7f00 condition_color_map = { 'control': 'rgb(100, 100, 100)', 'day 2': '#33a02c', 'day 4': '#1f78b4', 'day 6': '#e31a1c', 'day 8': 'magenta' } replicate_symbol_map = {'r1': 'circle', 'r2': 'diamond'} for replicate, replicate_symbol in replicate_symbol_map.items(): for condition, condition_color in condition_color_map.items(): mask = (df['condition'] == condition) & (df['replicate'] == replicate) fig.add_trace(go.Scatter3d(x=df.loc[mask, 'PC1'], y=df.loc[mask, 'PC2'], z=df.loc[mask, 'PC3'], mode='markers', name=f'{condition}' if replicate == 'r1' else None, legendgroup=f'{condition}', showlegend=True if replicate == 'r1' else False, marker=dict(size=6 if replicate_symbol in ['diamond'] else 10, opacity=0.8, color=condition_color, symbol=replicate_symbol))) for replicate, replicate_symbol in replicate_symbol_map.items(): fig.add_trace(go.Scatter3d(x=[None], y=[None], z=[None], mode='markers', name=replicate, legendgroup=f'{replicate}', showlegend=True, marker=dict(size=10, opacity=1, color='black', symbol=replicate_symbol), hoverinfo='none')) # Annotations for the legend blocks #TODO: calculate the PC values. #TODO: adjust the axis length according to the actual size of axis! fig.update_layout( annotations=[ dict(x=1.1, y=1.0, xref='paper', yref='paper', showarrow=False, text='Condition', font=dict(size=15)), dict(x=1.1, y=0.6, xref='paper', yref='paper', showarrow=False, text='Replicate', font=dict(size=15)) ], scene=dict( #aspectmode='cube', #xaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title='PC1: 53% v.', scaleratio=0.53), #yaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title='PC2: 21% v.', scaleratio=0.21), #zaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title='PC3: 7% variance', scaleratio=0.07), aspectmode='manual', aspectratio=dict(x=explained_variance_ratio[0], y=explained_variance_ratio[1], z=explained_variance_ratio[2]), xaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title='PC1: 53% v.', range=[min(df['PC1']), max(df['PC1'])]), yaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title='PC2: 21% v.', range=[min(df['PC2']), max(df['PC2'])]), zaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title='PC3: 7% variance', range=[min(df['PC3']), max(df['PC3'])]), bgcolor='white' ), margin=dict(l=5, r=5, b=5, t=0) # Adjust the margins to prevent clipping of axis titles ) #fig.show() fig.write_image("fig1.svg") fig.update_layout( annotations=[ dict(x=1.1, y=1.0, xref='paper', yref='paper', showarrow=False, text='Condition', font=dict(size=15)), dict(x=1.1, y=0.6, xref='paper', yref='paper', showarrow=False, text='Donor', font=dict(size=15)), dict(x=1.08, y=0.2, xref='paper', yref='paper', showarrow=False, text=f'PC3: {explained_variance_ratio[2]*100:.2f}% v.', font=dict(size=15), textangle=-90) ], scene=dict( aspectmode='manual', aspectratio=dict(x=explained_variance_ratio[0]*2, y=explained_variance_ratio[1]*2, z=explained_variance_ratio[2]*2), #, range=[min(df['PC1']), max(df['PC1'])] #, range=[min(df['PC2']), max(df['PC2'])] #, range=[min(df['PC3']), max(df['PC3'])] xaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title=f'PC1: {explained_variance_ratio[0]*100:.2f}% variance'), yaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title=f'PC2: {explained_variance_ratio[1]*100:.2f}% v.'), zaxis=dict(gridcolor='black', backgroundcolor='white', zerolinecolor='black', title=''), bgcolor='white' ), margin=dict(l=0, r=0, b=0, t=0) # Adjust the margins to prevent clipping of axis titles ) fig.write_image("PCA_3D.svg") #/usr/bin/convert g235.png -crop 3250x1680+1+750 PCA_3D_.png -
(optional) estimate size factors and dispersion values.
#Size Factors: These are used to normalize the read counts across different samples. The size factor for a sample accounts for differences in sequencing depth (i.e., the total number of reads) and other technical biases between samples. After normalization with size factors, the counts should be comparable across samples. Size factors are usually calculated in a way that they reflect the median or mean ratio of gene expression levels between samples, assuming that most genes are not differentially expressed. #Dispersion: This refers to the variability or spread of gene expression measurements. In RNA-seq data analysis, each gene has its own dispersion value, which reflects how much the counts for that gene vary between different samples, more than what would be expected just due to the Poisson variation inherent in counting. Dispersion is important for accurately modeling the data and for detecting differentially expressed genes. #So in summary, size factors are specific to samples (used to make counts comparable across samples), and dispersion values are specific to genes (reflecting variability in gene expression). sizeFactors(dds) #NULL # Estimate size factors dds <- estimateSizeFactors(dds) # Estimate dispersions dds <- estimateDispersions(dds) #> sizeFactors(dds) #control_r1 control_r2 HSV.d2_r1 HSV.d2_r2 HSV.d4_r1 HSV.d4_r2 HSV.d6_r1 #2.3282468 2.0251928 1.8036883 1.3767551 0.9341929 1.0911693 0.5454526 #HSV.d6_r2 HSV.d8_r1 HSV.d8_r2 #0.4604461 0.5799834 0.6803681 # If alone with virus data, the following BUG occured: #Still NULL --> BUG --> using manual calculation method for sizeFactor calculation! HeLa_TO_r1 HeLa_TO_r2 0.9978755 1.1092227 data.frame(genes = rownames(dds), dispersions = dispersions(dds)) #Given the raw counts, the control_r1 and control_r2 samples seem to have a much lower sequencing depth (total read count) than the other samples. Therefore, when normalization methods are applied, the normalization factors for these control samples will be relatively high, boosting the normalized counts. 1/0.9978755=1.002129023 1/1.1092227= #bamCoverage --bam ../markDuplicates/${sample}Aligned.sortedByCoord.out.bam -o ${sample}_norm.bw --binSize 10 --scaleFactor --effectiveGenomeSize 2864785220 bamCoverage --bam ../markDuplicates/HeLa_TO_r1Aligned.sortedByCoord.out.markDups.bam -o HeLa_TO_r1.bw --binSize 10 --scaleFactor 1.002129023 --effectiveGenomeSize 2864785220 bamCoverage --bam ../markDuplicates/HeLa_TO_r2Aligned.sortedByCoord.out.markDups.bam -o HeLa_TO_r2.bw --binSize 10 --scaleFactor 0.901532217 --effectiveGenomeSize 2864785220 raw_counts <- counts(dds) normalized_counts <- counts(dds, normalized=TRUE) #write.table(raw_counts, file="raw_counts.txt", sep="\t", quote=F, col.names=NA) #write.table(normalized_counts, file="normalized_counts.txt", sep="\t", quote=F, col.names=NA) #convert bam to bigwig using deepTools by feeding inverse of DESeq’s size Factor estimSf <- function (cds){ # Get the count matrix cts <- counts(cds) # Compute the geometric mean geomMean <- function(x) prod(x)^(1/length(x)) # Compute the geometric mean over the line gm.mean <- apply(cts, 1, geomMean) # Zero values are set to NA (avoid subsequentcdsdivision by 0) gm.mean[gm.mean == 0] <- NA # Divide each line by its corresponding geometric mean # sweep(x, MARGIN, STATS, FUN = "-", check.margin = TRUE, ...) # MARGIN: 1 or 2 (line or columns) # STATS: a vector of length nrow(x) or ncol(x), depending on MARGIN # FUN: the function to be applied cts <- sweep(cts, 1, gm.mean, FUN="/") # Compute the median over the columns med <- apply(cts, 2, median, na.rm=TRUE) # Return the scaling factor return(med) } #https://dputhier.github.io/ASG/practicals/rnaseq_diff_Snf2/rnaseq_diff_Snf2.html #http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#data-transformations-and-visualization #https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html #https://hbctraining.github.io/DGE_workshop/lessons/04_DGE_DESeq2_analysis.html #https://genviz.org/module-04-expression/0004/02/01/DifferentialExpression/ #DESeq2’s median of ratios [1] #EdgeR’s trimmed mean of M values (TMM) [2] #http://www.nathalievialaneix.eu/doc/html/TP1_normalization.html #very good website! test_normcount <- sweep(raw_counts, 2, sizeFactors(dds), "/") sum(test_normcount != normalized_counts) -
compare the normalization methods
#The provided code indeed simulates the normalization method used by DESeq, which is known as the "Relative Log Expression (RLE)" normalization. The core idea behind this method is to scale the counts in each sample (column) by a size factor derived from the geometric means across all samples. #Here's a step-by-step breakdown of the code: # 1. estimSf function: # - The counts matrix is retrieved from the DESeq dataset. # - For each gene (row), the geometric mean across all samples is computed. # - Counts are divided by their respective gene's geometric mean. # - For each sample (column), the median of these ratios is computed. This median serves as the size factor for the sample. # 2. Once size factors are computed, the counts in the original matrix are then divided by these size factors to get the normalized counts. ### cds is a countDataset estimSf <- function (cds){ # Get the count matrix cts <- counts(cds) # Compute the geometric mean geomMean <- function(x) prod(x)^(1/length(x)) # Compute the geometric mean over the line gm.mean <- apply(cts, 1, geomMean) # Zero values are set to NA (avoid subsequentcdsdivision by 0) gm.mean[gm.mean == 0] <- NA # Divide each line by its corresponding geometric mean # sweep(x, MARGIN, STATS, FUN = "-", check.margin = TRUE, ...) # MARGIN: 1 or 2 (line or columns) # STATS: a vector of length nrow(x) or ncol(x), depending on MARGIN # FUN: the function to be applied cts <- sweep(cts, 1, gm.mean, FUN="/") # Compute the median over the columns med <- apply(cts, 2, median, na.rm=TRUE) # Return the scaling factor return(med) } #https://dputhier.github.io/ASG/practicals/rnaseq_diff_Snf2/rnaseq_diff_Snf2.html #http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#data-transformations-and-visualization #https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html #https://hbctraining.github.io/DGE_workshop/lessons/04_DGE_DESeq2_analysis.html #https://genviz.org/module-04-expression/0004/02/01/DifferentialExpression/ #DESeq2’s median of ratios [1] #EdgeR’s trimmed mean of M values (TMM) [2] #http://www.nathalievialaneix.eu/doc/html/TP1_normalization.html #very good website! test_normcount <- sweep(raw_counts, 2, sizeFactors(dds), "/") sum(test_normcount != normalized_counts) # round(estimSf(dds),6) is manually calculated sizeFactors of head(estimSf(dds)) all(round(estimSf(dds),6) == round(sizeFactors(dds), 6)) ## Checking the normalization png("normalization_small.png", width=800, height=600) epsilon <- 1 # pseudo-count to avoid problems with log(0) par(mfrow=c(1,2),cex.lab=0.7) boxplot(log2(raw_counts+epsilon), cex.axis=0.7, las=1, xlab="log2(raw counts)", horizontal=TRUE, main="Raw counts") boxplot(log2(normalized_counts+epsilon), cex.axis=0.7, las=1, xlab="log2(normalized counts)", horizontal=TRUE, main="Normalized counts") #boxplot(log10(assays(dds)[["cooks"]]), range=0, las=2) #plotDensity(log2(counts(dds.norm)+epsilon), col=col.pheno.selected, # xlab="log2(counts)", cex.lab=0.7, panel.first=grid()) #plotDensity(log2(counts(dds.norm, normalized=TRUE)+epsilon), col=col.pheno.selected, # xlab="log2(normalized counts)", cex.lab=0.7, panel.first=grid()) dev.off() # since we Gene-level differential expression analysis with DESeq2, the splicing plays no role in the analysis! # 用nanopore 可以 compare transcript length distribution. 有可能Cellline很长，Extracellular vesicles (EVs)很短! #In the folloing code, we compare different normalization methods, however, the function estimateSizeFactors(dge) doesn't work. I want replace the methods with the method implemented above. Please update the code #uding five norm factors from the following five normalization methods. How to replace the last command " sizeFactors(dds) <- dge2" for the five method? library(ggplot2) library(gridExtra) library(reshape2) library(mixOmics) library(RColorBrewer) library(DESeq) library(edgeR) library(VennDiagram) library(devtools) raw_counts_wn <- raw_counts[rowSums(raw_counts) > 0, ] dim(raw_counts_wn) #--Raw counts-- pseudo_counts <- log2(raw_counts_wn + 1) head(pseudo_counts) df_raw <- melt(pseudo_counts, id = rownames(raw_counts_wn)) names(df_raw)[1:2]<- c("id", "sample") df_raw$method <- rep("Raw counts", nrow(df_raw)) head(df_raw) #--DESeq-- cData = data.frame(row.names=colnames(raw_counts_wn), condition=condition) dge<-DESeqDataSetFromMatrix(countData=raw_counts_wn, colData=cData, design=~condition) dge <- estimateSizeFactors(dge) sizeFactors(dge) # Use your function to get the size factors sizeFactorsDGE <- estimSf(dge) isTRUE(all.equal(sizeFactors(dge), sizeFactorsDGE)) # Test the explained method using sweep to simulate the internal process of dge deseq_normcount <- counts(dge, normalized = TRUE) test_normcount <- sweep(raw_counts_wn, 2, sizeFactors(dge), "/") sum(test_normcount != deseq_normcount) pseudo_deseq <- log2(deseq_normcount + 1) df_deseq <- melt(pseudo_deseq, id = rownames(raw_counts_wn)) names(df_deseq)[1:2]<- c("id", "sample") df_deseq$method <- rep("DESeq (RLE)", nrow(df_raw)) #--edgeR-- dge2 <- DGEList(raw_counts_wn) dge2 dge2$samples #--Total count-- pseudo_TC <- log2(cpm(dge2) + 1) df_TC <- melt(pseudo_TC, id = rownames(raw_counts_wn)) names(df_TC)[1:2] <- c ("id", "sample") df_TC$method <- rep("TC", nrow(df_TC)) ##--RPKM-- #gene_lengths_wn <- gene_lengths[rowSums(raw_counts) > 0] #pseudo_RPKM <- log2(rpkm(dge2, gene.length = gene_lengths_wn) + 1) #df_RPKM <- melt(pseudo_RPKM, id = rownames(raw_counts_wn)) #names(df_RPKM)[1:2] <- c ("id", "sample") #df_RPKM$method <- rep("RPKM", nrow(df_RPKM)) #--Upper quartile-- dge2 <- calcNormFactors(dge2, method = "upperquartile") dge2$samples test_normcount <- sweep(dge2$counts, 2, dge2$samples$lib.size*dge2$samples$norm.factors / 10^6, "/") range(as.vector(test_normcount - cpm(dge2))) pseudo_UQ <- log2(cpm(dge2) + 1) df_UQ <- melt(pseudo_UQ, id = rownames(raw_counts_wn)) names(df_UQ)[1:2] <- c ("id", "sample") df_UQ$method <- rep("UQ", nrow(df_UQ)) #--TMM-- dge2 <- calcNormFactors(dge2, method = "TMM") dge2$samples pseudo_TMM <- log2(cpm(dge2) + 1) df_TMM <- melt(pseudo_TMM, id = rownames(raw_counts_wn)) names(df_TMM)[1:2] <- c ("id", "sample") #MODIFIED! df_TMM$method <- rep("TMM", nrow(df_TMM)) #TMM #--Comparison-- png("normalization.png", width=800, height=600) df_allnorm <- rbind(df_raw, df_deseq, df_TC, df_UQ, df_TMM) df_allnorm$method <- factor(df_allnorm$method, levels = c("Raw counts", "DESeq (RLE)", "TC", "UQ", "TMM")) #df_allnorm <- rbind(df_raw, df_TMM) #df_allnorm$method <- factor(df_allnorm$method, levels = c("Raw counts", "DESeq (RLE)")) p <- ggplot(data=df_allnorm, aes(x=sample, y=value, fill=method)) p <- p + geom_boxplot() p <- p + theme_bw() p <- p + ggtitle("Boxplots of normalized pseudo counts\n for all samples by normalization methods") p <- p + facet_grid(. ~ method) p <- p + ylab(expression(log[2] ~ (normalized ~ count + 1))) + xlab("") p <- p + theme(title = element_text(size=10), axis.text.x = element_blank(), axis.ticks.x = element_blank()) print(p) dev.off() #Assign a normalization size factor for downstream analysis from "Raw counts", "DESeq (RLE)", "TC", "UQ", "TMM") #You will need to repeat the DESeq2 analysis steps (like DESeq()) for each set of size factors to see how the results change with different normalization methods. #Note: This approach allows you to apply normalization factors from various methods to DESeq2. But keep in mind that each normalization method was developed with a specific intent and assumptions. Combining normalization methods from one package (like edgeR's TMM) with differential analysis from another package (like DESeq2) might not always be theoretically sound. Always interpret results with caution and preferably in consultation with domain experts. #DESeq (RLE) normalization: sizeFactors(dds) <- sizeFactors(dge) #control_r1 control_r2 HSV.d2_r1 HSV.d2_r2 HSV.d4_r1 HSV.d4_r2 HSV.d6_r1 #0.01996676 0.03367626 0.33493034 0.65395381 4.96825596 4.14671012 4.07416461 #HSV.d6_r2 HSV.d8_r1 HSV.d8_r2 #4.43844964 5.96130192 3.75356239 #Total Count normalization: This is actually just scaling by the total library size. In DESeq2, this would be equivalent to: sizeFactors(dds) <- colSums(raw_counts_wn) #control_r1 control_r2 HSV.d2_r1 HSV.d2_r2 HSV.d4_r1 HSV.d4_r2 HSV.d6_r1 # 14594 18992 182152 374750 3404540 2521691 2646741 #HSV.d6_r2 HSV.d8_r1 HSV.d8_r2 # 3207708 3873254 2430625 identical(colnames(dds), rownames(dge2$samples)) #Upper Quartile (UQ) normalization: uq_factors <- dge2$samples$norm.factors * dge2$samples$lib.size / 10^6 names(uq_factors) <- rownames(dge2$samples) sizeFactors(dds) <- uq_factors #control_r1 control_r2 HSV.d2_r1 HSV.d2_r2 HSV.d4_r1 HSV.d4_r2 HSV.d6_r1 #0.01578156 0.02060450 0.19559938 0.39108853 3.27858657 2.38221969 2.61381808 # HSV.d6_r2 HSV.d8_r1 HSV.d8_r2 #3.03441572 3.65683812 2.30405047 #TMM normalization: tmm_factors <- dge2$samples$norm.factors names(tmm_factors) <- rownames(dge2$samples) sizeFactors(dds) <- tmm_factors #> tmm_factors #control_r1 control_r2 HSV.d2_r1 HSV.d2_r2 HSV.d4_r1 HSV.d4_r2 HSV.d6_r1 #1.0813729 1.0849043 1.0738250 1.0435985 0.9630043 0.9446914 0.9875610 #HSV.d6_r2 HSV.d8_r1 HSV.d8_r2 #0.9459763 0.9441256 0.9479251 #Given the raw counts, the control_r1 and control_r2 samples seem to have a much lower sequencing depth (total read count) than the other samples. Therefore, when normalization methods such as TMM (Trimmed Mean of M-values) or UQ (Upper Quartile normalization) are applied, the normalization factors for these control samples will be relatively high, boosting the normalized counts. #To better understand the situation: # - Check the sequencing depth of each sample. Sum the raw counts across all genes for each sample to get the total read count. If the control samples have a substantially lower total read count, then normalization methods will try to adjust for this discrepancy. # - Consider the normalization method: Different normalization methods might provide slightly different results. TMM normalization, for example, tries to adjust for the compositional differences between samples. It's common to observe larger normalization factors for samples with a lower total read count. # - Visualize the data: MA plots, box plots, or density plots of the raw and normalized counts can help understand the distribution of counts and the effect of normalization. # ---- Adapt the following code to print the normalized_counts using tmm_factors ---- # RLE size factors rle_factors <- sizeFactors(dge) # Print side by side data.frame(RLE = rle_factors, TMM = tmm_factors) sizeFactors(dds) <- rle_factors normalized_counts_rle <- counts(dds, normalized=TRUE) sizeFactors(dds) <- tmm_factors #is it possible in the following command, the calculate de novo normalization factors and ignore the given tmm normalization factors? normalized_counts_tmm <- counts(dds, normalized=TRUE) difference_matrix <- normalized_counts_tmm - normalized_counts_rle all(difference_matrix == 0) -
select the differentially expressed genes
#https://galaxyproject.eu/posts/2020/08/22/three-steps-to-galaxify-your-tool/ #https://www.biostars.org/p/282295/ #https://www.biostars.org/p/335751/ #> condition # [1] control control HSV.d2 HSV.d2 HSV.d4 HSV.d4 HSV.d6 HSV.d6 HSV.d8 HSV.d8 #Levels: control HSV.d2 HSV.d4 HSV.d6 HSV.d8 #CONSOLE: mkdir /home/jhuang/DATA/Data_Manja_RNAseq_Organoids_Virus/results_chrHsv1_downstream/star_salmon/degenes setwd("/home/jhuang/DATA/Data_Manja_RNAseq_Organoids_Merged") #---- relevel to control ---- dds$condition <- relevel(dds$condition, "control") dds = DESeq(dds, betaPrior=FALSE) resultsNames(dds) clist <- c("HSV.d2_vs_control","HSV.d4_vs_control", "HSV.d6_vs_control", "HSV.d8_vs_control") dds$condition <- relevel(dds$condition, "HSV.d2") dds = DESeq(dds, betaPrior=FALSE) resultsNames(dds) clist <- c("HSV.d4_vs_HSV.d2", "HSV.d6_vs_HSV.d2", "HSV.d8_vs_HSV.d2") dds$condition <- relevel(dds$condition, "HSV.d4") dds = DESeq(dds, betaPrior=FALSE) resultsNames(dds) clist <- c("HSV.d6_vs_HSV.d4", "HSV.d8_vs_HSV.d4") dds$condition <- relevel(dds$condition, "HSV.d6") dds = DESeq(dds, betaPrior=FALSE) resultsNames(dds) clist <- c("HSV.d8_vs_HSV.d6") library(biomaRt) listEnsembl() listMarts() #--> total 69, 27 GRCh38.p7 and 39 GRCm38.p4 ensembl <- useEnsembl(biomart = "ensembl", dataset = "hsapiens_gene_ensembl", version="104") datasets <- listDatasets(ensembl) # -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[1] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[1] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # Define the original and compressed ranges original_range <- c(20, 40) compressed_range <- c(20.0, 24.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 20, by=5) y_breaks_compressed <- c(20.0, 24.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 20, by=5) y_labels_compressed <- c(20, 40) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off() # -- HSV.d8_vs_control.png -- # Define the original and compressed ranges original_range <- c(80, 115) compressed_range <- c(80.0, 90.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 80, by=10) y_breaks_compressed <- c(80.0, 90.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 80, by=10) y_labels_compressed <- c(80, 115) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # -- HSV.d6_vs_control.png -- # Define the original and compressed ranges original_range <- c(80, 115) compressed_range <- c(80.0, 90.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 80, by=10) y_breaks_compressed <- c(80.0, 90.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 80, by=10) y_labels_compressed <- c(80, 115) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # -- HSV.d4_vs_control.png -- # Define the original and compressed ranges original_range <- c(80, 100) compressed_range <- c(80.0, 90.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 80, by=10) y_breaks_compressed <- c(80.0, 90.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 80, by=10) y_labels_compressed <- c(80, 100) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # -- HSV.d2_vs_control.png -- # Define the original and compressed ranges original_range <- c(20, 40) compressed_range <- c(20.0, 24.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 20, by=5) y_breaks_compressed <- c(20.0, 24.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 20, by=5) y_labels_compressed <- c(20, 40) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) #under DIR degenes under KONSOLE for i in HSV.d2_vs_control HSV.d4_vs_control HSV.d6_vs_control HSV.d8_vs_control HSV.d4_vs_HSV.d2 HSV.d6_vs_HSV.d2 HSV.d8_vs_HSV.d2 HSV.d6_vs_HSV.d4 HSV.d8_vs_HSV.d4 HSV.d8_vs_HSV.d6; do echo "~/Tools/csv2xls-0.4/csv_to_xls.py ${i}-all_annotated.txt ${i}-up_annotated.txt ${i}-down_annotated.txt -d$',' -o ${i}.xls;"; done -
clustering the genes and draw heatmap
for i in HSV.d2_vs_control HSV.d4_vs_control HSV.d6_vs_control HSV.d8_vs_control HSV.d4_vs_HSV.d2 HSV.d6_vs_HSV.d2 HSV.d8_vs_HSV.d2 HSV.d6_vs_HSV.d4 HSV.d8_vs_HSV.d4 HSV.d8_vs_HSV.d6; do echo "cut -d',' -f1-1 ${i}-up_annotated.txt > ${i}-up.id"; echo "cut -d',' -f1-1 ${i}-down_annotated.txt > ${i}-down.id"; done cat *.id | sort -u > ids #add Gene_Id in the first line, delete the "" GOI <- read.csv("ids")$Gene_Id #4647 RNASeq.NoCellLine <- assay(rld) #install.packages("gplots") library("gplots") #clustering methods: "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC). pearson or spearman datamat = RNASeq.NoCellLine[GOI, ] #datamat = RNASeq.NoCellLine write.csv(as.data.frame(datamat), file ="gene_expressions.txt") constant_rows <- apply(datamat, 1, function(row) var(row) == 0) if(any(constant_rows)) { cat("Removing", sum(constant_rows), "constant rows.\n") datamat <- datamat[!constant_rows, ] } hr <- hclust(as.dist(1-cor(t(datamat), method="pearson")), method="complete") hc <- hclust(as.dist(1-cor(datamat, method="spearman")), method="complete") mycl = cutree(hr, h=max(hr$height)/1.05) mycol = c("YELLOW", "BLUE", "ORANGE", "MAGENTA", "CYAN", "RED", "GREEN", "MAROON", "LIGHTBLUE", "PINK", "MAGENTA", "LIGHTCYAN", "LIGHTRED", "LIGHTGREEN"); mycol = mycol[as.vector(mycl)] #png("DEGs_heatmap.png", width=900, height=800) #cex.lab=10, labRow="", png("DEGs_heatmap.png", width=800, height=1000) heatmap.2(as.matrix(datamat),Rowv=as.dendrogram(hr),Colv = NA, dendrogram = 'row',labRow="", scale='row',trace='none',col=bluered(75), cexCol=1.8, RowSideColors = mycol, margins=c(10,2), cexRow=1.5, srtCol=30, lhei = c(1, 8), lwid=c(2, 8)) #rownames(datamat) #heatmap.2(datamat, Rowv=as.dendrogram(hr), col=bluered(75), scale="row", RowSideColors=mycol, trace="none", margin=c(5,5), sepwidth=c(0,0), dendrogram = 'row', Colv = 'false', density.info='none', labRow="", srtCol=30, lhei=c(0.1,2)) dev.off() #### cluster members ##### write.csv(names(subset(mycl, mycl == '1')),file='cluster1_YELLOW.txt') write.csv(names(subset(mycl, mycl == '2')),file='cluster2_DARKBLUE.txt') write.csv(names(subset(mycl, mycl == '3')),file='cluster3_DARKORANGE.txt') write.csv(names(subset(mycl, mycl == '4')),file='cluster4.txt') #~/Tools/csv2xls-0.4/csv_to_xls.py cluster*.txt -d',' -o DEGs_heatmap_cluster_members.xls ~/Tools/csv2xls-0.4/csv_to_xls.py \ significant_gene_expressions.txt \ -d',' -o DEGs_heatmap_expression_data.xls; #### cluster members (advanced) ##### subset_1<-names(subset(mycl, mycl == '1')) data <- as.data.frame(datamat[rownames(datamat) %in% subset_1, ]) #2575 subset_2<-names(subset(mycl, mycl == '2')) data <- as.data.frame(datamat[rownames(datamat) %in% subset_2, ]) #1855 subset_3<-names(subset(mycl, mycl == '3')) data <- as.data.frame(datamat[rownames(datamat) %in% subset_3, ]) #217 subset_4<-names(subset(mycl, mycl == '4')) data <- as.data.frame(datamat[rownames(datamat) %in% subset_4, ]) # subset_5<-names(subset(mycl, mycl == '5')) data <- as.data.frame(datamat[rownames(datamat) %in% subset_5, ]) # # Initialize an empty data frame for the annotated data annotated_data <- data.frame() # Determine total number of genes total_genes <- length(rownames(data)) # Loop through each gene to annotate for (i in 1:total_genes) { gene <- rownames(data)[i] result <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = gene, mart = ensembl) # If multiple rows are returned, take the first one if (nrow(result) > 1) { result <- result[1, ] } # Check if the result is empty if (nrow(result) == 0) { result <- data.frame(ensembl_gene_id = gene, external_gene_name = NA, gene_biotype = NA, entrezgene_id = NA, chromosome_name = NA, start_position = NA, end_position = NA, strand = NA, description = NA) } # Transpose expression values expression_values <- t(data.frame(t(data[gene, ]))) colnames(expression_values) <- colnames(data) # Combine gene information and expression data combined_result <- cbind(result, expression_values) # Append to the final dataframe annotated_data <- rbind(annotated_data, combined_result) # Print progress every 100 genes if (i %% 100 == 0) { cat(sprintf("Processed gene %d out of %d\n", i, total_genes)) } } # Save the annotated data to a new CSV file #write.csv(annotated_data, "cluster1_YELLOW.csv", row.names=FALSE) write.csv(annotated_data, "cluster2_DARKBLUE.csv", row.names=FALSE) write.csv(annotated_data, "cluster3_DARKORANGE.csv", row.names=FALSE) write.csv(annotated_data, "cluster4_DARKMAGENTA.csv", row.names=FALSE) write.csv(annotated_data, "cluster5_DARKCYAN.csv", row.names=FALSE) #~/Tools/csv2xls-0.4/csv_to_xls.py cluster*.csv -d',' -o DEGs_heatmap_clusters.xls -
code of differential gene analysis for clist <- c("HSV.d2_vs_control","HSV.d4_vs_control", "HSV.d6_vs_control", "HSV.d8_vs_control")
# -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[1] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # -- HSV.d2_vs_control.png -- # Define the original and compressed ranges original_range <- c(20, 40) compressed_range <- c(20.0, 24.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 20, by=5) y_breaks_compressed <- c(20.0, 24.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 20, by=5) y_labels_compressed <- c(20, 40) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off() # -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # -- HSV.d4_vs_control.png -- # Define the original and compressed ranges original_range <- c(80, 100) compressed_range <- c(80.0, 90.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 80, by=10) y_breaks_compressed <- c(80.0, 90.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 80, by=10) y_labels_compressed <- c(80, 100) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off() # -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[3] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # -- HSV.d6_vs_control.png -- # Define the original and compressed ranges original_range <- c(80, 115) compressed_range <- c(80.0, 90.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 80, by=10) y_breaks_compressed <- c(80.0, 90.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 80, by=10) y_labels_compressed <- c(80, 115) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off() # -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[4] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # -- HSV.d8_vs_control.png -- # Define the original and compressed ranges original_range <- c(80, 115) compressed_range <- c(80.0, 90.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 80, by=10) y_breaks_compressed <- c(80.0, 90.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 80, by=10) y_labels_compressed <- c(80, 115) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off() -
code of differential gene analysis for clist <- c("HSV.d4_vs_HSV.d2", "HSV.d6_vs_HSV.d2", "HSV.d8_vs_HSV.d2")
# -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[1] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # Define the original and compressed ranges original_range <- c(15, 20) compressed_range <- c(15.0, 20.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 15, by=5) y_breaks_compressed <- c(15.0, 20.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 15, by=5) y_labels_compressed <- c(15, 20) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + #annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + #annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off() # -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # Define the original and compressed ranges original_range <- c(40, 50) compressed_range <- c(40.0, 50.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 40, by=10) y_breaks_compressed <- c(40.0, 50.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 40, by=10) y_labels_compressed <- c(40, 50) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + #annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + #annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off() # -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[3] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # Define the original and compressed ranges original_range <- c(40, 50) compressed_range <- c(40.0, 50.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 40, by=10) y_breaks_compressed <- c(40.0, 50.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 40, by=10) y_labels_compressed <- c(40, 50) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + #annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + #annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off() -
code of differential gene analysis for clist <- c("HSV.d6_vs_HSV.d4", "HSV.d8_vs_HSV.d4")
# -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[1] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # Define the original and compressed ranges original_range <- c(15, 20) compressed_range <- c(15.0, 20.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 15, by=5) y_breaks_compressed <- c(15.0, 20.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 15, by=5) y_labels_compressed <- c(15, 20) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + #annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + #annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off() # -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # Define the original and compressed ranges original_range <- c(40, 50) compressed_range <- c(40.0, 50.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 40, by=5) y_breaks_compressed <- c(40.0, 50.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 40, by=5) y_labels_compressed <- c(40, 50) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + #annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + #annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off() -
code of differential gene analysis for clist <- c("HSV.d8_vs_HSV.d6")
# -- 1. export res_df containing both human and virus genes -- #for (i in clist) { i<-clist[1] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] res_df <- as.data.frame(res) write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-")) up <- subset(res_df, padj<=0.05 & log2FoldChange>=2) down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-")) #} # -- 2. annatete human genes 'geness_res', note that the virus genes in that is ignored in the process since they are not in the database -- #for (i in clist) { #i<-clist[2] contrast = paste("condition", i, sep="_") res = results(dds, name=contrast) res <- res[!is.na(res$log2FoldChange),] #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("ENTREZID","EXONID","GENEBIOTYPE","GENEID","GENENAME","PROTEINDOMAINSOURCE","PROTEINID","SEQNAME","SEQSTRAND","SYMBOL","TXBIOTYPE","TXID","TXNAME","UNIPROTID")) #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "ENTREZID", "SYMBOL", "GENENAME","GENEBIOTYPE","TXBIOTYPE","SEQSTRAND","UNIPROTID")) # In the ENSEMBL-database, GENEID is ENSEMBL-ID. #geness <- AnnotationDbi::select(edb86, keys = rownames(res), keytype = "GENEID", columns = c("GENEID", "SYMBOL", "GENEBIOTYPE")) # "ENTREZID", "TXID","TXBIOTYPE","TXSEQSTART","TXSEQEND" #geness <- geness[!duplicated(geness$GENEID), ] #using getBM replacing AnnotationDbi::select #filters = 'ensembl_gene_id' means the records should always have a valid ensembl_gene_ids. geness <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'), filters = 'ensembl_gene_id', values = rownames(res), mart = ensembl) geness_uniq <- distinct(geness, ensembl_gene_id, .keep_all= TRUE) #merge by column by common colunmn name, in the case "GENEID" res$ENSEMBL = rownames(res) identical(rownames(res), rownames(geness_uniq)) res_df <- as.data.frame(res) geness_res <- merge(geness_uniq, res_df, by.x="ensembl_gene_id", by.y="ENSEMBL") dim(geness_res) rownames(geness_res) <- geness_res$ensembl_gene_id geness_res$ensembl_gene_id <- NULL #} # -- 3. prepare annatete virus genes -- virus_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") # replace with the actual list of column names virus_rows <- res_df[rownames(res_df) %in% virus_genes, ] virus_rows$external_gene_name <- rownames(virus_rows) virus_rows$chromosome_name <- "chrHsv1" # Define default values based on data type default_values <- list( character = NULL, numeric = 0, integer = 0L, logical = FALSE ) # Ensure that virus_rows has the same columns as geness_res for (col in colnames(geness_res)) { if (!col %in% colnames(virus_rows)) { data_type <- class(geness_res[[col]])[1] default_value <- default_values[[data_type]] virus_rows[[col]] <- rep(default_value, nrow(virus_rows)) } } missing_cols <- setdiff(colnames(geness_res), colnames(virus_rows)) for (col in missing_cols) { virus_rows[[col]] <- NA # Or another default value as appropriate } # Reorder columns in virus_rows to match the order in geness_res virus_rows <- virus_rows[, colnames(geness_res), drop = FALSE] # -- 4. merge them together -- #for (i in clist) { merged_df <- rbind(geness_res, virus_rows) merged_df_sorted <- as.data.frame(merged_df[order(merged_df$padj),]) write.csv(merged_df_sorted, file = paste(i, "all_annotated.txt", sep="-")) up <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange>=2) down <- subset(merged_df_sorted, padj<=0.05 & log2FoldChange<=-2) write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up_annotated.txt", sep="-")) write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down_annotated.txt", sep="-")) #} # -- 5. draw graphics -- #geness_res <- read.csv(file = "HSV.d2_vs_control-all_annotated.txt", sep=",", row.names=1) geness_res <- merged_df_sorted # Color setting geness_res$Color <- ifelse(geness_res$padj > 0.05 | abs(geness_res$log2FoldChange) < 2, "gray", ifelse(geness_res$log2FoldChange > 0, "red", "blue")) # Predefined genes colored in green predefined_genes <- c("AL", "IRL1", "IRL2", "IRL3", "IRS1", "LAT", "ORF-O/P", "Ori-S", "pri-miRNA", "TRL1", "TRL2", "TRL3", "TRS1", "UL1", "UL10", "UL11", "UL12", "UL13", "UL14", "UL15", "UL16", "UL17", "UL18", "UL19", "UL2", "UL20", "UL21", "UL22", "UL23", "UL24", "UL25", "UL26", "UL27", "UL28", "UL29", "UL3", "UL30", "UL31", "UL32", "UL33", "UL34", "UL35", "UL36", "UL37", "UL38", "UL39", "UL4", "UL40", "UL41", "UL42", "UL43", "UL44", "UL45", "UL46", "UL47", "UL48", "UL49", "UL5", "UL50", "UL51", "UL52", "UL53", "UL54", "UL55", "UL56", "UL6", "UL7", "UL8", "UL9", "US1", "US10", "US11", "US12", "US2", "US3", "US4", "US5", "US6", "US7", "US8", "US9") geness_res$Color[geness_res$external_gene_name %in% predefined_genes] <- "green" geness_res$invert_Padj <- (-log10(geness_res$padj)) * sign(geness_res$log2FoldChange) top_g <- unique(c(geness_res[order(geness_res$invert_Padj, decreasing = TRUE), 'external_gene_name'][1:100], geness_res[order(geness_res$invert_Padj, decreasing = FALSE), 'external_gene_name'][1:100])) # Define the original and compressed ranges original_range <- c(40, 50) compressed_range <- c(40.0, 50.0) # Calculate breaks for the y-axis y_breaks_below <- seq(0, 40, by=5) y_breaks_compressed <- c(40.0, 50.0) y_breaks_above <- c() y_breaks <- c(y_breaks_below, y_breaks_compressed, y_breaks_above) y_labels_below <- seq(0, 40, by=5) y_labels_compressed <- c(40, 50) y_labels_above <- c() y_labels <- c(y_labels_below, y_labels_compressed, y_labels_above) # Adjust the p-values based on the ranges geness_res$adjusted_pvalue <- with(geness_res, ifelse(-log10(padj) > original_range[1] & -log10(padj) <= original_range[2], ((-log10(padj) - original_range[1]) / (original_range[2] - original_range[1])) * (compressed_range[2] - compressed_range[1]) + compressed_range[1], ifelse(-log10(padj) > original_range[2], -log10(padj) - (original_range[2] - original_range[1]) + (compressed_range[2] - compressed_range[1]), -log10(padj)))) # Create the plot png(paste(i, "png", sep="."), width=1000, height=1000) ggplot(geness_res, aes(x = log2FoldChange, y = adjusted_pvalue, color = Color, label = external_gene_name)) + geom_vline(xintercept = c(2, -2), lty = "dashed", size = 1.5) + geom_hline(yintercept = -log10(0.05), lty = "dashed", size = 1.5) + geom_point(size = 3) + labs(x = "log2(Fold change)", y = "-log10(P-adj)", color = "Significance") + scale_color_identity() + geom_text_repel(data = subset(geness_res, external_gene_name %in% top_g & padj < 0.05 & (abs(log2FoldChange) >= 2)), size = 7, point.padding = 0.15, color = "black", min.segment.length = .1, box.padding = .2, lwd = 2) + theme_bw(base_size = 24) + theme(legend.position = "bottom") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[1], ymax = compressed_range[1], linetype = "dashed", color = "grey") + #annotate("rect", xmin = -Inf, xmax = Inf, ymin = compressed_range[2], ymax = compressed_range[2], linetype = "dashed", color = "grey") + #annotate("text", x = -Inf, y = compressed_range[1], label = "/", hjust = 0, size = 10) + #annotate("text", x = -Inf, y = compressed_range[2], label = "/", hjust = 0, size = 10) + scale_y_continuous(breaks = sort(y_breaks), labels = sort(y_labels)) dev.off()