{"id":372,"date":"2022-03-29T22:29:46","date_gmt":"2022-03-29T14:29:46","guid":{"rendered":"https:\/\/linguopeng.top\/?p=372"},"modified":"2022-03-29T22:29:49","modified_gmt":"2022-03-29T14:29:49","slug":"predict_puls","status":"publish","type":"post","link":"https:\/\/linguopeng.top\/?p=372","title":{"rendered":"predict_puls"},"content":{"rendered":"\n
#!\/usr\/bin\/env Rscript\n\n#################################\n#\n# Functions\n#\n#################################\n\nprocess_pairs <- function(idx, td) {\n\n        mindbc <- dbcan[dbcan$protein_id %in% td$protein_id,]\n\n\tmindbc <- unique(mindbc[,c(\"protein_id\",\"hmm\")])\n\t\n\tgh <- data.frame(protein_id=\"\",hmm=\"\", stringsAsFactors=FALSE)\n\tif (nrow(mindbc) > 0) {\n        \tgh <- aggregate(mindbc$hmm, by=list(protein_id=mindbc$protein_id), paste, collapse=\";\")\n        }\n\tcolnames(gh) <- c(\"protein_id\",\"hmm\")\n\n        cmin = idx -1\n        while(cmin > 0) {\n\n                contains = 0\n                for (j in cmin:(cmin-4)) {\n\n                        if (j>0) {\n\n\t\t\t\tdist <- td$start[j+1] - td$end[j]\n\t\t\t\tif (dist > 500) {\n\t                                break\n        \t                }\n\n                                # check if it is susC\/D\n                                if (td$sus[j] == \"susC\" || td$sus[j] == \"susD\") {\n                                        contains = 1\n                                }\n\n                                # check if it is a GH\n                                prot_id = td$protein_id[j]\n                                ghvec <- gh$hmm[gh$protein_id==prot_id]\n                                if (length(ghvec) > 0) {\n                                        contains = 1\n                                }\n\n                        }\n\n                }\n\n                if (contains > 0) {\n                        cmin <- cmin -1\n                } else {\n                        break\n                }\n        }\n\n        cmin <- cmin + 1\n\n        cmax = idx + 1\n        while(1) {\n\n                contains = 0\n                for (j in cmax:(cmax+4)) {\n\n                        if (j <= nrow(td)) {\n\n\t\t\t\tdist <- td$start[j] - td$end[j-1]\n\t\t\t\tif (dist > 500) {\n\t\t\t\t\tbreak\n\t\t\t\t}\n\n                                # check if it is susC\/D\n                                if (td$sus[j] == \"susC\" || td$sus[j] == \"susD\") {\n                                        contains = 1\n                                }\n\n                                # check if it is a GH\n                                prot_id = td$protein_id[j]\n                                ghvec <- gh$hmm[gh$protein_id==prot_id]\n                                if (length(ghvec) > 0) {\n                                        contains = 1\n                                }\n                        }\n\n                }\n\n                if (contains > 0) {\n                        cmax <- cmax + 1\n                } else {\n                        break\n                }\n        }\n\n        cmax <- cmax - 1\n\n\tif (cmax > cmin) {\n\n        \tpultdf <- merge(td[cmin:cmax,], gh, all.x=TRUE, sort=FALSE)\n        \tpultdf <- pultdf[order(pultdf$order),]\n        \tpultdf[is.na(pultdf)] <- \"\"\n        \tpultdf$pulid <- rep(paste(\"PUL\",PULCOUNTER, sep=\"\"), nrow(pultdf))\n\n\t\tprint(paste(\"created PUL\", PULCOUNTER))\n\t\tprint(paste(cmin,\"-\",cmax))\n\n        \tALLPULS <<- rbind(ALLPULS, pultdf)\n\n        \tPULCOUNTER <<- PULCOUNTER + 1\n\t}\n        return(cmax)\n\n}\n\n#################################\n#\n# code\n#\n#################################\n\n# get command line arguments as an array\n# args <- commandArgs(trailingOnly = TRUE)\nid<-c(\"MAG609\")\nargs <-c(paste0(\"..\/feature_table\/\",id[1],\"_ft.txt\"),\n         paste0(\"..\/pfam\/\",id[1],\".pfam\"),\n         paste0(\"..\/dbcan\/\",id[1],\".out.dm.ps.filtered\"),\n         paste0(\"..\/puls\/\",id[1],\".puls.tsv\"),\n         paste0(\"..\/puls\/\",id[1],\".puls.sum.tsv\"),\n         id[1])\n\n# arguments\nft.name <- args[1]\npf.name <- args[2]\ndb.name <- args[3]\n\n# file base\nfile.out <- args[4]\nfile.sum <- args[5]\n\n# genome\ngenome <- args[6]\n\n# ft.name<-read.table(\"Prevotella_rodentium_ft.txt\",stringsAsFactors=FALSE)\n# colnames(ft.name)<-c(\"contig\",\"start\",\"end\",\"strand\",\"protein_id\",\"protein_name\")\ncn <- c(\"contig\",\"start\",\"end\",\"strand\",\"protein_id\",\"protein_name\")\nft <- read.table(ft.name, \n\t\theader=FALSE, \n\t\tsep=\"\\t\", \n\t\tstringsAsFactors=FALSE, \n\t\tcol.names=cn)\n\nft$order = 1:nrow(ft)\n\ncn <- c(\"protein_id\",\"as\",\"ae\",\"es\",\"ee\",\"ha\",\"hn\",\"type\",\"hs\",\"he\",\"hl\",\"bit\",\"e\",\"sig\",\"clan\",\"active\")\npfam <- read.table(pf.name, \n\t\tskip=29, header=FALSE, \n\t\tstringsAsFactors=FALSE, \n\t\tfill=TRUE, \n\t\tcol.names=cn)\n\ncn <- c(\"hmm\",\"hit_len\",\"protein_id\",\"query_len\",\"evalue\",\"hit_start\",\"hit_end\",\"query_start\",\"query_end\",\"cov\")\ndbcan <- read.table(db.name, \n\t\tsep=\"\\t\", \n\t\theader=FALSE, \n\t\tstringsAsFactors=FALSE, \n\t\tcol.names=cn)\ndbcan$hmm <- gsub(\".hmm\",\"\",dbcan$hmm)\n\n# GLOBAL VARIABLES\nPULCOUNTER <- 1\nALLPULS <- data.frame(protein_id=character(),\n\t\t      contig=character(),\n\t\t      start=numeric(),\n\t\t      end=numeric(),\n\t\t      strand=character(),\n\t\t      protein_name=character(),\n\t\t      order=numeric(), \n\t\t      active=character(), \n\t\t      sus=character(), \n\t\t      hmm=character(), \n\t\t      pulid=character(),\n\t\t      stringsAsFactors=FALSE)\n\n\n# filter pfam\nhalign_prop = (pfam$he - pfam$hs + 1) \/ pfam$hl\npfam <- pfam[halign_prop>=0.6,]\n\n# find and annotate susC\/susD\npfam$sus <- rep(\"none\", nrow(pfam))\npfam$sus[grep(\"PF00593\",pfam$ha)] <- \"susC\"\npfam$sus[grep(\"PF13715\",pfam$ha)] <- \"susC\"\npfam$sus[grep(\"PF07715\",pfam$ha)] <- \"susC\"\npfam$sus[grep(\"PF07980\",pfam$ha)] <- \"susD\"\npfam$sus[grep(\"PF12741\",pfam$ha)] <- \"susD\"\npfam$sus[grep(\"PF12771\",pfam$ha)] <- \"susD\"\npfam$sus[grep(\"PF14322\",pfam$ha)] <- \"susD\"\n\n# limit pfam to relevant rows\npfam <- pfam[pfam$sus != \"none\",]\n\n# merge with feature table\nftp <- merge(ft, pfam, by.x=\"protein_id\", by.y=\"protein_id\", all.x=TRUE, sort=FALSE)[,c(1:7,22,23)]\nftp <- unique(ftp)\nftp <- ftp[order(ftp$order),]\n\n# unique list of contigs\ncons <- unique(ftp$contig)\n\n# go through each contig\nfor (c in cons[order(cons)]) {\n\n\ttdf <- ftp[ftp$contig==c,]\n\ttdf <- tdf\n\ttdf[is.na(tdf)] <- \"\"\n\n\ti <- 1\n\twhile(i <= nrow(tdf)) {\n\t#for (i in 1:nrow(tdf)) {\n\n\t\t#print(c)\n\t\t#print(nrow(tdf))\n\t\t#print(i)\n\t\tif (i>=nrow(tdf)) {\n\t\t\tbreak\n\t\t}\n\t\tif (tdf$sus[i] == \"susC\" && tdf$sus[i+1] == \"susD\") {\n\t\t\t# we have a pair - do something\n\t\t\tprint(paste(\"pair at\",c,i))\n\t\t\ti <- process_pairs(i,tdf)\n\t\t\tprint(i)\n\t\t}\n\n\t\tif (i>=nrow(tdf)) {\n                        break\n                } \n\n\t\tif (tdf$sus[i] == \"susD\" && tdf$sus[i+1] == \"susC\") {\n\t\t\t# we have a pair - do something\n\t\t\tprint(paste(\"pair at\",c, i))\n\t\t\ti <- process_pairs(i,tdf)\n\t\t\tprint(i)\n\t\t}\n\t\t\n\t\ti <- i+1\n\t}\n\t\n}\n\nif (nrow(ALLPULS) >= 2) {\n\ttwo <- ALLPULS$end[1:(nrow(ALLPULS)-1)]\n\tone <- ALLPULS$start[2:nrow(ALLPULS)]\n\n\tALLPULS$dist <- c(0,one - two)\n\n\tALLPULS$genome <- rep(genome, nrow(ALLPULS))\n\n\tALLPULS <- ALLPULS[, c(\"genome\",\"pulid\",\"protein_id\",\"contig\",\"start\",\"end\",\"strand\",\"dist\",\"protein_name\",\"sus\",\"hmm\",\"active\")]\n\t\n\tALLPULS <- ALLPULS[order(ALLPULS$contig, ALLPULS$start),]\n\n\twrite.table(ALLPULS, file.out, col.names=TRUE, row.names=FALSE, quote=FALSE, sep=\"\\t\")\n\n} else {\n\tALLPULS$dist <- numeric()\n\tfile.create(file.out)\n}\n\nif (nrow(ALLPULS) >= 2) {\n\n\tALLPULS$GENE <- ALLPULS$hmm\n\tALLPULS$GENE[ALLPULS$sus!=\"\"] <- ALLPULS$sus[ALLPULS$sus!=\"\"]\n\tALLPULS$GENE[ALLPULS$GENE==\"\"] <- \"unk\"\n\n\tagg <- aggregate(ALLPULS$GENE, by=list(pulid=ALLPULS$pulid), paste, sep=\"-\", collapse=\"-\")\n\tuni <- unique(ALLPULS[,c(\"genome\",\"pulid\",\"contig\")])\n\tstart <- aggregate(ALLPULS$start,  by=list(pulid=ALLPULS$pulid), function(x) return(x[1]))\n\tend   <- aggregate(ALLPULS$end,  by=list(pulid=ALLPULS$pulid), function(x) return(x[length(x)]))\n\n\tpos <- merge(start, end, by=\"pulid\")\n\tdeets <- merge(uni, pos, by=\"pulid\")\n\n\tout <- merge(deets,agg,by=\"pulid\")\n\n\tcolnames(out) <- c(\"pulid\",\"genome\",\"contigid\",\"start\",\"end\",\"pattern\")\n\n\tout <- out[,c(\"genome\",\"pulid\",\"contigid\",\"start\",\"end\",\"pattern\")]\n\n\tout <- out[order(out$contigid, out$start),]\n\n\twrite.table(out, file.sum, col.names=TRUE, row.names=FALSE, quote=FALSE, sep=\"\\t\")\n} else {\n\tfile.create(file.sum)\n}\n\n\n<\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-372","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"_links":{"self":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/372","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=372"}],"version-history":[{"count":1,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/372\/revisions"}],"predecessor-version":[{"id":373,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/372\/revisions\/373"}],"wp:attachment":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=372"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=372"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=372"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}