# script to set up protein blast database # this is for CAFA 6 competition # do not run this script - it is put here to show what I did on my local machine # it assumes that the additional data created from running notebook Cafa6_Protein_Data_Shape is available # https://www.kaggle.com/code/redberry/cafa6-protein-data-shape/notebook # I have also uploaded the addtional data and created a kaggle data set for this output # but it is not linked into this script properly # https://www.kaggle.com/datasets/redberry/additional-data-for-cafa6-protein-data-shape/data/settings/settings/settings/settings/settings/settings/settings # Sam Lycett # 23 Nov 2025 print("Define custom functions") ############################################################################################################ # Custom functions ############################################################################################################ # start define custom functions # read in output format 6 results, add column names and make into dataframe get_outfmt6_results <- function(results_file) { res <- read.table(results_file) colnames(res) <- c("query_seq","db_seq","perc_identity","seq_len","mismatches", "gap_opens","q_start","q_end","s_start","s_end","evalue","bit_score") res <- as.data.frame(res) return(res) } # split the sequence name into | and get the 2nd element which is the accession number in this format get_accn_from_seq_names <- function(seq_names,index=2) { els <- strsplit(seq_names,"\\|") nparts <- length(els[[1]]) nnames <- length(els) els <- matrix(unlist(els),nparts,nnames) if (index>0) { return(els[index,]) } else { return(t(els)) } } # get the train_terms_with_info rows for the specified protein accession number (multiple rows) get_go_terms_info_from_accn <- function(accn,train_terms_with_info=train_terms_with_info) { inds <- which(train_terms_with_info$EntryID==accn) return(t(train_terms_with_info[inds,])) } # match the query_accn to the blast results # look up the blast results go terms # use the top go terms and use the blast return percent identity as the probability # propograte the go terms up the ontology # return the predictions, which is query_accn, go_term, probability BLAST_basic_terms_matching <- function(query_accn, train_terms_with_info=train_terms_with_info, blast_results=blast_results, ontology=ontology, perc_identity_thres=25,evalue_thres=-1, useTop=1, useMention=2, return.prediction=TRUE, blast_results_has_accn=TRUE) { # extract the blast results relating to the query accession if (blast_results_has_accn) { blast_query <- blast_results$query_seq } else { blast_query <- get_accn_from_seq_names(blast_results$query_seq) } query_ind <- which(query_accn==blast_query) sel_res <- blast_results[query_ind,] if (perc_identity_thres>0) { inds <- which(sel_res$perc_identity>=perc_identity_thres) sel_res <- sel_res[inds,] } if (evalue_thres>0) { inds <- which(sel_res$evalue <= evalue_thres) sel_res <- sel_res[inds,] } num_retrieved <- length(sel_res$db_seq) if (num_retrieved > 0) { # get the retrieved sequences accessions from the query sequence blast retrieved_seq_accn <- get_accn_from_seq_names(sel_res$db_seq) retrieved_seq_prob <- sel_res$perc_identity/100 # use these as probabilities # get the GO terms of the retrieved blast sequences from the selected protein ? retrieved_info <- apply(as.matrix(retrieved_seq_accn),1,get_go_terms_info_from_accn,train_terms_with_info=train_terms_with_info) #els1 <- retrieved_info[[1]] #cn <- rownames(els1) cn <- c("EntryID","term","aspect","term_info","term_name","Organism","Scientific_Name","Type") retrieved_info <- unlist(retrieved_info) retrieved_info <- t(matrix(retrieved_info,length(cn),length(retrieved_info)/length(cn))) colnames(retrieved_info) <- cn retrieved_info <- as.data.frame(retrieved_info) prob <- array(perc_identity_thres/100,length(retrieved_info$EntryID)) for (j in 1:length(retrieved_seq_accn)) { inds <- which(retrieved_info$EntryID==retrieved_seq_accn[j]) prob[inds] <- retrieved_seq_prob[j] } retrieved_info <- cbind(retrieved_info,prob) retrieved_info <- as.data.frame(retrieved_info) ri_order <- sort(retrieved_info$prob,decreasing=TRUE,index.return=TRUE)$ix retrieved_info <- retrieved_info[ri_order,] # which terms are mentioned more than once in retrieved proteins mention_info <- retrieved_info[match(names(which(table(retrieved_info$term)>=useMention)),retrieved_info$term),] # which terms are from the top hit proteins if (useTop==1) { top_hit_proteins <- retrieved_seq_accn[1] top_inds <- which(retrieved_info$EntryID==top_hit_proteins) } else { top_hit_proteins <- retrieved_seq_accn[1:useTop] top_inds <- c() for (j in 1:length(top_hit_proteins)) { temp <- which(retrieved_info$EntryID==top_hit_proteins[j]) top_inds <- c(top_inds,temp) } } top_info <- retrieved_info[top_inds,] # combine and take the unique ones # match will match to the first in the object combined_info <- rbind(top_info,mention_info) if (!return.prediction) { uinds <- match(unique(combined_info$term),combined_info$term) combined_info <- combined_info[uinds,] } # propagate to root all_terms <- c() all_probs <- c() for (j in 1:length(combined_info$term)) { go_ancs <- ontology$ancestors[[combined_info$term[j]]] go_probs<- array(format(combined_info$prob[j],digits=3),length(go_ancs)) all_terms <- c(all_terms,go_ancs) all_probs <- c(all_probs,go_probs) } npreds <- length(all_terms) protein <- array(query_accn,npreds) preds <- cbind(protein,all_terms,all_probs) # make sure have unique go terms uterms <- unique(all_terms) # order is important uinds <- match(uterms,all_terms) preds <- preds[uinds,] colnames(preds) <- c("protein","term","probability") } else { print(paste("No good BLAST matches for",query_accn)) #BPO: GO:0008150 #CCO: GO:0005575 #MFO: GO:0003674 protein <- array(query_accn,3) all_terms <- c("GO:0008150","GO:0005575","GO:0003674") all_probs <- c(1,1,1) preds <- cbind(protein,all_terms,all_probs) colnames(preds) <- c("protein","term","probability") } if (return.prediction) { return(preds) } else { return(combined_info) } } # For any protein ID in the test superset, you must list a set of GO terms and assign your estimated probability. # If a protein ID is not listed in your submitted file, the organizers will assume that all predictions are 0. # The file should not contain a header; columns must be tab-separated. An example submission file may look as follows: # P9WHI7 GO:0009274 0.931 # P9WHI7 GO:0071944 0.540 # P9WHI7 GO:0005575 0.324 # P04637 GO:1990837 0.23 # P04637 GO:0031625 0.989 # P04637 GO:0043565 0.64 # P04637 GO:0001091 0.49 # end define custom functions ############################################################################################################ print("End define custom functions") print("Dont run the rest of this script") print("All the paths are wrong and need to be fixed") print("And you need to run BLAST locally anyway") doLoadLibraries <- FALSE if (doLoadLibraries) { library(seqinr) library(ontologyIndex) } doStep0 <- FALSE if (doStep0) { #################################### # STEP 0 - install BLAST program #################################### # # https://www.ncbi.nlm.nih.gov/books/NBK52640/ # https://www.ncbi.nlm.nih.gov/books/NBK279690/ # # Introduction # # NCBI provides command line standalone BLAST+ programs (based on the NCBI C++ toolkit) as a single compressed package. The package is available for the Linux, Mac OSX, and Windows platforms at: # # https://ftp.ncbi.nlm.nih.gov/blast/executables/LATEST/ # # Have already got Homebrew installed, so can do this on Terminal: # brew install blast # # Check by # blastn -help # blastp -help print("----------------") print("** STEP 0 **") print("Download and install BLAST from NCBI on the command line") print("----------------") } doStep1 <- FALSE if (doStep1) { #################################### # STEP 1 - make a database #################################### # # ~/Documents/data/ProteinFunction/cafa-6-protein-function-prediction/blast_database # put the train_sequences.fasta in here # # makeblastdb -in train_sequences.fasta -dbtype prot -out cafa_train_db print("----------------") print("** STEP 1 **") print("Make a blast database") blast_path <- "~/Documents/data/ProteinFunction/cafa-6-protein-function-prediction/blast_database/" setwd(blast_path) print("Use command line: (note this is done on the terminal, not within R)") print("makeblastdb -in train_sequences.fasta -dbtype prot -out cafa_train_db") print("Copy of output from terminal - ") screen_out_txt <- c( "Building a new DB, current time: 11/23/2025 15:09:49", "New DB name: Documents/data/ProteinFunction/cafa-6-protein-function-prediction/blast_database/cafa_train_db", "New DB title: train_sequences.fasta", "Sequence type: Protein", "Keep MBits: T", "Maximum file size: 3000000000B", "Adding sequences from FASTA; added 82404 sequences in 0.642918 seconds.") print(screen_out_txt) print("Files in directory") print(dir()) print("----------------") # now can remove train_sequences.fasta from this directory and just keep the created index files #[1] "cafa_train_db.pdb" "cafa_train_db.phr" "cafa_train_db.pin" "cafa_train_db.pjs" #[5] "cafa_train_db.pot" "cafa_train_db.psq" "cafa_train_db.ptf" "cafa_train_db.pto" } doStep2 <- FALSE if (doStep2) { #################################### # STEP 2 - test BLASTP #################################### # extract a few sequences from the train_sequences.fasta for testing purposes # library(seqinr) trainPath <- "~/Documents/data/ProteinFunction/cafa-6-protein-function-prediction/Train/" seqs <- read.fasta(paste0(trainPath,"train_sequences.fasta"), seqtype="AA", as.string=TRUE, forceDNAtolower = FALSE) print(paste("Number of training sequences = ",length(seqs))) # Number of training sequences = 82404 seqs_lens <- unlist(lapply(seqs,nchar)) print(paste("Min sequence length=",min(seqs_lens))) print(paste("Max sequence length=",max(seqs_lens))) print(paste("Median sequence length=",median(seqs_lens))) print(paste("Mean sequence length=",round(mean(seqs_lens)))) # Min sequence length= 3 # Max sequence length= 35213 # Median sequence length= 409 # Mean sequence length= 526 additionalPath <- "~/Documents/data/ProteinFunction/cafa-6-protein-function-prediction/additional_data/" train_terms_with_info <- read.csv(paste0(additionalPath,"train_terms_with_info.csv"),sep=";") train_terms_with_info[1,] #EntryID term aspect term_info term_name Organism Scientific_Name Type #1 Q5W0B1 GO:0000785 C 7.836723 chromatin 9606 Homo sapiens E # the first protein in the train_terms_with_info example_accn <- "Q5W0B1" ind <- grep(example_accn,attributes(seqs)$names) write.fasta(seqs[ind],attributes(seqs)$names[ind], file.out=paste0(blast_path,"query_seq.fasta")) # this will write to the blast path # the first 10 sequences in train ind <- 1:10 write.fasta(seqs[ind],attributes(seqs)$names[ind], file.out=paste0(blast_path,"query_10.fasta")) print(paste("Written example sequence",example_accn,"from the training set to the blast path as query_seq.fasta")) print("Written the first 10 sequences from the training set to the blast path as query_10.fasta") print(dir()) #"Written example sequence Q5W0B1 from the training set to the blast path as query_seq.fasta" #"Written the first 10 sequences from the training set to the blast path as query_10.fasta" #[1] "cafa_train_db.pdb" "cafa_train_db.phr" "cafa_train_db.pin" "cafa_train_db.pjs" "cafa_train_db.pot" #[6] "cafa_train_db.psq" "cafa_train_db.ptf" "cafa_train_db.pto" "query_seq.fasta" "query_10.fasta" print("from ~/Documents/data/ProteinFunction/cafa-6-protein-function-prediction/blast_database/ on command line do this:") system_call_txt <- "blastp -db cafa_train_db -query query_seq.fasta -out query_results.out.txt -outfmt 7 -max_target_seqs 10" print(system_call_txt) print("This can also be done from the R system call, assuming that blast is installed on the command line") # test system call # system("blastp -help") system(system_call_txt) print("Also trying the first 10 proteins (all in one file) with output format 6") system_call_txt <- "blastp -db cafa_train_db -query query_10.fasta -out query_10.out6.txt -outfmt 6 -max_target_seqs 10" print(system_call_txt) #system(system_call_txt) # read in the single example with output format 7 res <- readLines("query_results.out.txt") cn <- res[4] cn <- gsub("# Fields: ","",cn) cn <- strsplit(cn,", ")[[1]] # column names are cn: #[1] "query acc.ver" "subject acc.ver" "% identity" "alignment length" "mismatches" "gap opens" #[7] "q. start" "q. end" "s. start" "s. end" "evalue" "bit score" } doStep3 <- FALSE if (doStep3) { #################################### # STEP 3 - analyse results #################################### # read in results from query_10 res <- get_outfmt6_results("query_10.out6.txt") # look at the 2nd protein (first protein is too small) inds <- grep("A0JNW5",res$query_seq) sel_res <- res[inds,] print(sel_res) # query_seq db_seq perc_identity seq_len mismatches gap_opens q_start q_end s_start s_end evalue bit_score # 2 sp|A0JNW5|BLT3B_HUMAN sp|A0JNW5|BLT3B_HUMAN 100.000 1464 0 0 1 1464 1 1464 0.00e+00 3040.0 # 3 sp|A0JNW5|BLT3B_HUMAN sp|Q6BDS2|BLT3A_HUMAN 41.568 1518 744 32 1 1460 1 1433 0.00e+00 1027.0 # 4 sp|A0JNW5|BLT3B_HUMAN sp|Q54KX3|VP13F_DICDI 32.584 89 57 2 4 92 5 90 1.08e-06 55.5 # 5 sp|A0JNW5|BLT3B_HUMAN sp|Q80TY5|VP13B_MOUSE 33.333 78 48 2 9 85 10 84 5.34e-05 49.7 # 6 sp|A0JNW5|BLT3B_HUMAN sp|Q7Z7G8|VP13B_HUMAN 31.169 77 51 1 9 85 10 84 6.32e-05 49.7 # 7 sp|A0JNW5|BLT3B_HUMAN sp|Q55FG3|VP13C_DICDI 30.233 86 57 2 11 96 12 94 4.47e-04 46.6 # 8 sp|A0JNW5|BLT3B_HUMAN sp|Q5THJ4|VP13D_HUMAN 27.381 84 58 2 1 84 2 82 2.00e-03 44.7 # 9 sp|A0JNW5|BLT3B_HUMAN sp|A1Z713|VPS13_DROME 20.652 184 105 6 11 160 13 189 9.00e-03 42.4 # 10 sp|A0JNW5|BLT3B_HUMAN sp|G0S3B8|VPS13_CHATD 24.719 89 64 2 1 89 2 87 6.00e-02 39.7 # 11 sp|A0JNW5|BLT3B_HUMAN sp|Q5H8C4|VP13A_MOUSE 27.848 79 54 2 11 89 13 88 6.40e-02 39.7 # what are the GO terms of the retrieved blast sequences from the selected protein ? retrieved_seq_accn <- get_accn_from_seq_names(sel_res$db_seq) retrieved_info <- apply(as.matrix(retrieved_seq_accn),1,get_go_terms_info_from_accn,train_terms_with_info=train_terms_with_info) els1 <- retrieved_info[[1]] cn <- rownames(els1) retrieved_info <- unlist(retrieved_info) retrieved_info <- t(matrix(retrieved_info,length(cn),length(retrieved_info)/length(cn))) colnames(retrieved_info) <- cn retrieved_info <- as.data.frame(retrieved_info) # which terms are mentioned more than once in retrieved proteins retrieved_info[match(names(which(table(retrieved_info$term)>=2)),retrieved_info$term),] # EntryID term aspect term_info term_name Organism Scientific_Name Type # 7 A0JNW5 GO:0005515 F 3.968817 protein binding 9606 Homo sapiens E # 44 Q5THJ4 GO:0005739 C 10.781162 mitochondrion 9606 Homo sapiens E # 37 Q7Z7G8 GO:0007399 P 10.781162 nervous system development 9606 Homo sapiens E # 16 Q80TY5 GO:0035176 P 10.781162 social behavior 10090 Mus musculus E # 17 Q80TY5 GO:0035264 P 10.781162 multicellular organism growth 10090 Mus musculus E } doStep4 <- FALSE if (doStep4) { ################################################## # STEP 4 - make one prediction using functions ################################################## # suggest that use terms mentioned more than once + all terms of the top hit as a base line prediction # as baseline prediction # choose a query accession that is in the blast results query_accn <- "A0JNW5" # read blast results blast_results <- get_outfmt6_results("query_10.out6.txt") # read train terms with added information (see notebook) train_terms_with_info <- read.csv(paste0(additionalPath,"train_terms_with_info.csv"),sep=";") # read ontology file ontology <- get_ontology(paste0(trainPath,"go-basic.obo")) # format results into baseline predictions preds <- BLAST_basic_terms_matching(query_accn, train_terms_with_info = train_terms_with_info, blast_results=blast_results, ontology=ontology, blast_results_has_accn = FALSE) # output results from one protein write.table(preds,file=paste0(query_accn,"_BLAST_basic_terms_matching_example.txt"),sep="\t",row.names=FALSE,col.names=FALSE,quote=FALSE) } doStep5 <- FALSE if (doStep5) { #################################################################################### # STEP 5 - do blast of some test proteins and check that can do predictions of one #################################################################################### print("from ~/Documents/data/ProteinFunction/cafa-6-protein-function-prediction/blast_database/ on command line do this:") system_call_txt <- "blastp -db cafa_train_db -query testsuperset.fasta -out testsuperset_blast_results.out6.txt -outfmt 6 -max_target_seqs 10" print(system_call_txt) print("This can also be done from the R system call, assuming that blast is installed on the command line") # test system call # system("blastp -help") # system(system_call_txt) print("Actually this was only run for some of the test proteins") # actually just did this partial - but have a look for testing # read blast results blast_results <- get_outfmt6_results("testsuperset_blast_results.out6.txt") # read train terms with added information (see notebook) train_terms_with_info <- read.csv(paste0(additionalPath,"train_terms_with_info.csv"),sep=";") # read ontology file ontology <- get_ontology(paste0(trainPath,"go-basic.obo")) # choose a query accession that is in the blast results query_accn <- "A2A2Y4" # format results into baseline predictions preds <- BLAST_basic_terms_matching(query_accn, train_terms_with_info = train_terms_with_info, blast_results=blast_results, ontology=ontology, blast_results_has_accn = TRUE) # output results from one protein write.table(preds,file=paste0(query_accn,"_BLAST_basic_terms_matching_example.txt"),sep="\t",row.names=FALSE,col.names=FALSE,quote=FALSE) } #################################################################################### # STEP 6 - do blast of test proteins in blocks #################################################################################### doStep6 <- FALSE if (doStep6) { # extract sequences from the testsuperset.fasta # library(seqinr) blast_path <- "~/Documents/data/ProteinFunction/cafa-6-protein-function-prediction/blast_database/" trainPath <- "~/Documents/data/ProteinFunction/cafa-6-protein-function-prediction/Train/" testPath <- "~/Documents/data/ProteinFunction/cafa-6-protein-function-prediction/Test/" additionalPath <- "~/Documents/data/ProteinFunction/cafa-6-protein-function-prediction/additional_data/" seqs <- read.fasta(paste0(testPath,"testsuperset.fasta"), seqtype="AA", as.string=TRUE, forceDNAtolower = FALSE) print(paste("Number of test sequences = ",length(seqs))) # Number of test sequences = 224309 seqs_lens <- unlist(lapply(seqs,nchar)) print(paste("Min sequence length=",min(seqs_lens))) print(paste("Max sequence length=",max(seqs_lens))) print(paste("Median sequence length=",median(seqs_lens))) print(paste("Mean sequence length=",round(mean(seqs_lens)))) # Min sequence length= 2 # Max sequence length= 35213 # Median sequence length= 342 # Mean sequence length= 429 very_short_seqs <- which(seqs_lens <= 11) short_seqs <- which(seqs_lens > 11 & seqs_lens <= 42 ) medium_seqs <- which(seqs_lens > 42 & seqs_lens <= 2000 ) long_seqs <- which(seqs_lens > 2000) print(paste("Number very short seqs (<=11) =",length(very_short_seqs))) print(paste("Number short seqs (>11 & <= 42) =",length(short_seqs))) print(paste("Number medium seqs (>42 & <= 2000) =",length(medium_seqs))) print(paste("Number long seqs (>2000) =",length(long_seqs))) # Number very short seqs (<=11) = 1599 # Number short seqs (>11 & <= 42) = 6629 # Number medium seqs (>42 & <= 2000) = 213752 # Number long seqs (>2000) = 2329 # write the very short, short and long sequences to one file each if (!file.exists(paste0(blast_path,"testsuperset_very_short.fasta"))) { write.fasta(seqs[very_short_seqs], attributes(seqs)$names[very_short_seqs], file.out=paste0(blast_path,"testsuperset_very_short.fasta")) } else { print("already written testsuperset_very_short.fasta") } if (!file.exists(paste0(blast_path,"testsuperset_short.fasta"))) { write.fasta(seqs[short_seqs], attributes(seqs)$names[short_seqs], file.out=paste0(blast_path,"testsuperset_short.fasta")) } else { print("already written testsuperset_short.fasta") } if (!file.exists(paste0(blast_path,"testsuperset_long.fasta"))) { write.fasta(seqs[long_seqs], attributes(seqs)$names[long_seqs], file.out=paste0(blast_path,"testsuperset_long.fasta")) } else { print("already written testsuperset_long.fasta") } # do blasts on the command line (not in R) with these commands system_call_txt_1 <- "blastp -db cafa_train_db -query testsuperset_very_short.fasta -out testsuperset_very_short_blast_results.out6.txt -outfmt 6 -max_target_seqs 10" system_call_txt_2 <- "blastp -db cafa_train_db -query testsuperset_short.fasta -out testsuperset_short_blast_results.out6.txt -outfmt 6 -max_target_seqs 10" system_call_txt_3 <- "blastp -db cafa_train_db -query testsuperset_long.fasta -out testsuperset_long_blast_results.out6.txt -outfmt 6 -max_target_seqs 10" print(paste("Number medium seqs (>42 & <= 2000) =",length(medium_seqs))) print("This is likely too many to do all in one go - so split into smaller files") nseqs_per_batch <- 5120 nbatch <- ceiling(length(medium_seqs)/nseqs_per_batch) print(paste("number of medium sequences per batch = ",nseqs_per_batch)) print(paste("number of batches =",nbatch)) #"Number medium seqs (>42 & <= 2000) = 213752" #"This is likely too many to do all in one go - so split into smaller files" #"number of medium sequences per batch = 5120" #"number of batches = 42" is <- seq(1,length(medium_seqs),nseqs_per_batch) ie <- c(is[2:length(is)]-1,length(medium_seqs)) medium_snames <- paste0("testsuperset_medium_batch",1:nbatch,".fasta") for (i in 1:nbatch) { med_inds <- medium_seqs[is[i]:ie[i]] if (!file.exists(paste0(blast_path,medium_snames[i]))) { write.fasta(seqs[med_inds], attributes(seqs)$names[med_inds], file.out=paste0(blast_path,medium_snames[i])) } else { print(paste("already written",medium_snames[i])) } } setwd(blast_path) # try blasting these with system call - this works but is very slow # so try instead making shell file shell_file <- "do_medium_batches_blast.sh" for (i in 1:nbatch) { out_file <- gsub(".fasta","_blast_results.out6.txt",medium_snames[i]) system_call_txt_m <- paste0( "blastp -db cafa_train_db -query ",medium_snames[i], " -out ",out_file, " -outfmt 6 -max_target_seqs 10") write(system_call_txt_m, file=shell_file, append=(i>1)) #if (!file.exists(out_file)) { # t1 <- Sys.time() # print(system_call_txt_m) # system(system_call_txt_m) # t2 <- Sys.time() # print(paste("BLAST batch",i," in ",t2-t1)) #} else { # print(paste("already done blast and have results in",out_file)) #} } print("Created shell file with medium blast commands in") print("Run on the terminal command line:") print(paste("sh",shell_file)) #"sh do_medium_batches_blast.sh" ##################################### # read train terms with added information (see notebook) train_terms_with_info <- read.csv(paste0(additionalPath,"train_terms_with_info.csv"),sep=";") # read ontology file ontology <- get_ontology(paste0(trainPath,"go-basic.obo")) blast_runs <- c("very_short","short","long") blast_results_filename <- paste0("testsuperset_",blast_runs,"_blast_results.out6.txt") preds_filename <- paste0("testsuperset_",blast_runs,"_BLAST_basic_terms_matching_preds.txt") setwd(blast_path) for (j in 1:length(blast_results_filename)) { if (file.exists(blast_results_filename[j])) { if (!file.exists(preds_filename[j])) { blast_results <- get_outfmt6_results(blast_results_filename[j]) uquery <- unique(blast_results$query_seq) print(paste("Processing",blast_results_filename[j],"number of queries=",length(uquery))) for (i in 1:length(uquery)) { # format results into baseline predictions query_accn <- uquery[i] preds <- BLAST_basic_terms_matching(query_accn, train_terms_with_info = train_terms_with_info, blast_results=blast_results, ontology=ontology, blast_results_has_accn = TRUE) write.table(preds, file=preds_filename[j], sep="\t",row.names=FALSE,col.names=FALSE,quote=FALSE, append=(i>1)) } } else { print(paste("Already processed blast results",blast_results_filename[j],"to predictions",preds_filename[j])) } } else { print(paste("Cannot find blast results file ",blast_results_filename[j])) } } }