3C_whippet_pseudobulk
2025-03-20
Calculating the TPM values for the P1 & P2 promoters and then calculating the knockdown efficiency for off-target and on-target targets, as well as generating visualisations of the data.
The provided R code serves as the core validation engine, quantifying transcript-level changes to prove that CRISPR interference (CRISPRi) can selectively repress one promoter without affecting others in the same gene.
Isoform Quantification: It utilises Whippet outputs to quantify expression in Transcripts Per Million (TPM), allowing for a high-resolution view of which specific promoter is active in which cell.
On-Target vs. Off-Target Precision: It calculates the percentage knockdown for transcripts initiated by the targeted promoter (On-target) versus transcripts from untargeted promoters of the same gene (Off-target)
Quality Filtering: The script establishes strict criteria for “Successful KD” (e.g., P2 knockdown < 40%) and ensures that background expression in non-targeting control (NTC) populations is sufficient (NTC TPM ≥ 0.5) to make reliable conclusions.
Load libraries and define functions.
Define the function and the files.
rel_dir <- '../../'
dir <- '../../files/whippet_dedup'
files <- list.files(dir, pattern = ".gz", full.names = TRUE)
table_name <- file.path(dir, "all_data.txt")
# Print the list of files
print(files[1:3])## [1] "../../files/whippet_dedup/ADAR_AP.deduplicated.isoform.tpm.gz"
## [2] "../../files/whippet_dedup/ADAR_MP.deduplicated.isoform.tpm.gz"
## [3] "../../files/whippet_dedup/AHCYL1_AP.deduplicated.isoform.tpm.gz"Read through entire whippet output and create gene table if it does not exist. Here the off-target and on-target P1 and P2 knockdown percentages are calculated. Other information such as the transcripts assigned to the P1 anad P2 promoters, the distance from the TSS and the negative control TPMs are also calculated. Finally, boolean values for if the KD is significant, if there is a significant NTC count (> 0.5 TPM in negative control population), if KD is looser (less than the 85th percentile) or strict (both P1 and P2 < 0% KD) are also calculated.
#if gene_table.csv exists then skip
if (file.exists(paste0(dir,"/gene_table.csv"))){
print("gene_table.csv exists")
#load the gene_table.csv
gene_table <- read.csv(paste0(dir,"/gene_table.csv"))
} else {
print("gene_table.csv does not exist")
# Read or create the table
if (file.exists(table_name)) {
all_data <- read.table(table_name, header = TRUE, sep = "\t", row.names = NULL)
} else {
#map the fiunction process to every element in list of files and merge the otuputs to simplest df use applyr
all_data <- lapply(files, process_file)
#flatten and merge on isforom column
all_data <- Reduce(function(x, y) merge(x, y, by = "Isoform", all = TRUE), all_data)
#instead of cbind use merge on Isoform
write.table(all_data, table_name, sep = "\t", row.names = FALSE, quote = FALSE)
}
#tx2gene has list of transcripts and their associated genes
tx2gene_loc <- "files/kallisto/transcriptome/tx2gene.txt"
tx2gene <- read_delim(tx2gene_loc, delim = "\t", col_names = c("transcript_id", "gene_id","gene_name","transcript_name","chromosome","start","end","strand"))
print(sum((tx2gene$start - tx2gene$end) > 0))
tx2gene$transcript_index <- 1:nrow(tx2gene)
tx2gene$tss <- ifelse(tx2gene$strand == "+", tx2gene$start, tx2gene$end)
#read in the cells that are assigned, place into metadata find
#read in the bed file which conrtains p1 and p2 locations
bed_loc <- "cellranger_input/candidate_AP_MP.bed"
bed <- read_delim(bed_loc, delim = "\t", col_names = c("chromosome", "start", "end", "name", "strand"))
print(sum((bed$start - bed$end) > 0))
bed$tss <- ifelse(bed$strand == "+", bed$start, bed$end)
#do an example for esr1 mp and ap subpopulations
#split _ and get last
bed$promoter_type <- str_split(bed$name, "_") %>% map_chr(~last(.x))
#get the second to last element
bed$gene_name <- str_split(bed$name, "_") %>% map_chr(~.x[length(.x)-1])
bed$promoter_gene <- paste0(bed$gene_name,"_",bed$promoter_type)
#check the number of gene name prior to merge
print(length(unique(bed$gene_name)))
print(length(unique(bed$promoter_gene)))
#merge the dataframe on gene name
#merge the bed and tx2gene on gene name
bed <- merge(bed, tx2gene, by = "gene_name", suffixes = c("_target","_assigned"), how="inner")
print(length(unique(bed$gene_name)))
print(length(unique(bed$promoter_gene)))
#filter for for rows that only have
#calculate column that has the distance from the tss
bed$distance <- abs(bed$tss_target - bed$tss_assigned)
#filter for rows that have a distance of less than 1000
bed <- bed %>% group_by(transcript_id) %>% mutate(min_distance = distance == min(distance))
print(length(unique(bed$transcript_id))==length(bed$transcript_id[bed$min_distance == TRUE]))
ggplot(bed) +
labs(x = "Distance from TSS", y = "Count") + #add denistty plot
geom_density(aes(x=distance, color=min_distance), alpha = 0.5)
#then find the numbers that assign to those transcripts
#get the number from the matrix number for those
bed <- bed %>% filter(min_distance == TRUE)
#only keep genes with both P1 and P2
all_data_names <- stringr::str_split_i(colnames(all_data), "Tpm_", 2)
#use these transcripts
bed <- bed %>% group_by(gene_name) %>% filter(promoter_gene %in% unique(all_data_names)) %>% filter(n_distinct(promoter_type) == 2)
bed <- bed %>% filter(distance < 1000)
#split transcript ID by . and get the first element name column ISoform
bed$Isoform <- str_split_fixed(bed$transcript_id, "\\.", 2)[,1]
#merge the bed and all_data on Isoform
all_data <- merge(all_data, bed, by = "Isoform")
head(all_data)
#create empty table gene_table with columsn for both mp and ap
gene_table_ontarget <- data.frame(gene = character(), percent_knockdown_mp = numeric(), distance_mp = numeric(), Isoform_mp = character(), transcript_name_mp = character(), neg_ctrl_tpm_mp= character(), percent_knockdown_ap = numeric(), distance_ap = numeric(), Isoform_ap = character(), transcript_name_ap = character(), neg_ctrl_tpm_ap= character(), stringsAsFactors = FALSE)
gene_table_offtarget <- data.frame(gene = character(), percent_knockdown_mp = numeric(), distance_mp = numeric(), Isoform_mp = character(), transcript_name_mp = character(), neg_ctrl_tpm_mp= character(),percent_knockdown_ap = numeric(), distance_ap = numeric(), Isoform_ap = character(), transcript_name_ap = character(), neg_ctrl_tpm_ap= character(),stringsAsFactors = FALSE)
gene_table_total <- data.frame(gene = character(), percent_knockdown_mp_total = numeric(), percent_knockdown_ap = numeric(), stringsAsFactors = FALSE)
#loop through every gene
#unique gene names
for (gene in unique(all_data$gene_name)){
print(gene)
gene_prom_mp <- paste0(gene, "_MP")
gene_prom_mp_tpm <- paste0("Tpm_", gene_prom_mp)
gene_prom_ap <- paste0(gene, "_AP")
gene_prom_ap_tpm <- paste0("Tpm_", gene_prom_ap)
negative_column <- "Tpm_Negative"
#subset the data for the gene
#firrst for main promoter
gene_data_mp <- all_data %>% filter(promoter_gene == gene_prom_mp)
#then for the alternative promoter
gene_data_ap <- all_data %>% filter(promoter_gene == gene_prom_ap)
#check if all the sum(gene_data_mp[[gene_prom_mp_tpm]]) or sum(gene_data_mp[[gene_prom_ap_tpm]]) for every column is greater than the minimum filter tpm
# if (( > filter) & (sum(gene_data_ap[[negative_column]]) > filter) ){
# print("skipped")
# next
# }
print(sum(gene_data_mp[[negative_column]]))
print(sum(gene_data_ap[[negative_column]]))
#calculate the % knockdown between negative control and ESR1_MP population by summing the Tpm and then dividing by the negative control
#fill NA wuth 0 in gene_data_ap[[gene_prom_ap_tpm]] and gene_data_mp[[gene_prom_mp_tpm]] and negative control and gene_data_mp[[gene_prom_ap_tpm]]
total_mp <- sum(gene_data_mp[[gene_prom_mp_tpm]]) + sum(gene_data_mp[[gene_prom_ap_tpm]])
total_ap <- sum(gene_data_ap[[gene_prom_mp_tpm]]) + sum(gene_data_ap[[gene_prom_ap_tpm]])
total_percentage_knockdown_mp <- (total_mp - sum(gene_data_mp[[negative_column]]))/sum(gene_data_mp[[negative_column]])*100
total_percentage_knockdown_ap <- (total_ap - sum(gene_data_ap[[negative_column]]))/sum(gene_data_ap[[negative_column]])*100
gene_data_mp_percent_knockdown_mp <- ( sum(gene_data_mp[[gene_prom_mp_tpm]])-sum(gene_data_mp[[negative_column]]))/sum(gene_data_mp[[negative_column]])*100
gene_data_ap_percent_knockdown_ap <- ( sum(gene_data_ap[[gene_prom_ap_tpm]])-sum(gene_data_ap[[negative_column]]))/sum(gene_data_ap[[negative_column]])*100
gene_data_mp_percent_knockdown_ap <- ( sum(gene_data_mp[[gene_prom_ap_tpm]])-sum(gene_data_mp[[negative_column]]))/sum(gene_data_mp[[negative_column]])*100
gene_data_ap_percent_knockdown_mp <- ( sum(gene_data_ap[[gene_prom_mp_tpm]])-sum(gene_data_ap[[negative_column]]))/sum(gene_data_ap[[negative_column]])*100
# whenever sum(gene_data_mp[[gene_prom_mp_tpm]]) == sum(gene_data_mp[[negative_column]])) make the value equal 0
gene_data_mp_percent_knockdown_mp[sum(gene_data_mp[[gene_prom_mp_tpm]]) == sum(gene_data_mp[[negative_column]])] <- 0
gene_data_ap_percent_knockdown_ap[sum(gene_data_ap[[gene_prom_ap_tpm]]) == sum(gene_data_ap[[negative_column]])] <- 0
gene_data_mp_percent_knockdown_ap[sum(gene_data_mp[[gene_prom_ap_tpm]]) == sum(gene_data_mp[[negative_column]])] <- 0
gene_data_ap_percent_knockdown_mp[sum(gene_data_ap[[gene_prom_mp_tpm]]) == sum(gene_data_ap[[negative_column]])] <- 0
gene_list_kd <- list(gene,total_percentage_knockdown_mp,total_percentage_knockdown_ap)
gene_list_ontarget <- list(gene, sum(gene_data_mp[[gene_prom_mp_tpm]]),
gene_data_mp_percent_knockdown_mp,median(gene_data_mp$distance),str_flatten_comma(gene_data_mp$Isoform),str_flatten_comma(gene_data_mp$transcript_name),sum(gene_data_mp[[negative_column]]),
sum(gene_data_ap[[gene_prom_ap_tpm]]),
gene_data_ap_percent_knockdown_ap,median(gene_data_ap$distance),str_flatten_comma(gene_data_ap$Isoform),str_flatten_comma(gene_data_ap$transcript_name),sum(gene_data_ap[[negative_column]]))
gene_list_offtarget <- list(gene, sum(gene_data_mp[[gene_prom_ap_tpm]]),
gene_data_mp_percent_knockdown_ap,median(gene_data_mp$distance),str_flatten_comma(gene_data_mp$Isoform),str_flatten_comma(gene_data_mp$transcript_name),sum(gene_data_mp[[negative_column]]),
sum(gene_data_ap[[gene_prom_mp_tpm]]),
gene_data_ap_percent_knockdown_mp,median(gene_data_ap$distance),str_flatten_comma(gene_data_ap$Isoform),str_flatten_comma(gene_data_ap$transcript_name),sum(gene_data_ap[[negative_column]]))
#add list to table keep column names
print(gene_prom_ap_tpm)
gene_table_ontarget <- rbind(gene_table_ontarget, gene_list_ontarget)
gene_table_offtarget <- rbind(gene_table_offtarget, gene_list_offtarget)
gene_table_total <- rbind(gene_table_total, gene_list_kd)
}
#change gene_table names
column_names <- c("gene", "sum_MP","percent_knockdown_MP", "distance_MP", "Isoform_MP", "transcript_name_MP", "neg_ctrl_tpm_MP",
"sum_AP","percent_knockdown_AP", "distance_AP", "Isoform_AP", "transcript_name_AP","neg_ctrl_tpm_AP")
colnames(gene_table_ontarget) <-column_names
colnames(gene_table_offtarget) <- column_names
colnames(gene_table_total) <- c("gene", "percent_knockdown_MP_total", "percent_knockdown_AP_total")
gene_table <- merge(gene_table_ontarget, gene_table_offtarget, by = "gene", suffixes = c("_ontarget", "_offtarget"))
gene_table <- merge(gene_table, gene_table_total, by = "gene")
head(gene_table)
#plot the data
gene_table <- as.data.frame(gene_table)
dim(gene_table)
gene_table$successfulKD <- ifelse(gene_table$percent_knockdown_AP_ontarget < 40,TRUE,FALSE)
quantile_mp_on <- quantile(gene_table$percent_knockdown_MP_ontarget, c(0.1,0.9))
quantile_ap_on <- quantile(gene_table$percent_knockdown_AP_ontarget, c(0.1,0.85))
gene_table$sig_ntc_count <- ifelse(gene_table$neg_ctrl_tpm_MP_ontarget >= 0.5 & gene_table$neg_ctrl_tpm_AP_ontarget >= 0.5,TRUE,FALSE)
gene_table$looser_KD <- ifelse(( (gene_table$percent_knockdown_AP_ontarget < quantile_ap_on[2])), TRUE,FALSE)
gene_table$strict_KD <- ifelse(((gene_table$percent_knockdown_MP_ontarget < 0 )& (gene_table$percent_knockdown_AP_ontarget < 0) ), TRUE,FALSE)
#save gene_table as a file
write.csv(gene_table,paste0(dir,"/gene_table.csv"))
}## [1] "gene_table.csv exists"Create heatmap of the percent knockdown for both promoters with background of on-target and off-target effects.
gene_table_heatmap <-gene_table %>%
dplyr::select(gene,percent_knockdown_MP_ontarget, percent_knockdown_AP_ontarget, percent_knockdown_MP_offtarget, percent_knockdown_AP_offtarget,percent_knockdown_MP_total,percent_knockdown_AP_total)
#put column to rowname
rownames(gene_table_heatmap) <- gene_table_heatmap$gene
gene_table_heatmap <- gene_table_heatmap %>% dplyr::select(-gene)
breaks2 <- seq(-100, 100, length.out = 150)
pheatmap(gene_table_heatmap %>%
dplyr::select(percent_knockdown_MP_ontarget,percent_knockdown_MP_offtarget,percent_knockdown_MP_total),
cluster_rows = TRUE, cluster_cols = FALSE,
show_rownames = TRUE, show_colnames = TRUE, main = "P1",
fontsize = 8, cellwidth = 8, cellheight = 8 , scale ="none" ,cutree_rows=5 ,colorRampPalette(brewer.pal(9,"Spectral"))(150),
breaks = breaks2)
#explain the number present
print(paste0("Number of genes in heatmap: ", nrow(gene_table_heatmap %>%
dplyr::select(percent_knockdown_MP_ontarget,percent_knockdown_MP_offtarget,percent_knockdown_MP_total) %>% drop_na()) ))## [1] "Number of genes in heatmap: 42"Repeat with P2-targeted guides. This is to show the difference in knockdown efficiency between the two promoters.
breaks2 <- seq(-100, 100, length.out = 150)
pheatmap(gene_table_heatmap %>%
dplyr::select(percent_knockdown_AP_ontarget,percent_knockdown_AP_offtarget,percent_knockdown_AP_total) %>% drop_na(),
cluster_rows = TRUE, cluster_cols = FALSE,
show_rownames = TRUE, show_colnames = TRUE, main = "P2",
fontsize = 8, cellwidth = 8, cellheight = 10 , scale ="none" ,cutree_rows=5 ,colorRampPalette(brewer.pal(9,"Spectral"))(150),
breaks = breaks2)
Comparison of the boolean assignments, looser KD, strict KD and significant NTC expression.
ggven_plot <- gene_table[c("looser_KD","strict_KD","sig_ntc_count")] %>%
dplyr::rename(
"Successful \n KD" = looser_KD,
"Strict \n Successful KD" = strict_KD,
"NTC Count\n >= 0.5 TPM" = sig_ntc_count
)
venn_diagram <- ggvenn(ggven_plot, show_stats="cp", text_size=3)
#save venndiagram
ggsave(paste0(rel_dir,"figures/venn_diagram.pdf"), plot = venn_diagram, width = 10, height = 10, units = "cm")
venn_diagram
Visualisation of the knockdown percentages for both promoters and both on-target and off-target effects.
gene_table_melt <- gene_table %>%
filter((looser_KD ==TRUE)) %>%
melt(id.vars = c("gene","neg_ctrl_tpm_AP_ontarget","neg_ctrl_tpm_MP_ontarget"),
measure.vars = c("percent_knockdown_MP_ontarget", "percent_knockdown_AP_ontarget",
"percent_knockdown_MP_offtarget", "percent_knockdown_AP_offtarget"))
gene_table_melt$target <- ifelse(grepl("ontarget", gene_table_melt$variable), "ontarget", "offtarget")
gene_table_melt$promoter <- ifelse(grepl("MP", gene_table_melt$variable), "P1", "P2")
ggplot(gene_table_melt) +
geom_boxplot(aes(x=variable, y=value, fill=variable)) +
theme_minimal() + labs(x = "Promoter", y = "Percent Knockdown") + #rotate
theme(axis.text.x = element_text(angle = 45, hjust = 1)) 
Plot the knockdown percentages for both promoters split by on-target and off-target for all 35 genes with NTC > 0.5.
a <- gene_table_melt %>% filter(promoter=="P1") %>%
ggviolin( x = "target", y = "value",fill = "target",add = "boxplot", ylab = "Percent Knockdown", xlab = "Promoter",
palette = c(palletee1[1],palletee1[2]),add.params = list(fill = "white") ) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means( paired = TRUE)
b <- gene_table_melt %>% filter(promoter=="P2") %>%
ggviolin( x = "target", y = "value",fill = "target",add = "boxplot", ylab = "Percent Knockdown", xlab = "Promoter",
palette = c(palletee2[1],palletee2[2]) ,add.params = list(fill = "white") ) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means( paired = TRUE)
plot_grid(a, b, ncol = 1)
Visualisation of the knockdown percentages for both promoters and both on-target and off-target effects.
gene_table_melt <- gene_table %>%
filter((looser_KD ==TRUE) & (sig_ntc_count==TRUE) ) %>%
melt(id.vars = c("gene","neg_ctrl_tpm_AP_ontarget","neg_ctrl_tpm_MP_ontarget"),
measure.vars = c("percent_knockdown_MP_ontarget", "percent_knockdown_AP_ontarget",
"percent_knockdown_MP_offtarget", "percent_knockdown_AP_offtarget"))
gene_table_melt$target <- ifelse(grepl("ontarget", gene_table_melt$variable), "ontarget", "offtarget")
gene_table_melt$promoter <- ifelse(grepl("MP", gene_table_melt$variable), "P1", "P2")
ggplot(gene_table_melt) +
geom_boxplot(aes(x=variable, y=value, fill=variable)) +
theme_minimal() + labs(x = "Promoter", y = "Percent Knockdown") + #rotate
theme(axis.text.x = element_text(angle = 45, hjust = 1)) 
Plot the knockdown percentages for both promoters split by on-target and off-target of the 31 genes that show “successful KD” which removes the top 15 quantiles.
a <- gene_table_melt %>% filter(promoter=="P1") %>%
ggviolin( x = "target", y = "value",fill = "target",add = "boxplot", ylab = "Percent Knockdown", xlab = "Promoter",
palette = c(palletee1[1],palletee1[2]),add.params = list(fill = "white") ) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means( paired = TRUE)
b <- gene_table_melt %>% filter(promoter=="P2") %>%
ggviolin( x = "target", y = "value",fill = "target",add = "boxplot", ylab = "Percent Knockdown", xlab = "Promoter",
palette = c(palletee2[1],palletee2[2]) ,add.params = list(fill = "white") ) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means( paired = TRUE)
plot_grid(a, b, ncol = 1)
Plot the knockdown percentages against the negative control TPM for both promoters split by on-target and off-target.
p<-ggscatter(gene_table_melt %>% filter(target=="ontarget"), x = "neg_ctrl_tpm_MP_ontarget", y = "value",
color = "promoter", line.color = "gray", line.size = 0.4,
palette = "npg")+
ggtitle("P1") + theme_minimal() + labs(y = "Percent Knockdown") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
a <-ggscatter(gene_table_melt %>% filter(target=="ontarget"), x = "neg_ctrl_tpm_MP_ontarget", y = "value",
color = "promoter", line.color = "gray", line.size = 0.4,
palette = "npg")+
ggtitle("P1") + theme_minimal() + labs(y = "Percent Knockdown") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
gene_table_melt %>% #filter for value between -100 and 100
melt(id.vars = c("gene"),measure.vars=c("neg_ctrl_tpm_AP_ontarget","neg_ctrl_tpm_MP_ontarget")) %>%
unique() %>%
ggplot(aes(x=value, fill=variable)) + #log scale
geom_density() +
labs(x = "Negative Control TPM", y = "Density") + theme_minimal() +
scale_x_continuous(trans='log10') + scale_fill_manual(values = c("#00AFBB", "#E7B800"))
Plot heatmap of the negative control TPMs for both promoters. Negative control TPMs represent the background expression of the targeted transcripts in the NTC population.
#create heatmap
matrix_ng_cntrl <- gene_table_melt %>% #filter for value between -100 and 100
dplyr::select( c("gene","neg_ctrl_tpm_MP_ontarget", "neg_ctrl_tpm_AP_ontarget")) %>%
unique()
rownames(matrix_ng_cntrl) <- matrix_ng_cntrl$gene
matrix_ng_cntrl <- as.matrix(matrix_ng_cntrl %>% dplyr::select(-c("gene")))
#make it so tthis is TPM
matrix_ng_cntrl <- log10(matrix_ng_cntrl + 1)
break2 <- seq(0, max(matrix_ng_cntrl), length.out = 150)
pheatmap(matrix_ng_cntrl, cluster_rows = TRUE, cluster_cols = FALSE,
show_rownames = TRUE, show_colnames = TRUE, main = "Negative Control",
fontsize = 8, cellwidth = 8, cellheight = 8 , scale ="none" ,cutree_rows=5 ,colorRampPalette(brewer.pal(9,"YlGnBu"))(150),
breaks = break2)
Plot the knockdown percentages for both promoters split by on-target only.
gene_table_melt %>% filter(target=="ontarget") %>%
ggviolin( x = "promoter", y = "value",fill = "promoter",add = "boxplot", ylab = "Percent Knockdown", xlab = "Promoter",
palette = c(palletee1[1],palletee2[1]),add.params = list(fill = "white") ) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) 
t.test(gene_table$sum_MP_ontarget, gene_table$neg_ctrl_tpm_MP_ontarget, paired = TRUE)##
## Paired t-test
##
## data: gene_table$sum_MP_ontarget and gene_table$neg_ctrl_tpm_MP_ontarget
## t = -3.8011, df = 41, p-value = 0.00047
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -270.44973 -82.77884
## sample estimates:
## mean difference
## -176.6143t.test(gene_table$sum_AP_ontarget, gene_table$neg_ctrl_tpm_AP_ontarget, paired = TRUE)##
## Paired t-test
##
## data: gene_table$sum_AP_ontarget and gene_table$neg_ctrl_tpm_AP_ontarget
## t = -1.8942, df = 41, p-value = 0.06526
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -30.7019295 0.9828819
## sample estimates:
## mean difference
## -14.85952Calculation of the effect size for the on-target P1 knockdown.
eff_size_mp <- gene_table %>%
dplyr::filter((looser_KD ==TRUE) & (sig_ntc_count==TRUE) ) %>%
wilcox_effsize( data = ., formula = sum_MP_ontarget ~ neg_ctrl_tpm_MP_ontarget, paired = TRUE)
eff_size_mp## # A tibble: 465 × 7
## .y. group1 group2 effsize n1 n2 magnitude
## * <chr> <chr> <chr> <dbl> <int> <int> <ord>
## 1 sum_MP_ontarget 1 14.2 1 1 1 large
## 2 sum_MP_ontarget 1 17.8 1 1 1 large
## 3 sum_MP_ontarget 1 17.9 1 1 1 large
## 4 sum_MP_ontarget 1 22 1 1 1 large
## 5 sum_MP_ontarget 1 25.2 1 1 1 large
## 6 sum_MP_ontarget 1 30.3 1 1 1 large
## 7 sum_MP_ontarget 1 45.8 1 1 1 large
## 8 sum_MP_ontarget 1 46.4 1 1 1 large
## 9 sum_MP_ontarget 1 59.7 1 1 1 large
## 10 sum_MP_ontarget 1 68.2 1 1 1 large
## # ℹ 455 more rowsAnother method to confirm base TPM expression in th enegative control population and the proportion of expression.
gene_table_melt_all <- gene_table %>%
melt(id.vars = c("gene","neg_ctrl_tpm_AP_ontarget","neg_ctrl_tpm_MP_ontarget","successfulKD"),
measure.vars = c("percent_knockdown_MP_ontarget", "percent_knockdown_AP_ontarget",
"percent_knockdown_MP_offtarget", "percent_knockdown_AP_offtarget"))
#plot the data negative control tpm
a <- gene_table_melt_all %>%
filter(successfulKD==TRUE) %>%
dplyr::select( c("gene","neg_ctrl_tpm_MP_ontarget", "neg_ctrl_tpm_AP_ontarget")) %>%
unique() %>%
melt(id.vars = c("gene"),measure.vars=c("neg_ctrl_tpm_MP_ontarget","neg_ctrl_tpm_AP_ontarget")) %>% #sort P1 to P2
mutate(promoter_type = ifelse(variable == "neg_ctrl_tpm_MP_ontarget", "P1", "P2"))
a$promoter_type <- factor(a$promoter_type, levels=c( "P2","P1"))
#change the scientific
a_plot_a <- a %>%
ggplot(aes(y=value, x=promoter_type, fill = promoter_type)) + #log scale
geom_violin() +
labs(y = "Sum TPM of Transcripts Under Promoter \n in NTC Population \n in Genes Assigned as Successful KD",x="Targeted Promoter") + yscale("log10")+
theme(axis.text.x = element_text(angle = 45, hjust = 1)) + #hide legend
theme(legend.position = "none")+
scale_fill_manual(values = c("#00AFBB", "#E7B800")) + #chaange labels for x axis
geom_boxplot(width = 0.2, fill = "white")+ #flip
coord_flip()
#plot the data negative control tpm
a <- gene_table_melt_all %>%
dplyr::select( c("gene", "successfulKD","neg_ctrl_tpm_MP_ontarget", "neg_ctrl_tpm_AP_ontarget")) %>%
unique() %>%
melt(id.vars = c("gene","successfulKD"),measure.vars=c("neg_ctrl_tpm_MP_ontarget","neg_ctrl_tpm_AP_ontarget")) %>% #sort P1 to P2
mutate(promoter_type = ifelse(variable == "neg_ctrl_tpm_MP_ontarget", "P1", "P2"))
a$promoter_type <- factor(a$promoter_type, levels=c( "P2","P1"))
# Count observations per group
counts <- a %>%
group_by(promoter_type, successfulKD) %>%
summarise(n = n(), .groups = "drop")
a_plot <- a %>%
ggplot(aes(y=value, x=promoter_type, fill = successfulKD)) + #log scale
geom_violin( position = position_dodge(0.9)) +
labs(y = "Sum TPM of Transcripts Under Promoter \n in NTC Population \n in Genes Assigned as Successful KD",x="Targeted Promoter") + yscale("log10")+
theme(axis.text.x = element_text(angle = 45, hjust = 1)) + #hide legend
scale_fill_manual(values = c(palletee1[1], palletee1[2],palletee2[1], palletee2[2])) + #chaange labels for x axis
geom_boxplot(width = 0.2, position = position_dodge(0.9) )+ #flip
coord_flip() +
geom_text(data = counts,
aes(x = promoter_type, y = 10000, label = paste0("n=", n)),
position = position_dodge(0.9), size = 3)
#change from scientific notation to visual
options(scipen=999)
p4 <- a %>% #filter for value between -100 and 100
ggdotchart( x = "gene", y = "value",
color = "promoter_type", # change fill color by mpg_level
palette = c("#00AFBB", "#E7B800"), # flip
orientation = "vertical", # change the orientation
yscale="log10" , xlab = "TPM of the Targeted Transcripts \n in NTC Population",
ylab="" ,
)
plot_grid(a_plot_a,a_plot, p4, ncol = 1, labels = c("A", "B", "C"))
Plot the cell counts and UMI counts per pseudobulk population to visualise the distribution of these values across all populations.
Plot the sum of TPMs per isoform across all samples to visualise the distribution of expression levels. This is to confirm that the data has enough depth to capture lowly expressed isoforms.
# Read or create the table
if (file.exists(table_name)) {
all_data <- read.table(table_name, header = TRUE, sep = "\t", row.names = NULL)
} else {
#map the fiunction process to every element in list of files and merge the otuputs to simplest df use applyr
all_data <- lapply(files, process_file)
#flatten and merge on isforom column
all_data <- Reduce(function(x, y) merge(x, y, by = "Isoform", all = TRUE), all_data)
#instead of cbind use merge on Isoform
write.table(all_data, table_name, sep = "\t", row.names = FALSE, quote = FALSE)
}
rowsums <- all_data %>% select(-Isoform) %>% rowSums() %>% as.data.frame() %>%
rename("Rowsums" = ".") %>% arrange(desc(Rowsums)) %>% #add 1
mutate(Rowsums = Rowsums + 1, logRowsums = log10(Rowsums))
rowsums <- rowsums %>%
mutate(index_col = as.integer(rownames(rowsums)), logIndex= log10(index_col))
#plot the logrowsums
ggplot(rowsums, aes(x = logIndex, y = logRowsums)) +
geom_line() +
theme_minimal() +
labs(x = "Log10 Index of Isoform (Descending Sum TPM)", y = "log10 Sum Isoform TPM") 