setwd("E:/ONT/RUN20_beeITS_refDat_XLONG/") library(ShortRead) library(Biostrings) library(magrittr) library(ape) library(umap) library(dbscan) vsearch_exec <- "D:/R/vsearch.exe" number_of_cores <- 20 vmatchPattern2 <- function(pattern, subject, max.mismatch=0, min.mismatch=0, with.indels=FALSE, fixed=TRUE, algorithm="auto") { if (!is(subject, "XStringSet")) subject <- Biostrings:::XStringSet(NULL, subject) algo <- Biostrings:::normargAlgorithm(algorithm) if (Biostrings:::isCharacterAlgo(algo)) stop("'subject' must be a single (non-empty) string ", "for this algorithm") pattern <- Biostrings:::normargPattern(pattern, subject) max.mismatch <- Biostrings:::normargMaxMismatch(max.mismatch) min.mismatch <- Biostrings:::normargMinMismatch(min.mismatch, max.mismatch) with.indels <- Biostrings:::normargWithIndels(with.indels) fixed <- Biostrings:::normargFixed(fixed, subject) algo <- Biostrings:::selectAlgo(algo, pattern, max.mismatch, min.mismatch, with.indels, fixed) C_ans <- .Call2("XStringSet_vmatch_pattern", pattern, subject, max.mismatch, min.mismatch, with.indels, fixed, algo, "MATCHES_AS_RANGES", PACKAGE="Biostrings") unlisted_ans <- IRanges(start=unlist(C_ans[[1L]], use.names=FALSE), width=unlist(C_ans[[2L]], use.names=FALSE)) relist(unlisted_ans, C_ans[[1L]]) } x <- readFastq(dirPath = "./fastq", full=FALSE) #length: 685038 reads; width: 3..61636 cycles x <- x[width(x) > 260] PB1F <- "GCAGTCGAACATGTAGCTGACTCAGGTCAC" PB1R <- "TGGATCACTTGTGCAAGCATCACATCGTAG" # ITSp5 <- "CCTTATCAYTTAGAGGAAGGAG" # R2876 <- "GTCAGACGAACGATTTGCACGTCAGTATC" LSU_end_Apocrita_R <- "GGTTTCCAGGGAACTCGAACGCTTATAC" #REVERSE SSU_end_Ecdysozoa_F <- "TAGGTGAACCTGCGGAAGGATCATTA" #FORWARD #start_Hym_R <- "TATGCTTAAATTCAGCGGGTRRTCTCG" #ITS2_1522_mid_consF <- "GRCAYYWTGGYGVGGTR" #ExEu_ITS2_1522F_PB1F GAGRCACTTTGGYGRGGTGM #Euf_ITS2_1522F_PB1F --GGCACCTTGGTGCGGTAC #Eug_ITS2_1522F_PB1F AAGRCAYTATGGTGCGGTG- x <- x[!x@id %>% duplicated] #if primer changes, check fixed = T temp_forw_match <- vmatchPattern2(SSU_end_Ecdysozoa_F %>% DNAString %>% reverseComplement, x@sread %>% narrow(., start = (x@sread %>% width - 250), end = (x@sread %>% width)), max.mismatch = 4, fixed = F, with.indels= T) fwd_matches_reverse <- as.data.frame(temp_forw_match)[,"group"] temp_rev_match <- vmatchPattern2(LSU_end_Apocrita_R %>% DNAString, x@sread %>% narrow(., start = 1, width = 250), max.mismatch = 4, fixed = F, with.indels= T) rev_matches_forw <- as.data.frame(temp_rev_match)[,"group"] rev_candidates <- fwd_matches_reverse[fwd_matches_reverse %in% rev_matches_forw] temp_forw_match <- vmatchPattern2(SSU_end_Ecdysozoa_F %>% DNAString, x@sread %>% narrow(., start = 1, width = 250), max.mismatch = 4, fixed = F, with.indels= T) forw_matches_forw <- as.data.frame(temp_forw_match)[,"group"] temp_rev_match <- vmatchPattern2(LSU_end_Apocrita_R %>% DNAString %>% reverseComplement, x@sread %>% narrow(., start = (x@sread %>% width - 250), end = (x@sread %>% width)), max.mismatch = 4, fixed = F, with.indels= T) rev_matches_revcomp <- as.data.frame(temp_rev_match)[,"group"] forw_candidates <- forw_matches_forw[forw_matches_forw %in% rev_matches_revcomp] rev_candidates <- rev_candidates[rev_candidates %in% forw_candidates %>% not] y <- append(x[forw_candidates %>% unique], x[rev_candidates %>% unique] %>% reverseComplement) y <- y[( alphabetScore(y) / width(y) ) > 10] #length: 481057 reads; width: 261..30914 cycles #allow only reads with exactly ONE PB1R and ONE PB1F tail temp_forw_match <- vmatchPattern2(PB1F %>% DNAString, y@sread, max.mismatch = 4, fixed = F, with.indels= T) forw_matches_forw <- as.data.frame(temp_forw_match)[,"group"] y[forw_matches_forw[forw_matches_forw %>% duplicated]] y <- y[-forw_matches_forw[forw_matches_forw %>% duplicated]] temp_rev_match <- vmatchPattern2(PB1R %>% DNAString %>% reverseComplement, y@sread, max.mismatch = 4, fixed = F, with.indels= T) rev_matches_revcomp <- as.data.frame(temp_rev_match)[,"group"] y[rev_matches_revcomp[rev_matches_revcomp %>% duplicated]] y <- y[-rev_matches_revcomp[rev_matches_revcomp %>% duplicated]] y #length: 469254 reads; width: 261..16264 cycles width(y) %>% plot if(file.exists("no_demult.fastq")) file.remove("no_demult.fastq") writeFastq(y, "no_demult.fastq", compress=F) temp_index <- read.csv("demult24x96.csv", sep=";") forw_list <- temp_index[,c("Forward.UMI","Forward.Seq")] %>% unique rev_list <- temp_index[,c("Reverse.UMI","Reverse.Seq")] %>% unique save <- y for (i in 1:nrow(forw_list)) { if (i %% 20 == 0) print(i) temp_forw_match <- vmatchPattern2(forw_list[i,"Forward.Seq"] %>% DNAString, y@sread %>% narrow(., start=1, end=100), max.mismatch = 2, fixed = T, with.indels= T) forw_match <- as.data.frame(temp_forw_match)[,"group"] forw_match %<>% unique y@id[forw_match] <- y@id[forw_match] %>% paste0(.,"_forwbc:",forw_list[i,"Forward.UMI"]) %>% BStringSet } for (i in 1:nrow(rev_list)) { temp_rev_match <- vmatchPattern2(rev_list[i,"Reverse.Seq"] %>% DNAString, y@sread %>% narrow(., start = (y@sread %>% width - 100), end = (y@sread %>% width) ), max.mismatch = 2, fixed = T, with.indels= T) rev_match <- as.data.frame(temp_rev_match)[,"group"] rev_match %<>% unique y@id[rev_match] <- y@id[rev_match] %>% paste0(.,"_revbc:",rev_list[i,"Reverse.UMI"]) %>% BStringSet } #filter for two barcodes found y <- y[( gregexpr("forwbc", y@id) %>% lengths == 1 ) & ( gregexpr("revbc", y@id) %>% lengths == 1 )] #exactly one rev AND one forw barcode found y <- y[( gregexpr("forwbc", y@id) %>% lapply(.,`[[`,1) %>% unlist != -1 ) & ( gregexpr("revbc", y@id) %>% lapply(.,`[[`,1) %>% unlist != -1 ) ] #at least one reverse and one forward barcode y #length: 237168 reads; width: 261..14957 cycles temp_forw <- y@id %>% gsub(".*forwbc:","",.) %>% gsub("_revbc.*|_exp.*","",.) %>% gsub("^[0-9a-z].*-.*","",.) %>% table %>% as.data.frame temp_rev <- y@id %>% gsub(".*_revbc:","",.) %>% gsub("_forwbc:.*","",.) %>% gsub("^[0-9a-z].*-.*","",.) # reverse_barcodes <- c("bc1108_rc") # # y <- y[temp_rev %in% reverse_barcodes] y <- y[( y@id %>% gsub(".*_revbc:","",.) %>% gsub("_forwbc:.*","",.) %>% gsub("^[0-9a-z].*-.*","",.) ) %in% ( c(paste0("bc110",1:9,"_rc"), paste0("bc11",10:24,"_rc")) )] temp_rev <- y@id %>% gsub(".*_revbc:","",.) %>% gsub("_forwbc:.*","",.) %>% gsub("^[0-9a-z].*-.*","",.) %>% table %>% as.data.frame merge(temp_forw %>% as.data.frame, temp_rev %>% as.data.frame, all=T, by=".") # . Freq.x Freq.y # 1 bc1001 1839 NA # 2 bc1002 1143 NA # 3 bc1003 4696 NA # 4 bc1004 3203 NA # 5 bc1005 5065 NA # 6 bc1006 2254 NA # 7 bc1007 1865 NA # 8 bc1008 2947 NA # 9 bc1009 1959 NA # 10 bc1010 1716 NA # 11 bc1011 1752 NA # 12 bc1012 1365 NA # 13 bc1013 1801 NA # 14 bc1014 1887 NA # 15 bc1015 2772 NA # 16 bc1016 2969 NA # 17 bc1017 3276 NA # 18 bc1018 2766 NA # 19 bc1019 2060 NA # 20 bc1020 2792 NA # 21 bc1021 2540 NA # 22 bc1022 2856 NA # 23 bc1023 2612 NA # 24 bc1024 2003 NA # 25 bc1025 4321 NA # 26 bc1026 2545 NA # 27 bc1027 4088 NA # 28 bc1028 2664 NA # 29 bc1029 2212 NA # 30 bc1030 1965 NA # 31 bc1031 1816 NA # 32 bc1032 1038 NA # 33 bc1033 2798 NA # 34 bc1034 4754 NA # 35 bc1035 2682 NA # 36 bc1036 2974 NA # 37 bc1037 2635 NA # 38 bc1038 2574 NA # 39 bc1039 4290 NA # 40 bc1040 2133 NA # 41 bc1041 2209 NA # 42 bc1042 1681 NA # 43 bc1043 2648 NA # 44 bc1044 4696 NA # 45 bc1045 1283 NA # 46 bc1046 2640 NA # 47 bc1047 818 NA # 48 bc1048 3064 NA # 49 bc1049 284 NA # 50 bc1050 1754 NA # 51 bc1051 1749 NA # 52 bc1052 2460 NA # 53 bc1053 2101 NA # 54 bc1054 2926 NA # 55 bc1055 2972 NA # 56 bc1056 2344 NA # 57 bc1057 2425 NA # 58 bc1058 2538 NA # 59 bc1059 2332 NA # 60 bc1060 1709 NA # 61 bc1061 3098 NA # 62 bc1062 2545 NA # 63 bc1063 1875 NA # 64 bc1064 3049 NA # 65 bc1065 245 NA # 66 bc1066 3535 NA # 67 bc1067 1575 NA # 68 bc1068 3011 NA # 69 bc1069 2129 NA # 70 bc1070 2497 NA # 71 bc1071 2705 NA # 72 bc1072 2951 NA # 73 bc1073 2755 NA # 74 bc1074 2606 NA # 75 bc1075 2959 NA # 76 bc1076 3554 NA # 77 bc1077 2771 NA # 78 bc1078 3225 NA # 79 bc1079 1838 NA # 80 bc1080 2686 NA # 81 bc1081 2245 NA # 82 bc1082 2402 NA # 83 bc1083 2039 NA # 84 bc1084 2632 NA # 85 bc1085 2064 NA # 86 bc1086 2856 NA # 87 bc1087 2694 NA # 88 bc1088 3624 NA # 89 bc1089 2728 NA # 90 bc1090 1808 NA # 91 bc1091 2277 NA # 92 bc1092 1528 NA # 93 bc1093 280 NA # 94 bc1094 1646 NA # 95 bc1095 1978 NA # 96 bc1096 2187 NA # 97 bc1101_rc NA 28549 # 98 bc1102_rc NA 19078 # 99 bc1103_rc NA 25729 # 100 bc1104_rc NA 24880 # 101 bc1105_rc NA 116304 # 102 bc1106_rc NA 22284 #PLATE1 UMI 12 #PLATE2 UMI 7 #PLATE3 UMI 8 y <- y[grep("bc1101_rc|bc1102_rc|bc1103_rc|bc1104_rc|bc1105_rc|bc1106_rc", y@id)] #length: 236824 reads; width: 261..14957 cycles save <- y temp_forw_match <- vmatchPattern2(SSU_end_Ecdysozoa_F %>% DNAString, y@sread %>% narrow(., start = 1, width = 200), max.mismatch = 4, fixed = F, with.indels= T) forw_matches_forw <- as.data.frame(temp_forw_match) forw_matches_forw <- forw_matches_forw[forw_matches_forw[,"group"] %>% duplicated %>% not,] y <- y[forw_matches_forw[,"group"]] y <- narrow(y, start = pmax(1,forw_matches_forw[,"start"]+nchar(SSU_end_Ecdysozoa_F)), end = width(y)) y <- y[width(y) > 200] temp_rev_match <- vmatchPattern2(LSU_end_Apocrita_R %>% DNAString %>% reverseComplement, y@sread %>% narrow(., start = (y@sread %>% width - 200), end = (y@sread %>% width)), max.mismatch = 4, fixed = F, with.indels= T) rev_matches_revcomp <- as.data.frame(temp_rev_match) rev_matches_revcomp <- rev_matches_revcomp[rev_matches_revcomp[,"group"] %>% duplicated %>% not,] y <- y[rev_matches_revcomp[,"group"]] y <- narrow(y, start = 1, end = pmin ( rev_matches_revcomp[,"end"] - nchar(LSU_end_Apocrita_R) + (y@sread %>% width) - 200, (y@sread %>% width) ) ) y@id %<>% gsub("\t.*_forw","_forw",.) %>% BStringSet y %>% width %>% sort %>% plot #y[grep("bc1119_rc", y@id)] %>% width %>% plot if(file.exists("orient.fastq")) file.remove("orient.fastq") writeFastq(y, file="orient.fastq", compress=F) if(file.exists("orient.fasta")) file.remove("orient.fasta") writeFasta(y, file="orient.fasta", compress=F) #( alphabetScore(y) / width(y) ) %>% sort %>% barplot y <- readFasta(dirPath = ".", pattern = "orient.fasta$") y@id <- y@id %>% gsub("-","",.) %>% gsub("^(.*?) .*forwbc:","\\1_forwbc:",.) %>% BStringSet y <- y[y@id %>% gsub(".*_forwbc:","",.) %>% order(y@id %>% gsub(".*_revbc:","",.),.)] #check exclusion critera y %>% width %>% plot abline(h=2500, col="red") abline(h=10000, col="red") y <- y[width(y)>5000 & width(y)< 10000] #check specific plate success #y[grep("revbc:bc1101_rc", y@id)]@id %>% gsub(".*_forwbc:","",.) %>% table temp_forw <- y@id %>% gsub(".*forwbc:","",.) %>% gsub("_revbc.*|_exp.*","",.) %>% gsub("^[0-9a-z].*-.*","",.) temp_rev <- y@id %>% gsub(".*_revbc:","",.) %>% gsub("_forwbc:.*","",.) %>% gsub("^[0-9a-z].*-.*","",.) umis_all <- paste0(temp_forw, "__",temp_rev) umis <- umis_all %>% unique umis <- umis[umis!=""] umis <- umis[order( umis %>% gsub(".*__","",.) , umis ) ] # i <- which(umis =="bc1031__bc1112_rc") #over split # i <- which(umis =="bc1090__bc1112_rc") #difficult homopolymer region for (i in 1:length(umis)) { temp_hits <- y[grepl(umis_all, pattern=umis[i])] if (temp_hits %>% length %>% sum < 5) next temp_hits <- temp_hits[order(temp_hits %>% .@sread %>% alphabetFrequency(.,as.prob=T) %>% .[,c(2,3)] %>% rowSums, temp_hits %>% width)] if (temp_hits %>% length > 500) temp_hits <- temp_hits %>% sample(500) writeFasta(temp_hits, file=paste0("./smp2/", umis[i], ".fasta")) #if (file.exists(paste0("./smp2/", umis[i], "_aln.fasta"))) next system(command = paste0(vsearch_exec, " --cluster_size ", paste0("./smp2/", umis[i], ".fasta"), " --id 0.9 --threads ",number_of_cores, " --consout ", paste0("./smp2/", umis[i], "_cons.fasta")," --uc ", paste0("./smp2/", umis[i], "_cluster.csv"), " --iddef 1 --clusterout_id" )) temp_cons <- readDNAStringSet(paste0("./smp2/", umis[i], "_cons.fasta")) temp_seq <- readDNAStringSet(paste0("./smp2/", umis[i], ".fasta")) temp_clusterinfo <- read.csv(paste0("./smp2/", umis[i], "_cluster.csv"), sep="\t", header = F) names(temp_seq) %<>% gsub(" .*","",.) temp_keepclusterid <- which(temp_clusterinfo[,"V2"] %>% table > 3) %>% names %>% as.integer #each cluster counts twice minimum if (length(temp_seq) < 10) temp_keepclusterid <- which(temp_clusterinfo[,"V2"] %>% table > 2) %>% names %>% as.integer if (temp_keepclusterid %>% length == 0) print(paste0(i, " No sequences left after clustering.")) temp_keepseqnames <- temp_clusterinfo[temp_clusterinfo[,"V2"] %in% temp_keepclusterid,"V9"] temp_select <- temp_seq[temp_seq %>% names %in% temp_keepseqnames] writeFasta(temp_select, file=paste0("./smp2/", umis[i], "_refined.fasta")) system.time( res <- system2(command="D:/R/mafft-win/mafft.bat", args=paste0(" --retree 2 --maxiterate 0 --thread 12 ", "--threadtb 12 --threadit 12 --op 1 --ep 0.12 --nomemsave --kimura 1 ", "--reorder ", paste0("./smp2/", umis[i], "_refined.fasta")," > ",paste0("./smp2/", umis[i], "_aln.fasta"))) ) } for (i in 1:length(umis)) { if (! file.exists(paste0("./smp2/", umis[i], "_aln.fasta")) ) next temp_hits <- readDNAStringSet( paste0("./smp2/", umis[i], "_aln.fasta") ) temp_start <- length(temp_hits) #check for conserved regions, only keep sequences which show these regions mid <- "gaactgcaggacacatgaacatcgacatttcgaacgcacattgcggtccacggatacaattcctggaccacgcctggctgagggtcgttcacgtaacccatatactgcttgcgttgc" end <- "ttgagacgacctcagagtaggcgagattacccgctgaatttaagcatattawtaagcggaggaaaagaaactaacaaggatttccttagtagcg" temp_forw_match <- vmatchPattern2(mid %>% DNAString, temp_hits %>% gsub("-","",.) %>% DNAStringSet, max.mismatch = round(nchar(mid)/10), fixed = F, with.indels= T) forw_matches_forw <- as.data.frame(temp_forw_match) forw_matches_forw <- forw_matches_forw[forw_matches_forw[,"group"] %>% duplicated %>% not,] temp_hits <- temp_hits[forw_matches_forw[,"group"]] temp_forw_match <- vmatchPattern2(end %>% DNAString, temp_hits %>% gsub("-","",.) %>% DNAStringSet, max.mismatch = round(nchar(end)/10), fixed = F, with.indels= T) forw_matches_forw <- as.data.frame(temp_forw_match) forw_matches_forw <- forw_matches_forw[forw_matches_forw[,"group"] %>% duplicated %>% not,] temp_hits <- temp_hits[forw_matches_forw[,"group"]] temp_stop <- length(temp_hits) paste0( round( (1-temp_stop/temp_start) *100), "% reads removed, umi: ", i) %>% print if (length(temp_hits)<5) next writeXStringSet(temp_hits, filepath= paste0("./smp2/", umis[i], "_alnfilt.fasta")) } repeat_info <- data.frame(matrix(ncol=10)) colnames(repeat_info) <- c("plate","sample", "coords", "seq", "start", "stop", "score", "period", "copies", "cons") for (i in 1:length(umis)) { sum( list.files(path = "./temp", pattern= ".*html$", full.names = T) %>% file.remove(.) ) res <- system2(command="./temp/trf409.dos64.exe", args= paste0( " ./smp2/", umis[i], "_alnfilt.fasta 2 7 5 80 10 250 1500"), stderr = NA, stdout = NULL ) sum( list.files(path = "./smp2", pattern= ".*html$", full.names = T) %>% file.copy(from=., to = "./temp/") ) sum( list.files(path = "./smp2", pattern= ".*html$", full.names = T) %>% file.remove(.) ) temp_files <- list.files(path = "./temp", pattern= ".txt.html$") if (length(temp_files)==0) next for (j in 1:length(temp_files)) { temp <- readLines( paste0("./temp/", temp_files[j]), skipNul = TRUE,warn = FALSE ,encoding = "UTF-8") temp_hits <- grep("Found at i", temp) if (length(temp_hits)==0) next temp_hits <- c(temp_hits, length(temp)) for (w in 1:( length(temp_hits)-1) ) { temp_sub <- temp[temp_hits[w]:temp_hits[w+1]] temp_cons <- (grep("Consensus pattern", temp_sub)+1 ):(temp_sub %>% length -2) %>% temp_sub[.] %>% paste0(.,collapse="") temp_info <- temp_sub[grep("Indices:", temp_sub):(grep("Indices:", temp_sub)+1)] %>% paste0(.,collapse=" ") %>% gsub(" {2,}", " ", .) %>% strsplit(.,"[A-Z]") %>% .[[1]] %>% .[-1] %>% gsub(".*: ","",.) %>% strsplit(.,"--") %>% unlist %>% gsub(" {1,}","",.) %>% as.numeric temp_info <- temp_info[-length(temp_info)] temp_info <- temp_info %>% c(.,temp_cons) repeat_info <- rbind.data.frame(repeat_info, c(umis[i] %>% gsub(".*__|bc11|_rc","",.), umis[i] %>% gsub("__.*|bc10","",.), umis[i], temp[grep("Sequence", temp)] %>% gsub(".* ","",.),temp_info)) } } } repeat_info <- repeat_info[-1,] write.csv(repeat_info ,"repeat_info.csv") repeat_info <- read.csv("repeat_info.csv") repeat_info <- cbind.data.frame(repeat_info, "aln_start"= NA) repeat_info <- cbind.data.frame(repeat_info, "aln_stop"= NA) #remove dinuc repeats kickme <- repeat_info[( as.numeric(repeat_info[,"period"])<=2 & as.numeric(repeat_info[,"copies"]) > 20) ,"seq"] for (i in 1:length(umis)) { temp <- repeat_info[repeat_info[,"coords"]==umis[i],] if (!file.exists(paste0("./smp2/", umis[i], "_alnfilt.fasta"))) next if (nrow(temp)==0) next temp_aln <- readDNAStringSet(paste0("./smp2/", umis[i], "_alnfilt.fasta")) temp_aln <- temp_aln[!temp_aln %>% names %in% kickme] #map base pairs to unaligned sequences temp_base_coords <- temp_aln %>% gsub("[A-Z]","1",.) %>% gsub("-","0",.) %>% strsplit(.,"") %>% lapply(.,"as.numeric") %>% lapply(.,cumsum) for (j in 1:nrow(temp)) { temp_seq <- temp_base_coords[names(temp_aln) %in% temp[j,"seq"]] if (length(temp_seq)==0) next repeat_info[repeat_info[,"seq"]==temp[j,"seq"],"aln_start"] <- which((temp_seq %>% unlist) == ( temp[j,"start"] %>% as.numeric))[1] repeat_info[repeat_info[,"seq"]==temp[j,"seq"],"aln_stop"] <- which((temp_seq %>% unlist) == ( temp[j,"stop"] %>% as.numeric))[1] } } # # plot(round(as.numeric(repeat_info[,"period"])/10)) # abline(h=16, col="red") # abline(h=26, col="red") #i<- which(umis=="bc1065__bc1103_rc") repeat_info <- cbind.data.frame(repeat_info, "actual_start"= NA) repeat_info <- cbind.data.frame(repeat_info, "actual_stop"= NA) #Algae i <- 78 #truncate alignments to delete repeats for (i in 1:length(umis)) { if (!file.exists(paste0("./smp2/", umis[i], "_alnfilt.fasta"))) next temp_aln <- readDNAStringSet(paste0("./smp2/", umis[i], "_alnfilt.fasta")) temp_aln <- temp_aln[!temp_aln %>% names %in% kickme] temp <- repeat_info[repeat_info[,"coords"]==umis[i],] temp <- temp[complete.cases(temp[,1:12]),] temp <- temp[as.numeric(temp[,"period"]) > 160 & as.numeric(temp[,"period"]) < 260,] if (nrow(temp)<1) { writeLines(text = umis[i], con = paste0("./smp2/", umis[i], "_NOREPEATS.fasta") ) writeXStringSet(temp_aln, filepath= paste0("./smp2/", umis[i], "_alncut.fasta")) next } #keep only strongest per sequence temp <- temp[order(temp[,"plate"]%>%as.numeric,temp[,"sample"]%>%as.numeric,temp[,"coords"],temp[,"seq"], temp[,"score"]%>%as.numeric%>%rev),] temp <- temp[temp[,"seq"] %>% duplicated %>% not,] #keep only best repeat ids if ( nrow(temp)>10 & sum(as.numeric(temp[,"score"])>1000)>10 ) temp <- temp[as.numeric(temp[,"score"])>1000,] if (nrow(temp)==0) { writeLines(text = umis[i], con = paste0("./smp2/", umis[i], "_NOREPEATS.fasta") ) writeXStringSet(temp_aln, filepath= paste0("./smp2/", umis[i], "_alncut.fasta")) next } #remove outliers temp_upper_IQR <- temp[,"aln_stop"] %>% IQR(.,na.rm=T) if (temp_upper_IQR==0) temp_upper_IQR <- 1 temp_upper_median <- temp[,"aln_stop"] %>% median(.,na.rm=T) temp <- temp[temp[,"aln_stop"] < (temp_upper_median + temp_upper_IQR * 2) & temp[,"aln_stop"] > (temp_upper_median - temp_upper_IQR * 2),] temp_lower_IQR <- temp[,"aln_start"] %>% IQR(.,na.rm=T) if (temp_lower_IQR==0) temp_lower_IQR <- 1 temp_lower_median <- temp[,"aln_start"] %>% median(.,na.rm=T) temp <- temp[temp[,"aln_start"] < (temp_lower_median + temp_lower_IQR * 2) & temp[,"aln_start"] > (temp_lower_median - temp_lower_IQR * 2),] if (nrow(temp) < 3 ) { writeLines(text = umis[i], con = paste0("./smp2/", umis[i], "_NOREPEATS.fasta") ) writeXStringSet(temp_aln, filepath= paste0("./smp2/", umis[i], "_alncut.fasta")) next } repeat_info[repeat_info[,"coords"]==umis[i],"actual_start"] <- temp[,"aln_start"] %>% min repeat_info[repeat_info[,"coords"]==umis[i],"actual_stop"] <- temp[,"aln_stop"] %>% max temp_cut <- xscat( subseq(temp_aln, start=1, end= temp[,"aln_start"] %>% min ), subseq(temp_aln, start=temp[,"aln_stop"] %>% max, end= width(temp_aln) ) ) temp_repeats <- subseq(temp_aln, start=temp[,"aln_start"] %>% min, end= temp[,"aln_stop"] %>% max ) names(temp_cut) <- names(temp_aln) names(temp_repeats) <- names(temp_aln) writeXStringSet(temp_cut, filepath= paste0("./smp2/", umis[i], "_alncut.fasta")) writeXStringSet(temp_repeats, filepath= paste0("./smp2/", umis[i], "_removed_repeats.fasta")) } #analyse repeats # temp <- repeat_info # temp <- temp[complete.cases(temp[,1:12]),] # temp <- temp[as.numeric(temp[,"period"]) > 160 & as.numeric(temp[,"period"]) < 260,] # temp2 <- temp[,"cons"] %>% DNAStringSet # names(temp2) <- paste0( temp[,"coords"], "___", temp[,"seq"]) # writeXStringSet(temp2, "all_repeats.fasta") # # system(command = paste0(vsearch_exec, " --cluster_size ", # paste0("all_repeats.fasta"), # " --id 0.97 --threads ",number_of_cores, # " --consout ", # paste0("all_repeats_cons.fasta")," --uc ", # paste0("all_repeats_cluster.csv"), # " --iddef 1 --clusterout_id" # )) # # for (i in 1:length(umis)) # { # if (!file.exists(paste0("./smp2/", umis[i], "_alncut.fasta"))) next # # x <- temp_hits <- readDNAStringSet( paste0("./smp2/", umis[i], "_alncut.fasta") ) # # keepme_thresh <- 0.05 # # coords <- x %>% alphabetByCycle # # for (k in 1:ncol(coords)) # { # temp_max <- which.max(coords[,k]) %>% coords[.,k] # temp_dom <- which.max(coords[,k]) %>% coords[.,k] %>% names # coords[temp_dom,k] <- coords[temp_dom,k] + coords["N",k] # coords["N",k] <- 0 # coords[coords[,k] / temp_max < keepme_thresh ,k] <- 0 # } # # for (w in 1:length(x)) # { # temp_seq <- x[w,] %>% as.character %>% strsplit(.,"") %>% unlist # for (j in 1:length(temp_seq)) # { # temp_dom <- which.max(coords[,j]) %>% coords[.,j] %>% names # temp_should_be <- ( ( coords[,j] >= keepme_thresh ) %>% coords[.,j] %>% names ) # if (temp_seq[j] %in% temp_should_be) { next } else { temp_seq[j] <- temp_dom } # } # x[w,] <- temp_seq %>% paste0(.,collapse="") %>% DNAStringSet # } # # writeXStringSet(x, file="temp.fasta") # # temp <- ape::read.dna("temp.fasta", format="fasta") # temp <- ape::del.colgapsonly(temp, threshold = 0.9) # ape::write.dna(temp, file= "temp.fasta", format="fasta", nbcol=10^6, colsep="", colw=10^6, blocksep=0) # # temp_hits <- readDNAStringSet( "temp.fasta" ) # # #mask homopolymeric regions # temp_homo <- matrix(nrow=0, ncol=0) # for (j in c("AAAAAAA", "CCCCCCC", "GGGGGGG", "TTTTTTT")) # temp_homo <- rbind.data.frame(vmatchPattern(j, temp_hits, max.mismatch=0, fixed=T, with.indels = T) %>% as.data.frame, temp_homo) # # for (j in c("AAAAAAAAA", "CCCCCCCCC", "GGGGGGGGG", "TTTTTTTTT")) # temp_homo <- rbind.data.frame(vmatchPattern(j, temp_hits, max.mismatch=2, fixed=T) %>% as.data.frame, temp_homo) # # # mask_pos <- NULL # for (j in 1:nrow(temp_homo)) # mask_pos <- c((temp_homo[j,"start"]-2):(temp_homo[j,"end"]+2), mask_pos) # # mask_pos %<>% unique %>% sort # # mask_pos <- mask_pos[mask_pos < width(temp_hits[1]) & mask_pos > 1] # # for (j in 1:length(temp_hits)) # temp_hits[j] <- replaceLetterAt(x = temp_hits[j] %>% unlist, at = mask_pos , letter = rep("N", length(mask_pos))) %>% DNAStringSet # # writeXStringSet(temp_hits, paste0("./smp2/", umis[i], "_nohomo.fasta")) # } # # # for (i in 1:length(umis)) # { # if (!file.exists(paste0("./smp2/", umis[i], "_nohomo.fasta"))) next # # temp_hits <- readDNAStringSet(paste0("./smp2/", umis[i], "_nohomo.fasta")) # # kmer_counts <- oligonucleotideFrequency(temp_hits, with.labels = TRUE, step = 1 , width = 6, simplify.as="matrix") # # #DNAbin conversion for comparison # # temp_label <- temp_hits@id # # temp <- temp_hits # # temp <- temp@sread %>% as.character %>% sapply(.,strsplit, split="") %>% as.DNAbin # # names(temp) <- temp_label # # kmer_counts <- kcount(temp, k=6) # # umres <- umap(kmer_counts, n_neighbors = min(temp_hits %>% length, 10), n_components=2, min.dist=0) # # temp_points <- temp_hits %>% length # if (temp_hits %>% length >= 6) temp_points <- 6 # if (temp_hits %>% length > 40) temp_points <- 10 # if (temp_hits %>% length > 100) temp_points <- round(temp_hits %>% length/10)-1 # temp_clust <- hdbscan(umres$layout, minPts = temp_points) # # if ( temp_clust$membership_prob %>% mean < 0.1 & temp_hits %>% length > 20) # temp_clust <- hdbscan(umres$layout, minPts = temp_points * 2) # # # png(filename= paste0("./smp2/",umis[i],"_cluster_outlier_img.png"), width=1560, height=1560, res = 300) # #create image of clustering, all black = no clusters # colors <- mapply(function(col, i) adjustcolor(col, alpha.f = temp_clust$outlier_scores[i]), # palette()[temp_clust$cluster+1], seq_along(temp_clust$cluster)) # plot(umres$layout, col=colors, pch=20) # dev.off() # # # png(filename= paste0("./smp2/",umis[i],"_cluster_membership_img.png"), width=1560, height=1560, res = 300) # #create image of clustering, all black = no clusters # colors <- mapply(function(col, i) adjustcolor(col, alpha.f = temp_clust$membership_prob[i]), # palette()[temp_clust$cluster+1], seq_along(temp_clust$cluster)) # plot(umres$layout, col=colors, pch=20) # dev.off() # # # temp_hits@id <- temp_hits@id %>% paste0(temp_clust$cluster,"_",.) %>% BStringSet() # # writeFasta(temp_hits, "temp.fasta") # # #if no cluster, define all points as cluster 1 # if (temp_clust$cluster %>% unique %>% length == 1 & temp_clust$cluster[1] == 0) # temp_clust$cluster <- rep(1, temp_clust$cluster %>% length) # # #if cluster size > 5 then apply membership filter of 0.2 # temp_assign <- temp_clust$cluster # temp_clus_ids <- temp_assign %>% unique # temp_clus_number <- temp_clus_ids %>% length # # for (j in 1:temp_clus_number) # { # temp_coords <- temp_clust$cluster == temp_clus_ids[j] # if ( sum(temp_coords) < 6 ) next # if ( sum(temp_clust$membership_prob[temp_coords] >= 0.2) < 6) next # temp_clust$cluster[temp_coords][temp_clust$membership_prob[temp_coords] < 0.2] <- -1 # } # # if (!file.exists(paste0("./smp2/", umis[i], "_alncut.fasta"))) next # x <- readDNAStringSet( paste0("./smp2/", umis[i], "_alncut.fasta") ) # # keepme_thresh <- 0.1 # # coords <- x %>% alphabetByCycle # # for (k in 1:ncol(coords)) # { # temp_max <- which.max(coords[,k]) %>% coords[.,k] # temp_dom <- which.max(coords[,k]) %>% coords[.,k] %>% names # coords[temp_dom,k] <- coords[temp_dom,k] + coords["N",k] # coords["N",k] <- 0 # coords[coords[,k] / temp_max < keepme_thresh ,k] <- 0 # } # # for (w in 1:length(x)) # { # temp_seq <- x[w,] %>% as.character %>% strsplit(.,"") %>% unlist # for (j in 1:length(temp_seq)) # { # temp_dom <- which.max(coords[,j]) %>% coords[.,j] %>% names # temp_should_be <- ( ( coords[,j] >= keepme_thresh ) %>% coords[.,j] %>% names ) # if (temp_seq[j] %in% temp_should_be) { next } else { temp_seq[j] <- temp_dom } # } # x[w,] <- temp_seq %>% paste0(.,collapse="") %>% DNAStringSet # } # # writeXStringSet(x, file="temp2.fasta") # # temp <- ape::read.dna("temp2.fasta", format="fasta") # temp <- ape::del.colgapsonly(temp, threshold = 0.5) # ape::write.dna(temp, file= paste0("./smp2/", umis[i], "_refined_cleaned.fasta"), format="fasta", nbcol=10^6, colsep="", colw=10^6, blocksep=0) # # # # if (!file.exists(paste0("./smp2/", umis[i], "_refined_cleaned.fasta"))) next # temp <- readDNAStringSet(paste0("./smp2/", umis[i], "_refined_cleaned.fasta")) # # temp_N <- temp %>% gsub("-","N",.) %>% DNAStringSet # # #get number of clusters # name_clus <- temp_clust$cluster %>% unique # # for (j in 1:length(name_clus)) # { # if (name_clus[j] < 1) next # temp_subset <- temp_N[temp_clust$cluster == name_clus[j]] # # IUPAC_CODE_MAP[16:31] <- c("A-","C-","G-","T-","AC-","AG-","AT-","CG-","CT-","GT-","ACG-","ACT-","AGT-","CGT-","ACGT-","-") # attributes(IUPAC_CODE_MAP)$names[16:31] <- c("A","C","G","T","M","R","W","S","Y","K","V","H","D","B","N","-") # # if (temp_subset %>% length < 5) next # cons_thresh <- 0.2 # #if (temp_subset %>% length > 50) cons_thresh <- 0.1 # # cons <- consensusString(temp_subset, threshold=cons_thresh, ambiguityMap=IUPAC_CODE_MAP) # cons %<>% DNAStringSet # # cons %<>% gsub("N","",.) %>% DNAStringSet # # names(cons) <- paste0(umis[i],"_cluster",j,"_reads",temp_subset %>% length) # # writeFasta(temp_subset, file=paste0("./smp2/", umis[i],"_cluster",j,"_reads",temp_subset %>% length,".fasta")) # writeFasta(cons , file=paste0("./smp2/", umis[i],"_cluster",j,"_reads",temp_subset %>% length,"_cons.fasta")) # } # # # } # # # # # # # # # # # # # # # # # #