PageRenderTime 26ms CodeModel.GetById 27ms RepoModel.GetById 1ms app.codeStats 0ms

/data/pcc6803/cyanobase.R

https://gitlab.com/raim/genomeBrowser
R | 324 lines | 210 code | 46 blank | 68 comment | 7 complexity | 54eb395bb73c7998701b50015f5baf25 MD5 | raw file
    

  1. 

  2. require(segmenTools)
    3. 

  3. require(readxl)
    4. 

  4. 

  5. ## load functions
    6. 

  6. browser.path <- sub("GENBRO=","",system("env|grep GENBRO",intern=TRUE))
    7. 

  7. data.path <- sub("GENDAT=","",system("env|grep GENDAT",intern=TRUE))
    8. 

  8. source(file.path(browser.path,"src/genomeBrowser_utils.R"))
    9. 

  9. ## data path
    10. 

  10. cyano.path <- sub("PCC6803DAT=","",system("env|grep PCC6803DAT",intern=TRUE))
    11. 

  11. 

  12. ## working with IDs, set this to false!
    13. 

  13. op <- options(stringsAsFactors=FALSE)
    14. 

  14. 

  15. 

  16. ### CYANOBASE ANNOTATION DATA
    17. 

  17. ## NOTE: this assumes equal sequence lengths
    18. 

  18. ## between cyanobase annotation and ncbi releases!
    19. 

  19. ## checked for cyanobase 20170923 vs. .1 versions of NCBI
    20. 

  20. 

  21. ## TODO behle21 and kopf14 data are parsed here AFTER this
    22. 

  22. ## file was required
    23. 

  23. 

  24. ## chromosome Length Index
    25. 

  25. chrIdx <- read.delim(file.path(cyano.path,"chromosomes","pcc6803.csv"),
    26. 

  26.                      header=TRUE,row.names=1)
    27. 

  27. chrL <- chrIdx[,2]
    28. 

  28. chrS <- segmenTools::getChrSum(chrL) # chromosome index
    29. 

  29. ## mapping of different chromosome IDs
    30. 

  30. chrMap <- read.delim(file.path(browser.path,"data","pcc6803","parameters",
    31. 

  31.                                "chromosomemapping.txt"),
    32. 

  32.                      header=FALSE,row.names=1)
    33. 

  33. 

  34. ## parse basic gene annotation
    35. 

  35. file <- file.path(cyano.path,"originalData", paste("genes.txt",sep=""))
    36. 

  36. cyan <- read.delim(file,header=FALSE)
    37. 

  37. colnames(cyan) <- c("assembly","ID","chr","start","end","strand","function")
    38. 

  38. cyan[,"chr"] <- as.numeric(chrMap[cyan[,"chr"],1])
    39. 

  39. ## correct negative coordinates of circular genes
    40. 

  40. ## USING NCBI chromosome lengths!!
    41. 

  41. circ <- which(cyan[,"start"]<0)
    42. 

  42. cyan[circ,"start"] <- chrL[cyan[circ,"chr"]]+cyan[circ,"start"]
    43. 

  43. ## rm unused columns
    44. 

  44. cyan <- cyan[,-1]
    45. 

  45. 

  46. ## NAMES
    47. 

  47. 

  48. ## use cyanobase GENE SYMBOLS as names/alias
    49. 

  49. genes <- read.delim(file.path(cyano.path,"originalData","gene_symbol.txt"),
    50. 

  50.                     header=FALSE)
    51. 

  51. cyan <- cbind.data.frame(cyan,
    52. 

  52.                          alias=genes[match(cyan[,"ID"],genes[,"V2"]),"V4"])
    53. 

  53. nms <- sub("/.*","",sub(",.*","",cyan[,"alias"]))
    54. 

  54. cyan <- cbind.data.frame(cyan, name=nms)
    55. 

  55. 

  56. ## use IDs where no name is present
    57. 

  57. cyan[is.na(cyan[,"name"]),"name"] <- cyan[is.na(cyan[,"name"]),"ID"]
    58. 

  58. 

  59. ## TODO: use first of ,- or /-separated lists
    60. 

  60. ## store full names in "alias" or similar
    61. 

  61. ## in annotation.RData - add kopf14 transcripts and cyanobase, but
    62. 

  62. ## only for plotFeatures.R reference (rm redundant from plot types)
    63. 

  63. 

  64. 

  65. ## add annotation data?
    66. 

  66. ## record and write index file!
    67. 

  67. cat <- read.delim(file.path(cyano.path,"originalData","category.txt"),
    68. 

  68.                   header=FALSE)
    69. 

  69. cat <- apply(cat, 2, function(x) gsub("&amp;","&",gsub("&#x27;","'",x)))
    70. 

  70. 

  71. catstr <- unlist(gsub(";$","",apply(cat[match(cyan[,"ID"],cat[,"V1"]),
    72. 

  72.                                         c("V3","V5")],1,
    73. 

  73.                                     paste, collapse=";")))
    74. 

  74. cyan <- cbind.data.frame(cyan, category=catstr)
    75. 

  75. 

  76. ## ADD UNIPROT GO
    77. 

  77. go.file <- file.path(cyano.path,"originalData","pcc6803_uniprot_go.txt")
    78. 

  78. go <- read.delim(go.file)
    79. 

  79. colnames(go) <- c("uniprot","cyanobaseID","go")
    80. 

  80. go$cyanobaseID <- gsub("syn:","",go$cyanobaseID)
    81. 

  81. ## extract all GO terms & descriptions
    82. 

  82. got <- strsplit(go$go,";")
    83. 

  83. god <- got <- unique(unlist(got))
    84. 

  84. got <- sapply(got, function(x) gsub("\\]","",gsub(".*\\[","",x)))
    85. 

  85. god <- sapply(god, function(x) trimws(gsub("\\[.*","",x)))
    86. 

  86. names(god) <- names(got) <- NULL
    87. 

  87. god <- cbind(term=got, description=god)
    88. 

  88. god <- god[!duplicated(god[,1]),]
    89. 

  89. ## WRITE OUT GO MAPPING
    90. 

  90. go.map <- file.path(cyano.path,"processedData","pcc6803_uniprot_go_terms.tsv")
    91. 

  91. write.table(god, file=go.map, quote=FALSE, row.names=FALSE, sep="\t")
    92. 

  92. 

  93. ## generate pure GO term column
    94. 

  94. got <- strsplit(go$go,";")
    95. 

  95. got <- lapply(got, function(x) gsub("\\]","",gsub(".*\\[","",x)))
    96. 

  96. got <- unlist(lapply(got, paste0, collapse=";"))
    97. 

  97. ## expand genes
    98. 

  98. ggenes <- strsplit(go$cyanobaseID,";")
    99. 

  99. glen <- unlist(lapply(ggenes, length))
    100. 

  100. ## loose two empty!!
    101. 

  101. rm <- which(glen==0)
    102. 

  102. go <- go[-rm,]
    103. 

  103. ggenes <- ggenes[-rm]
    104. 

  104. glen <- glen[-rm]
    105. 

  105. got <- got[-rm]
    106. 

  106. ## expand multiple cyanobase genes
    107. 

  107. tmp <- sapply(1:length(glen), function(x) rep(got[x], glen[x]))
    108. 

  108. ok <- sum(unlist(lapply(tmp,length)) != unlist(lapply(ggenes,length)))
    109. 

  109. if ( ok>0 ) stop("wrong length of GO gene expansion")
    110. 

  110. ## gene - go table!
    111. 

  111. allgo <- cbind.data.frame(ID=unlist(ggenes), GO=unlist(tmp))
    112. 

  112. 

  113. ##tmp <- merge(cyan, allgo, by="ID", all=TRUE)
    114. 

  114. 

  115. cyan <- cbind.data.frame(cyan, GO=allgo$GO[match(cyan$ID, allgo$ID)])
    116. 

  116. 

  117. ## load lehmann14 transcriptome data and add clusterings and colors!
    118. 

  118. leh14 <- read.delim(file.path(cyano.path,"originalData",
    119. 

  119.                               "supp_gku641_SupportingDataFile2.csv"))
    120. 

  120. leh14cols <- c(CL_prakash09_sc="supercoiling_transcriptome",
    121. 

  121.                CL_leh14="diurnal_transcriptome",
    122. 

  122.                CL_leh14_at2="AT2_periodicity")
    123. 

  123. leh14dat <- leh14[match(cyan[,"ID"],leh14[,"ID"]),leh14cols]
    124. 

  124. colnames(leh14dat) <- names(leh14cols)
    125. 

  125. cyan <- cbind(cyan, leh14dat)
    126. 

  126.   
    127. 

  127. ## TODO: load this from file!
    128. 

  128. diur.cls <- as.character(cyan[,"CL_leh14"])
    129. 

  129. cls.srt <- sort(as.numeric(unique(diur.cls)))
    130. 

  130. ## set cluster colors
    131. 

  131. cls.col <- rep("gray", length(cls.srt))
    132. 

  132. names(cls.col) <- cls.srt
    133. 

  133. ### main
    134. 

  134. cls.col[c("2")] <- "#8493CA" # pastel blue
    135. 

  135. cls.col[c("7")] <- "#0000FF"   # blue
    136. 

  136. cls.col[c("3")] <- "#00E6FF" # cyan"
    137. 

  137. cls.col[c("8")] <- "#FF0000"   # red
    138. 

  138. cls.col[c("6")] <- rgb(1,.8,0) # dark yellow
    139. 

  139. cls.col[c("1")] <- "#82CA9D" # pastel green
    140. 

  140. ## background
    141. 

  141. cls.col[c("5")] <- "#626262" # gray: foreign genes?? transposons!
    142. 

  142. cls.col[c("4")] <- rgb(0,0,.6) # dark blue
    143. 

  143. cls.col[c("9")] <- "#6ECFF6" # pastel cyan
    144. 

  144. cls.col[c("10")]<- "#F7977A" # pastel red
    145. 

  145. ## not present on array
    146. 

  146. cls.col[c("11")]<- "#000000" # black #"#E1E1E1" # light gray
    147. 

  147. 

  148. ftcols <- cls.col[as.character(cyan[,"CL_leh14"])]
    149. 

  149. ftcols[is.na(ftcols)] <- "#000000"
    150. 

  150. cyan <- cbind(cyan,CL_leh14_col=ftcols)
    151. 

  151. 

  152. ## load saha16 transcriptome clustering (segmenTier/segmenTools)
    153. 

  153. load(file.path(cyano.path,"saha16","saha16.RData"))
    154. 

  154. cset <- cdat$cset
    155. 

  155. k <- selected(cset)
    156. 

  156. cls <- cset$clusters[,k]
    157. 

  157. cls.col <-  cset$colors[[k]]
    158. 

  158. ftcols <- cls.col[as.character(cls[cyan[,"ID"]])]
    159. 

  159. cyan <- cbind(cyan,
    160. 

  160.               CL_saha16=as.character(cls[cyan[,"ID"]]),
    161. 

  161.               CL_saha16_col=ftcols)
    162. 

  162. 

  163. ## load Zavrel et al. 2019 proteome clustering
    164. 

  164. zav <- read.delim(file.path(cyano.path,"processedData","zavrel19_proteome.tsv"))
    165. 

  165. 

  166. ## search gene IDs
    167. 

  167. zlst <- strsplit(zav$kegg,";")
    168. 

  168. zids <- lapply(zlst, function(x) {
    169. 

  169.     y <- c(match(x, cyan$ID))
    170. 

  170.     y[!is.na(y)]})
    171. 

  171. znum <- unlist(lapply(zids, length))
    172. 

  172. 

  173. cat(paste(sum(znum==2), "proteins with multiple genes:\n",
    174. 

  174.           paste(zav[znum==2,"protein"], collapse=" "),"\n",
    175. 

  175.           paste(zav[znum==2,"gene"], collapse=" "),
    176. 

  176.           "\nonly counting first\n"))
    177. 

  177. zids <- unlist(lapply(zlst, function(x) x[1]))
    178. 

  178. 

  179. length(match(zids, cyan$ID))
    180. 

  180. 

  181. cyan <- cbind(cyan, zav[match(cyan$ID, zids),
    182. 

  182.                         c("CL_zavrel19","CL_zavrel19_super")])
    183. 

  183. 

  184. ## load Behle et al. 2021 topA-OX clusters and colors
    185. 

  185. ## TODO: REPLACE THIS BY SUPPLEMENT OF THE PAPER (add to download list!)
    186. 

  186. behle21.file <- "~/work/CoilHack/experiments/RNAseq/20201204_RM_topA/analysis/genome_0/20201204_RM_topA_results.tsv"
    187. 

  187. behle21.batch <- "~/work/CoilHack/experiments/RNAseq/20200602_WM_coilhack_endpoint/analysis_4/20200602_WM_coilhack_endpoint_results.tsv"
    188. 

  188. if ( file.exists(behle21.file) ) {
    189. 

  189.     b21 <- read.delim(behle21.file)
    190. 

  190.     b21 <- b21[match(cyan[,"ID"], b21$cyanobaseID),]
    191. 

  191.     b21b <- read.delim(behle21.batch)
    192. 

  192.     b21b <- b21b[match(cyan[,"ID"], b21b$cyanobaseID),]
    193. 

  193. 

  194.     cyan <- cbind(cyan, b21[,c("CL_behle21","CL_behle21_col")])
    195. 

  195. 

  196.     ## gray scale for log2 FC
    197. 

  197.     num2col <- function(x,n=100, limits, pal, colf=viridis::viridis){
    198. 

  198.         if ( missing(pal) ) pal <- colf(n)
    199. 

  199.         if ( missing(limits) ) limits <- range(x, na.rm=TRUE)
    200. 

  200.         pal[findInterval(x,seq(limits[1],limits[2],
    201. 

  201.                                length.out=length(pal)+1), all.inside=TRUE)]
    202. 

  202.     }
    203. 

  203.     cyan <- cbind(cyan, b21b[,grep("log2FoldChange",colnames(b21b))])
    204. 

  204.     topa.col <- num2col(cyan$TOPA_log2FoldChange, limits=c(-1,1))
    205. 

  205.     gyrb.col <- num2col(cyan$GyrB_log2FoldChange, limits=c(-2,2))
    206. 

  206.     gyra.col <- num2col(cyan$GyrA_log2FoldChange, limits=c(-2,2))
    207. 

  207.     cyan <- cbind(cyan, TOPA_color=topa.col,
    208. 

  208.                   GyrA_color=gyra.col, GyrB_color=gyrb.col)
    209. 

  209. }
    210. 

  210. 

  211. ## add kopf14 TU and add CL_kopf14, if it exists already,
    212. 

  212. ## TODO: generate kopf14 file in preprocessing step!
    213. 

  213. kopf14.file <-  file.path(cyano.path,"processedData","kopf14_tu_genes.tsv")
    214. 

  214. if ( file.exists(kopf14.file) ) {
    215. 

  215.     kpf <- read.delim(kopf14.file)
    216. 

  216.     kpf <- kpf[match(cyan[,"ID"], kpf$gene),]
    217. 

  217.     cyan <- cbind(cyan, TU_kopf14=kpf[,"TU.ID"],
    218. 

  218.                   CL_kopf14=kpf[,"Max.cond."])
    219. 

  219. }
    220. 

  220. 

  221. ## add Koskinen et al. 2016 transcriptome data
    222. 

  222. kosk16.file <- file.path(cyano.path,"originalData",
    223. 

  223.                          "MMI_13214_supp-0001-Dataset_S1.xlsx")
    224. 

  224. kosk16 <- readxl::read_xlsx(kosk16.file, sheet=1, skip=1)
    225. 

  225. kosk16 <- kosk16[match(cyan[,"ID"], kosk16$`Gene ID`),2:3]
    226. 

  226. colnames(kosk16) <- paste0("sigBCDE_",gsub(" ","",colnames(kosk16)))
    227. 

  227. cyan <- cbind(cyan, kosk16)
    228. 

  228. 

  229. ## FUNCTIONAL CATEGORIES
    230. 

  230. ## extract categories and store index!
    231. 

  231. maincat <- cbind(unique(cat[,3]), unique(cat[,4]))
    232. 

  232. tmp <- unique(cat[,5])
    233. 

  233. tmp <- tmp[tmp!=""]
    234. 

  234. subcat <- cbind(tmp,rep(NA,length(tmp)))
    235. 

  235. for ( i in 1:length(tmp) ) 
    236. 

  236.     subcat[i,2] <- c(cat[which(cat[,5]==tmp[i])[1],6])
    237. 

  237. catidx <- rbind(maincat,subcat)
    238. 

  238. catidx <- catidx[order(catidx[,1]),]
    239. 

  239. file.name <- file.path(cyano.path,"processedData","annotation_index.txt")
    240. 

  240. write.table(catidx, file=file.name, quote=F, col.names=F, row.names=F,sep="\t")
    241. 

  241. 

  242. ## GO TERMS - TODO
    243. 

  243. #go <- read.delim(file.path(cyano.path,"originalData","goterm.txt"),header=FALSE)
    244. 

  244. 

  245. ## add type column
    246. 

  246. cyan <- cbind(cyan,type=rep("gene",nrow(cyan)))
    247. 

  247. 

  248. ## reorder
    249. 

  249. mainCols <- c("ID","name","type","chr","start","end","strand")
    250. 

  250. cyan <- cbind(cyan[,mainCols],cyan[,!colnames(cyan)%in%mainCols])
    251. 

  251. 

  252. ## add chromosomes!
    253. 

  253. chr <- as.data.frame(matrix(NA, nrow=length(chrL), ncol=ncol(cyan)))
    254. 

  254. colnames(chr) <- colnames(cyan)
    255. 

  255. chr[,"chr"] <- 1:length(chrL)
    256. 

  256. chr[,"start"] <- 1
    257. 

  257. chr[,"end"] <- chrL
    258. 

  258. chr[,"strand"] <- "."
    259. 

  259. chr[,"ID"] <- 1:length(chrL)
    260. 

  260. chr[,"name"] <- as.character(chrIdx[,1])
    261. 

  261. chr[,"type"] <- "chromosome"
    262. 

  262. 

  263. cyan <- rbind(cyan,chr)
    264. 

  264. 

  265. ## * RENAME GENES!
    266. 

  266. ## gyrA/gyrA:
    267. 

  267. ## rename sll1941 to gyrA2 (48.61% AA identity to S.coccus PCC7942)
    268. 

  268. ## keep slr0417 as gyrA (62.4% AA identity to S.coccus PCC7942)
    269. 

  269. idx <- which(cyan[,"ID"]=="sll1941")
    270. 

  270. cyan[idx,"name"] <- "gyrA2"
    271. 

  271. 

  272. ## WRITE-OUT TABLE FILE
    273. 

  273. file.name <- file.path(cyano.path,"features_cyanobase.csv")
    274. 

  274. write.table(cyan, file=file.name, sep="\t", row.names=F, quote=F, na="")
    275. 

  275. 

  276. ## EXPAND CIRCULAR GENES for overlap analyses and genomeBrowser
    277. 

  277. ## keep color annotation!
    278. 

  278. cyan.circ <- segmenTools::expandCircularFeatures(cyan, chrL, copyCols=TRUE)
    279. 

  279. 

  280. 

  281. ## genomeBrowser Data!!
    282. 

  282. columns <- c(ID="ID", name="name", type="type", chr="chr", strand="strand",
    283. 

  283.              start="start", end="end", color="CL_behle21_col")
    284. 

  284. types <- c("gene","gene-circ2")
    285. 

  285. out.path <- file.path(data.path,"pcc6803")
    286. 

  286. odat <- list(ID="cyanobase",
    287. 

  287.              description="annotated genome features",
    288. 

  288.              data=cyan.circ,
    289. 

  289.              settings=list(type="plotFeatures",
    290. 

  290.                types=types,columns=columns,
    291. 

  291.                typord=FALSE,cuttypes=TRUE,names=TRUE,
    292. 

  292.                ylab="cyano-\nbase",hgt=.75,xaxis=FALSE))
    293. 

  293. ## store data
    294. 

  294. file.name <- file.path(out.path,paste(odat$ID,".RData",sep=""))
    295. 

  295. save(odat,file=file.name)
    296. 

  296. 

  297. ## features with color gradient by log2-fold change in gyrA-kd
    298. 

  298. ## features with color gradient by log2-fold change in topA-OX
    299. 

  299. ## copy each feature
    300. 

  300. 

  301. columns["color"] <- "behle21_batch_color"
    302. 

  302. cyx <- rbind(cyan.circ,
    303. 

  303.              cyan.circ,
    304. 

  304.              cyan.circ)
    305. 

  305. cyx <- cbind(cyx,
    306. 

  306.              "behle21_batch_color"=c(cyan.circ$TOPA_color,
    307. 

  307.                                      cyan.circ$GyrA_color,
    308. 

  308.                                      cyan.circ$GyrB_color))
    309. 

  309. cyx$type <- c(rep("TopA-OX",nrow(cyan.circ)),
    310. 

  310.               rep("GyrA-kd",nrow(cyan.circ)),
    311. 

  311.               rep("GyrB-kd",nrow(cyan.circ)))
    312. 

  312. types <- unique(cyx$type)
    313. 

  313. odat <- list(ID="behle21_batch",
    314. 

  314.              description="features by expression level in topoisomerase experiments, Behle et al. 2021",
    315. 

  315.              data=cyx,
    316. 

  316.              settings=list(type="plotFeatures",
    317. 

  317.                            arrow = list(length=0, code=0, angle=0,
    318. 

  318.                                         pch=NA, lwd=2),
    319. 

  319.                            types=types,columns=columns,
    320. 

  320.                            typord=TRUE,cuttypes=FALSE,names=FALSE,
    321. 

  321.                            ylab="",hgt=.25,xaxis=FALSE))
    322. 

  322. ## store data
    323. 

  323. file.name <- file.path(out.path,paste(odat$ID,".RData",sep=""))
    324. 

  324. save(odat,file=file.name)
    325. 

  325.