-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path19_total_expansion_counts.Rmd
141 lines (100 loc) · 5.42 KB
/
19_total_expansion_counts.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
---
title: "Untitled"
output: html_document
date: "2024-04-11"
---
```{r load data}
library(ggplot2)
library(data.table)
library(stringr)
library(viridis)
base_dir="/path/to/your/files"
setwd(base_dir)
#load mappings from mochi weight positions to pfam alignment positions
mochi_pos_to_pfam_alnpos<-fread("analysis_files/homolog_mochi_input_files/mochi_alnpos_to_pfam_alnpos.txt")
colnames(mochi_pos_to_pfam_alnpos)<-c("mochi_pos","aln_pos","PFAM_ID")
mochi_pos_to_pfam_alnpos<-mochi_pos_to_pfam_alnpos[!duplicated(mochi_pos_to_pfam_alnpos),]
#load pfam alignments to map position in protein to position in alignment
pfam_alignments<-fread("analysis_files/Pfam-A.human.seqpos_to_alnpos_predicted_mainiso")
colnames(pfam_alignments)<-c("PFAM_entry","pos","aln_pos","wt_aa","gene_ID","uniprot_ID","PFAM_ID",
"PFAM_ID.n","uniprot_ID_pos_in_uniprot")
#merge with mochi weights
mochi_pos_to_pfam_alnpos[,pfamid_alnpos:=paste(PFAM_ID,aln_pos,sep="_")]
pfamids<-fread("analysis_files/homolog_mochi_input_files/PFAM_IDs",header=FALSE)$V1
read_weights<-function(family,fit_type){
weights<-fread(paste("analysis_files/homolog_mochi_input_files/",family,"_",fit_type,"/weights/weights_Folding.txt",sep=""))
weights[,mut_aa:=substr(id,nchar(id),nchar(id))]
weights[id=="WT",mut_aa:=NA]
weights$mut_aa<-factor(weights$mut_aa,
levels=str_split("QNSTDEKRHGPCMAILVFYW","")[[1]])
features<-fread(paste("analysis_files/homolog_mochi_input_files/",family,"_features_solu.txt",sep=""))
wt_weights<-features[SoluWeight!="",]$SoluWeight
mut_weights<-weights[!is.na(`mean_kcal/mol`) & !(id %in% wt_weights),]
if (fit_type=="folding_linear"){
mut_weights$`mean_kcal/mol`<-mut_weights$`mean_kcal/mol`*(-1)
}
unfolded<-quantile(mut_weights$`mean_kcal/mol`,probs=0.975)
mut_weights[,mean_kcalmol_scaled:=`mean_kcal/mol`/unfolded]
mut_weights[,std_kcalmol_scaled:=`std_kcal/mol`/unfolded]
mut_weights<-mut_weights[,PFAM_ID:=family]
return(mut_weights)
}
homochi_weights<-data.table()
for (pfamid in pfamids){
homochi_weights<-rbind(homochi_weights,read_weights(pfamid,"folding"))
}
colnames(mochi_pos_to_pfam_alnpos)[1]<-"Pos"
homochi_weights_alnpos<-merge(homochi_weights,mochi_pos_to_pfam_alnpos,by=c("PFAM_ID","Pos"),
all.x=TRUE)
pfam_alignments[,pfamid_alnpos:=paste(PFAM_ID,aln_pos,sep="_")]
#merge
pfam_alignments_homochipredictions<-merge(pfam_alignments,homochi_weights_alnpos,by="pfamid_alnpos",
allow.cartesian = TRUE)
pfam_alignments_homochipredictions$mut_aa<-unlist(lapply(pfam_alignments_homochipredictions$id_ref,
FUN=function(string){
return(substr(string,nchar(string),nchar(string)))
}))
pfam_alignments_homochipredictions<-pfam_alignments_homochipredictions[mut_aa!=wt_aa,]
pfam_alignments_homochipredictions[,mut_ID:=paste(PFAM_entry,aln_pos.x,mut_aa,sep="_")]
length(unique(pfam_alignments_homochipredictions$mut_ID))
length(unique(pfam_alignments_homochipredictions$PFAM_entry))
length(unique(pfam_alignments_homochipredictions$mut_ID))/length(unique(pfam_alignments_homochipredictions$PFAM_entry))
expanded_counts<-data.table(table(pfam_alignments_homochipredictions$PFAM_ID.x))
#measured counts
mutated_domainome_filtered<-fread("analysis_files/mutated_domainome_merged_filtered.txt")
measured_counts<-data.table(table(mutated_domainome_filtered$PFAM_ID))
measured_counts<-measured_counts[V1 %in% expanded_counts$V1,]
expanded_counts[,source:="expanded"]
measured_counts[,source:="measured"]
all_counts<-rbind(expanded_counts,measured_counts)
colnames(all_counts)<-c("PFAM_ID","variant count","source")
all_counts_sum<-all_counts[,.(ct=sum(`variant count`)),by="PFAM_ID"]
all_counts_sum<-all_counts_sum[order(ct),]
all_counts$PFAM_ID<-factor(all_counts$PFAM_ID,
levels=all_counts_sum$PFAM_ID)
ggplot(all_counts)+
geom_col(aes(y=PFAM_ID,x=`variant count`,fill=source))
ggsave("output_files/ED_Figure9b_totalexpansioncount_variants.pdf")
#expansion for domains
expanded_counts_domains<-data.table(table(pfam_alignments_homochipredictions[!duplicated(PFAM_entry),]$PFAM_ID.x))
measured_counts_domains<-data.table(table(mutated_domainome_filtered[!duplicated(dom_ID),]$PFAM_ID))
measured_counts_domains<-measured_counts_domains[V1 %in% expanded_counts_domains$V1,]
expanded_counts_domains[,source:="expanded"]
measured_counts_domains[,source:="measured"]
all_counts_domains<-rbind(expanded_counts_domains,measured_counts_domains)
colnames(all_counts_domains)<-c("PFAM_ID","domain count","source")
all_counts_sum_domains<-all_counts_domains[,.(ct=sum(`domain count`)),by="PFAM_ID"]
all_counts_sum_domains<-all_counts_sum_domains[order(ct),]
all_counts_domains$PFAM_ID<-factor(all_counts_domains$PFAM_ID,
levels=all_counts_sum_domains$PFAM_ID)
ggplot(all_counts_domains)+
geom_col(aes(y=PFAM_ID,x=`domain count`,fill=source))
ggsave("output_files/ED_Figure9c_totalexpansioncount_domains.pdf")
#write table of predictions in homologs
supp_table_4<-pfam_alignments_homochipredictions[,c("PFAM_entry","uniprot_ID","gene_ID","pos","wt_aa","mut_aa","mean_kcalmol_scaled","std_kcalmol_scaled")]
write.table(supp_table_4,
file="output_files/ED_Table_4_homolog_predictions.txt",
quote=FALSE,
sep = "\t",
row.names=FALSE)
```