# TODO: Add comment # # Author: E.Korsching 2024 ############################################################################### # convert certain gff,gtf information to a Ensemble compatible gtf file # e.g. miRNA or tRNA # customize fasta for comparison reasons #fasta_keep2(fin="/home/korschi/1/mature.fa", fout="/home/korschi/1/hsa_mature.fa", tp1=">", tp2="hsa") # rows in 97770 / out 5312 # /2 = # 2656 ## miR ## # mirbase original gff important MI0022705 , hsa-mir-6859-1 # chr1 . miRNA_primary_transcript 17369 17436 . - . ID=MI0022705;Alias=MI0022705;Name=hsa-mir-6859-1 # chr1 . miRNA 17409 17431 . - . ID=MIMAT0027618;Alias=MIMAT0027618;Name=hsa-miR-6859-5p;Derives_from=MI0022705 # chr1 . miRNA 17369 17391 . - . ID=MIMAT0027619;Alias=MIMAT0027619;Name=hsa-miR-6859-3p;Derives_from=MI0022705 # # gtf translated from gff by gffread # chr1 . transcript 17369 17391 . - . transcript_id "MIMAT0027619"; gene_id "MIMAT0027619" # chr1 . exon 17369 17391 . - . transcript_id "MIMAT0027619"; # # Ensemble 111 miR format (but close to regular gene format) # miR #21 mirbase gene 25573980 25574044 . + . gene_id "ENSG00000283904"; gene_version "1"; gene_name "MIR155"; gene_source "mirbase"; gene_biotype "miRNA"; #21 mirbase transcript 25573980 25574044 . + . gene_id "ENSG00000283904"; gene_version "1"; transcript_id "ENST00000385060"; transcript_version "1"; gene_name "MIR155"; gene_source "mirbase"; gene_biotype "miRNA"; transcript_name "MIR155-201"; transcript_source "mirbase"; transcript_biotype "miRNA"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; #21 mirbase exon 25573980 25574044 . + . gene_id "ENSG00000283904"; gene_version "1"; transcript_id "ENST00000385060"; transcript_version "1"; exon_number "1"; gene_name "MIR155"; gene_source "mirbase"; gene_biotype "miRNA"; transcript_name "MIR155-201"; transcript_source "mirbase"; transcript_biotype "miRNA"; exon_id "ENSE00001500066"; exon_version "1"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; # # gene # 1 havana gene 182696 184174 . + . gene_id "ENSG00000279928"; gene_version "2"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; # 1 havana transcript 182696 184174 . + . gene_id "ENSG00000279928"; gene_version "2"; transcript_id "ENST00000624431"; transcript_version "2"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "DDX11L17-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; # 1 havana exon 182696 182746 . + . gene_id "ENSG00000279928"; gene_version "2"; transcript_id "ENST00000624431"; transcript_version "2"; exon_number "1"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "DDX11L17-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; exon_id "ENSE00003759020"; exon_version "2"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; # 1 havana exon 183132 183216 . + . gene_id "ENSG00000279928"; gene_version "2"; transcript_id "ENST00000624431"; transcript_version "2"; exon_number "2"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "DDX11L17-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; exon_id "ENSE00003759581"; exon_version "2"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; # 1 havana exon 183494 183571 . + . gene_id "ENSG00000279928"; gene_version "2"; transcript_id "ENST00000624431"; transcript_version "2"; exon_number "3"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "DDX11L17-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; exon_id "ENSE00003804405"; exon_version "1"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; # 1 havana exon 183740 183901 . + . gene_id "ENSG00000279928"; gene_version "2"; transcript_id "ENST00000624431"; transcript_version "2"; exon_number "4"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "DDX11L17-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; exon_id "ENSE00003807458"; exon_version "1"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; # 1 havana exon 183981 184174 . + . gene_id "ENSG00000279928"; gene_version "2"; transcript_id "ENST00000624431"; transcript_version "2"; exon_number "5"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "DDX11L17-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; exon_id "ENSE00003760199"; exon_version "2"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; # remove header first # eventually recreate: # #text -> #!text convert.gff.gtf.mir <- function(fin, fout, mfout, ENSGstart, ENSTstart){ # miRNA specific converter gff to gtf (for mirbase gff file) # fin, fout, mfout: gff path/file, gtf path/file, map path/file # ENSGstart, ENSTstart: number with 11 positions # see analytics in gtf_info.txt # e.g. at the moment: ENSG 298000 ENST 720000 arbitraily chosen - should work # because only the id links, name and accession/chr - increase/decrease of digit range will raise code problems (!) gff <- read.table(fin, header=F, sep="\t") # remove lines 'primary_transcript' logi1 <- grepl("primary_transcript", x=gff[,3]) gff <- gff[!logi1, ] gff.len <- nrow(gff) # customize second matrix block second <- data.frame(matrix("",gff.len,4)) for(i in 1:gff.len){ second[i,] <- unlist( strsplit(gff[i,9], split=";", fixed=T ) ) # char df } # shrink basic matrix now gff <- gff[ , -9] # remove duplicated in 'Name' dup <- duplicated(second[,3]) second <- second[!dup, ] gff <- gff[!dup, ] gff.len <- nrow(gff) # create empty Ensemble-column-9 data structure with dummy data new.df.len <- 3*gff.len new.df <- data.frame(rep("x",new.df.len), rep("x",new.df.len), rep("x",new.df.len), rep("x",new.df.len), rep("x",new.df.len), rep(".",new.df.len), rep("+-",new.df.len), rep(".",new.df.len), rep("x",new.df.len) ) # create additionally a final mapping table for the salmon analysis # ENST ENSG symbol map.df <- data.frame(ensembl_transcript_id=rep("x",gff.len), ensembl_gene_id=rep("x",gff.len), hgnc_symbol=rep("x",gff.len)) # fill k <- 1 for(i in 1:gff.len){ tmp1 <- sub('chr', '', gff[i, 1], fixed=T) # first row new.df[k, 1] <- tmp1 new.df[k, 2] <- "mirbase" new.df[k, 3] <- "gene" new.df[k, 4] <- gff[i, 4] new.df[k, 5] <- gff[i, 5] new.df[k, 6] <- "." new.df[k, 7] <- gff[i, 7] new.df[k, 8] <- "." tmp2 <- sub('Name=', '', second[i, 3], fixed=T) new.df[k, 9] <- paste('gene_id "ENSG00000',ENSGstart,'"; gene_name "',tmp2,'"; gene_source "mirbase"; gene_biotype "miRNA";', sep="") # mapping map.df[i, 1] <- paste("ENST00000",ENSTstart, sep="") map.df[i, 2] <- paste("ENSG00000",ENSGstart, sep="") map.df[i, 3] <- tmp2 k <- k+1 # second row new.df[k, 1] <- tmp1 new.df[k, 2] <- "mirbase" new.df[k, 3] <- "transcript" new.df[k, 4] <- gff[i, 4] new.df[k, 5] <- gff[i, 5] new.df[k, 6] <- "." new.df[k, 7] <- gff[i, 7] new.df[k, 8] <- "." # not the full Ensemble info new.df[k, 9] <- paste('gene_id "ENSG00000',ENSGstart,'"; transcript_id "ENST00000',ENSTstart,'"; gene_name "',tmp2,'"; gene_source "mirbase"; gene_biotype "miRNA";', sep="") k <- k+1 # third row new.df[k, 1] <- tmp1 new.df[k, 2] <- "mirbase" new.df[k, 3] <- "exon" new.df[k, 4] <- gff[i, 4] new.df[k, 5] <- gff[i, 5] new.df[k, 6] <- "." new.df[k, 7] <- gff[i, 7] new.df[k, 8] <- "." # not the full Ensemble info new.df[k, 9] <- paste('gene_id "ENSG00000',ENSGstart,'"; transcript_id "ENST00000',ENSTstart,'"; exon_number "1"; gene_name "',tmp2,'"; gene_source "mirbase"; gene_biotype "miRNA";', sep="") k <- k+1 ENSGstart <- ENSGstart+1 ENSTstart <- ENSTstart+1 } write.table(new.df, fout, append=F, quote=F, sep="\t", row.names=F, col.names=F) write.table(map.df, mfout, append=F, quote=F, sep="\t", row.names=F, col.names=T) cat("\n gff rows in ",gff.len," gtf rows out ",new.df.len) cat("\n ENSG final ",ENSGstart," ENST final ",ENSTstart," control and for next chunck of data...\n") return() } #convert.gff.gtf.mir(fin="/home/korschi/1/hsa_mir.gff3", # fout="/home/korschi/1/hsa_mir.gtf", # mfout="/home/korschi/1/hsa_mir_gtf_map.txt", # ENSGstart=298000, ENSTstart=720000) # gff rows in 2649 gtf rows out 7947 # ENSG final 300649 ENST final 722649 control and for next chunck of data... ## (!) add our 4 miRs manually at chr 16 ## ## Ensembl gtf to map file ## # Ensemble # gene # 1 havana gene 182696 184174 . + . gene_id "ENSG00000279928"; gene_version "2"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; # 1 havana transcript 182696 184174 . + . gene_id "ENSG00000279928"; gene_version "2"; transcript_id "ENST00000624431"; transcript_version "2"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "DDX11L17-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; # 1 havana exon 182696 182746 . + . gene_id "ENSG00000279928"; gene_version "2"; transcript_id "ENST00000624431"; transcript_version "2"; exon_number "1"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "DDX11L17-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; exon_id "ENSE00003759020"; exon_version "2"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; create.gtf.ens.map <- function(fin, mfout){ # Ensembl gtf specific creator of a map file ENST ENSG gene_name # fin, mfout: gtf path/file, map path/file gtf <- read.table(fin, header=F, sep="\t") gtf.len <- nrow(gtf) # customize second matrix block second <- data.frame(matrix("",gtf.len,2)) for(i in 1:gtf.len){ second[i,] <- unlist( strsplit(gtf[i,9], split=";", fixed=T ))[c(1,3)] # char df } # create a mapping table go <- seq(1,gtf.len,3) go.len <- length(go) # ENST ENSG symbol map.df <- data.frame(ensembl_transcript_id=rep("x",go.len), ensembl_gene_id=rep("x",go.len), hgnc_symbol=rep("x",go.len)) # fill symerr <- 0 k <- 1 for(i in go){ # mapping # first row map.df[k, 2] <- sub('gene_id ', '', second[i, 1], fixed=T) # gene_source mirbase gene_source mirbase tmp <- sub(' gene_name ', '', second[i, 2], fixed=T) if(tmp==" gene_source mirbase"){ map.df[k, 3] <- "no_name" symerr <- symerr+1 }else{ map.df[k, 3] <- tmp } # second row map.df[k, 1] <- sub(' transcript_id ', '', second[(i+1), 2], fixed=T) k <- k+1 } tmp1 <- map.df[order(map.df[, 1]), 1][c(1,go.len)] tmp2 <- map.df[order(map.df[, 2]), 2][c(1,go.len)] write.table(map.df, mfout, append=F, quote=F, sep="\t", row.names=F, col.names=T) cat("\n gtf rows in ",gtf.len," map rows out ",go.len) cat("\n ENSG:",tmp2," ENST:",tmp1," symbol_err:",symerr,"'no_name'\n") return() } #a <- create.gtf.ens.map(fin="/home/korschi/1/Homo_sapiens.GRCh38.111.mirbase.gtf", mfout="/home/korschi/1/m2_ens.gtf.map.txt") # gtf rows in 5637 map rows out 1879 # ENSG: ENSG00000194717 ENSG00000292355 ENST: ENST00000290239 ENST00000711167 symbol_err: 25 'no_name' map.enst.symbol <- function(x, y){ # Ensembl specific ENST to hgnc symbol translator # x: vector with ENST ids # y: mapping table: df: ENST ENSG symbol [ ensembl_transcript_id ensembl_gene_id hgnc_symbol ] nrx <- nrow(x) nry <- nrow(y) # customize second matrix block second <- data.frame(matrix("",gtf.len,2)) for(i in 1:gtf.len){ second[i,] <- unlist( strsplit(gtf[i,9], split=";", fixed=T ))[c(1,3)] # char df } # create a mapping table go <- seq(1,gtf.len,3) go.len <- length(go) # ENST ENSG symbol map.df <- data.frame(ensembl_transcript_id=rep("x",go.len), ensembl_gene_id=rep("x",go.len), hgnc_symbol=rep("x",go.len)) # fill symerr <- 0 k <- 1 for(i in go){ # mapping # first row map.df[k, 2] <- sub('gene_id ', '', second[i, 1], fixed=T) # gene_source mirbase gene_source mirbase tmp <- sub(' gene_name ', '', second[i, 2], fixed=T) if(tmp==" gene_source mirbase"){ map.df[k, 3] <- "no_name" symerr <- symerr+1 }else{ map.df[k, 3] <- tmp } # second row map.df[k, 1] <- sub(' transcript_id ', '', second[(i+1), 2], fixed=T) k <- k+1 } tmp1 <- map.df[order(map.df[, 1]), 1][c(1,go.len)] tmp2 <- map.df[order(map.df[, 2]), 2][c(1,go.len)] write.table(map.df, mfout, append=F, quote=F, sep="\t", row.names=F, col.names=T) cat("\n gtf rows in ",gtf.len," map rows out ",go.len) cat("\n ENSG:",tmp2," ENST:",tmp1," symbol_err:",symerr,"'no_name'\n") return() } #a <- map.enst.symbol(x, y=map.m2.ens.gtf.map) # m2_ens.gtf.map.txt, hsa_mir_gtf_map.txt, m2_tRNA.gtf.map.txt #map.m2.ens.gtf.map <- read.table("m2_ens.gtf.map.txt", header=F, sep="\t") #map.hsa.mir.gtf.map <- read.table("hsa_mir_gtf_map.txt", header=F, sep="\t") #map.m2.tRNA.gtf.map <- read.table("m2_tRNA.gtf.map.txt", header=F, sep="\t") ## tRNA ## # from gencode 45 based on Ensemble 111 # chr1 ENSEMBL tRNA 16520585 16520658 . - . gene_id "20658"; transcript_id "20658"; gene_type "Asn_tRNA"; gene_name "20658"; transcript_type "Asn_tRNA"; transcript_name "20658"; level 3; # chr1 ENSEMBL tRNA 16532398 16532471 . - . gene_id "20657"; transcript_id "20657"; gene_type "Asn_tRNA"; gene_name "20657"; transcript_type "Asn_tRNA"; transcript_name "20657"; level 3; # # Ensemble # gene # 1 havana gene 182696 184174 . + . gene_id "ENSG00000279928"; gene_version "2"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; # 1 havana transcript 182696 184174 . + . gene_id "ENSG00000279928"; gene_version "2"; transcript_id "ENST00000624431"; transcript_version "2"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "DDX11L17-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; # 1 havana exon 182696 182746 . + . gene_id "ENSG00000279928"; gene_version "2"; transcript_id "ENST00000624431"; transcript_version "2"; exon_number "1"; gene_name "DDX11L17"; gene_source "havana"; gene_biotype "unprocessed_pseudogene"; transcript_name "DDX11L17-201"; transcript_source "havana"; transcript_biotype "unprocessed_pseudogene"; exon_id "ENSE00003759020"; exon_version "2"; tag "basic"; tag "Ensembl_canonical"; transcript_support_level "NA"; # remove header first convert.gtf.gtf.trna <- function(fin, fout, mfout, ENSGstart, ENSTstart){ # tRNA specific converter gtf to gtf-ensemble (for gencode gtf file) # fin, fout, mfout: gtf path/file, gtf path/file, map path/file # ENSGstart, ENSTstart: number with 11 positions # see analytics in gtf_info.txt # e.g. at the moment: ENSG 302000 ENST 724000 arbitrarily chosen - should work # because only the id links, name and accession/chr is used - increase/decrease of digit range will raise code problems (!) gtf <- read.table(fin, header=F, sep="\t") gtf.len <- nrow(gtf) # customize second matrix block second <- data.frame(matrix("",gtf.len,7)) for(i in 1:gtf.len){ second[i,] <- unlist( strsplit(gtf[i,9], split=";", fixed=T ) ) # char df } # shrink basic matrix now gtf <- gtf[ , -9] # create empty data structure with dummy data new.df.len <- 3*gtf.len new.df <- data.frame(rep("x",new.df.len), rep("x",new.df.len), rep("x",new.df.len), rep("x",new.df.len), rep("x",new.df.len), rep(".",new.df.len), rep("+-",new.df.len), rep(".",new.df.len), rep("x",new.df.len) ) # create additionally a final mapping table for the salmon analysis # ENST ENSG symbol map.df <- data.frame(ensembl_transcript_id=rep("x",gtf.len), ensembl_gene_id=rep("x",gtf.len), hgnc_symbol=rep("x",gtf.len)) # fill k <- 1 for(i in 1:gtf.len){ tmp1 <- sub('chr', '', gtf[i, 1], fixed=T) # first row new.df[k, 1] <- tmp1 new.df[k, 2] <- "GENCODE" new.df[k, 3] <- "gene" new.df[k, 4] <- gtf[i, 4] new.df[k, 5] <- gtf[i, 5] new.df[k, 6] <- "." new.df[k, 7] <- gtf[i, 7] new.df[k, 8] <- "." tmp2 <- sub(' gene_name ', '', second[i, 4], fixed=T) new.df[k, 9] <- paste('gene_id "ENSG00000',ENSGstart,'"; gene_name "',tmp2,'"; gene_source "GENCODE"; gene_biotype "tRNA";', sep="") # tRNA not ENsemble official # mapping tmp3 <- sub('gene_id ', '', second[i, 1], fixed=T) map.df[i, 1] <- paste("ENST00000",ENSTstart, sep="") map.df[i, 2] <- paste("ENSG00000",ENSGstart, sep="") map.df[i, 3] <- paste("tRNA_",tmp3,"_",tmp2, sep="") k <- k+1 # second row new.df[k, 1] <- tmp1 new.df[k, 2] <- "GENCODE" new.df[k, 3] <- "transcript" new.df[k, 4] <- gtf[i, 4] new.df[k, 5] <- gtf[i, 5] new.df[k, 6] <- "." new.df[k, 7] <- gtf[i, 7] new.df[k, 8] <- "." # not the full Ensemble info new.df[k, 9] <- paste('gene_id "ENSG00000',ENSGstart,'"; transcript_id "ENST00000',ENSTstart,'"; gene_name "',tmp2,'"; gene_source "GENCODE"; gene_biotype "tRNA";', sep="") k <- k+1 # third row new.df[k, 1] <- tmp1 new.df[k, 2] <- "GENCODE" new.df[k, 3] <- "exon" new.df[k, 4] <- gtf[i, 4] new.df[k, 5] <- gtf[i, 5] new.df[k, 6] <- "." new.df[k, 7] <- gtf[i, 7] new.df[k, 8] <- "." # not the full Ensemble info new.df[k, 9] <- paste('gene_id "ENSG00000',ENSGstart,'"; transcript_id "ENST00000',ENSTstart,'"; exon_number "1"; gene_name "',tmp2,'"; gene_source "GENCODE"; gene_biotype "tRNA";', sep="") k <- k+1 ENSGstart <- ENSGstart+1 ENSTstart <- ENSTstart+1 } write.table(new.df, fout, append=F, quote=F, sep="\t", row.names=F, col.names=F) write.table(map.df, mfout, append=F, quote=F, sep="\t", row.names=F, col.names=T) cat("\n gtf rows in ",gtf.len," gtf rows out ",new.df.len) cat("\n ENSG final ",ENSGstart," ENST final ",ENSTstart," control and for next chunck of data...\n") return() } #a <- convert.gtf.gtf.trna(fin="/home/korschi/1/gencode.v45.tRNAs.gtf", # fout="/home/korschi/1/gencode.v45.tRNAs.ens.gtf", # mfout="/home/korschi/1/gencode.v45.tRNAs.map.txt", # ENSGstart=302000, ENSTstart=724000) # gtf rows in 649 gtf rows out 1947 # ENSG final 302649 ENST final 724649 control and for next chunck of data... # run fasta_keep3() (with m_tRNA.fasta) # correct manually m_tRNA_info.txt (fasta description line) # seq lines removed, additionally two lines adjusted: ..KI270713v1.. -> KI270713.1, chr removed before number # run convert.faH.gtf.trna() (on adjusted m_tRNA_info.txt) # >Homo_sapiens_tRNA-Ala-AGC-1-1 (tRNAscan-SE ID: chr6.trna116) Ala (AGC) 72 bp mature sequence Sc: 84.9 chr6:28795964-28796035 (-) # >Homo_sapiens_tRNA-Ala-AGC-10-1 (tRNAscan-SE ID: chr6.trna20) Ala (AGC) 73 bp mature sequence Sc: 60.1 chr6:26687257-26687329 (+) # >Homo_sapiens_tRNA-Ala-AGC-10-2 (tRNAscan-SE ID: chr6.trna169) Ala (AGC) 73 bp mature sequence Sc: 60.1 chr6:26814339-26814411 (-) # >Homo_sapiens_tRNA-Ala-AGC-11-1 (tRNAscan-SE ID: chr6.trna179) Ala (AGC) 73 bp mature sequence Sc: 58.9 chr6:26571864-26571936 (-) convert.faH.gtf.trna <- function(fin, fout, mfout, ENSGstart, ENSTstart){ # tRNA specific converter, fasta description line to gtf-ensemble (based on m_miRtRNA.fasta) # fin, fout, mfout: gtf path/file, gtf path/file, map path/file # ENSGstart, ENSTstart: number with 11 positions # see analytics in gtf_info.txt # e.g. at the moment: ENSG 302000 ENST 724000 arbitrarily chosen - should work # because only the id links, name and accession/chr is used - increase/decrease of digit range will raise code problems (!) gtf <- read.table(fin, header=F, sep=" ") gtf.len <- nrow(gtf) # customize V13, V14 second <- data.frame(matrix("",gtf.len,4)) for(i in 1:gtf.len){ tmpCC <- unlist( strsplit(gtf[i,13], split=":", fixed=T ) ) second[i,1:3] <- c( tmpCC[1], unlist(strsplit(tmpCC[2], split="-", fixed=T)) ) second[i,4] <- substr(gtf[i,14], 2, 2) } # create empty data structure with dummy data new.df.len <- 3*gtf.len new.df <- data.frame(rep("x",new.df.len), rep("x",new.df.len), rep("x",new.df.len), rep("x",new.df.len), rep("x",new.df.len), rep(".",new.df.len), rep("+-",new.df.len), rep(".",new.df.len), rep("x",new.df.len) ) # create additionally a final mapping table for the salmon analysis # ENST ENSG symbol map.df <- data.frame(ensembl_transcript_id=rep("x",gtf.len), ensembl_gene_id=rep("x",gtf.len), hgnc_symbol=rep("x",gtf.len)) # fill k <- 1 for(i in 1:gtf.len){ tmp1 <- second[i, 1] # first row new.df[k, 1] <- tmp1 new.df[k, 2] <- "Ensembl" new.df[k, 3] <- "gene" new.df[k, 4] <- second[i, 2] new.df[k, 5] <- second[i, 3] new.df[k, 6] <- "." new.df[k, 7] <- second[i, 4] new.df[k, 8] <- "." tmp2 <- sub('>', '', gtf[i, 1], fixed=T) new.df[k, 9] <- paste('gene_id "ENSG00000',ENSGstart,'"; gene_name "',tmp2,'"; gene_source "Ensembl"; gene_biotype "tRNA";', sep="") # tRNA not ENsemble official # mapping tmp3 <- sub('gene_id ', '', second[i, 1], fixed=T) map.df[i, 1] <- paste("ENST00000",ENSTstart, sep="") map.df[i, 2] <- paste("ENSG00000",ENSGstart, sep="") map.df[i, 3] <- tmp2 k <- k+1 # second row new.df[k, 1] <- tmp1 new.df[k, 2] <- "Ensembl" new.df[k, 3] <- "transcript" new.df[k, 4] <- second[i, 2] new.df[k, 5] <- second[i, 3] new.df[k, 6] <- "." new.df[k, 7] <- second[i, 4] new.df[k, 8] <- "." # not the full Ensemble info new.df[k, 9] <- paste('gene_id "ENSG00000',ENSGstart,'"; transcript_id "ENST00000',ENSTstart,'"; gene_name "',tmp2,'"; gene_source "Ensembl"; gene_biotype "tRNA";', sep="") k <- k+1 # third row new.df[k, 1] <- tmp1 new.df[k, 2] <- "Ensembl" new.df[k, 3] <- "exon" new.df[k, 4] <- second[i, 2] new.df[k, 5] <- second[i, 3] new.df[k, 6] <- "." new.df[k, 7] <- second[i, 4] new.df[k, 8] <- "." # not the full Ensemble info new.df[k, 9] <- paste('gene_id "ENSG00000',ENSGstart,'"; transcript_id "ENST00000',ENSTstart,'"; exon_number "1"; gene_name "',tmp2,'"; gene_source "Ensembl"; gene_biotype "tRNA";', sep="") k <- k+1 ENSGstart <- ENSGstart+1 ENSTstart <- ENSTstart+1 } write.table(new.df, fout, append=F, quote=F, sep="\t", row.names=F, col.names=F) write.table(map.df, mfout, append=F, quote=F, sep="\t", row.names=F, col.names=T) cat("\n gtf rows in ",gtf.len," gtf rows out ",new.df.len) cat("\n ENSG final ",ENSGstart," ENST final ",ENSTstart," control and for next chunck of data...\n") return() } #convert.faH.gtf.trna(fin="/home/korschi/1/m_tRNA_info.txt", # fout="/home/korschi/1/m_tRNA.gtf", # mfout="/home/korschi/1/m_tRNA.gtf.map.txt", # ENSGstart=302000, ENSTstart=724000) # gtf rows in 432 gtf rows out 1296 # ENSG final 302432 ENST final 724432 control and for next chunck of data... #convert.faH.gtf.trna(fin="/home/korschi/1/m2_tRNA_info.txt", # fout="/home/korschi/1/m2_tRNA.gtf", # mfout="/home/korschi/1/m2_tRNA.gtf.map.txt", # ENSGstart=302000, ENSTstart=724000) # gtf rows in 260 gtf rows out 780 # ENSG final 302260 ENST final 724260 control and for next chunck of data... ## miR tRNA correct pp pysical positions # hsa-miR-12136 8 100.000 16 0 0 1 16 95751483 95751498 1.2 32.2 # hsa-miR-12136 7 100.000 16 0 0 2 17 132208378 132208393 1.2 32.2 # Homo_sapiens_tRNA-Ala-AGC-1-1 6 100.000 72 0 0 1 72 28796035 28795964 1.26e-32 143 # Homo_sapiens_tRNA-Ala-AGC-1-1 5 89.091 55 6 0 8 62 181206875 181206929 3.70e-08 61.9 # blast - header: default tabular format "6" shown with "7": # Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score # note: correction of pp is not working for tRNA with multiple identical sequences # in different locations in close neighborhood !! # - blast : is not finding in genome fasta 111 # these multiple identical sequences in different locations in close neighborhood - only one (10 hit list) correct.pp.gtf.ens <- function(fio, fin, v=F){ # gtf.ens specific - correct physical position by blast run information on newest Ensemble111 fasta # fio, fin : gtf.ens path/file, blast path/file # blast with 10 return hits, otherwise no hit is possible # v: verbose: T,F gtf.ens <- read.table(fio, header=F, sep="\t", quote="'" ) # start end always increasing, quote ok blast <- read.table(fin, header=F, sep="\t", quote="" ) # start/end reversed according to (+-), quote ok gtf.ens.l <- nrow(gtf.ens) blast.l <- nrow(blast) # flip some entries so that start end is the right order h <- 0 dd <- (blast[,10]-blast[,9]) <0 for(i in 1:blast.l){ if(dd[i]==T){ a <- blast[i,9] blast[i,9] <- blast[i,10] blast[i,10] <- a h <- h+1 } } # start with gtf.ens coll <- NULL # collect differences coll1 <- NULL # collect problems s1 <- 0 # no changes s2a <- 0 # changes s2b <- 0 # changes s3 <- 0 # no blast hit gel <- seq(1, gtf.ens.l, 3) # check 3 for(i in gel){ # get chr,name,present positions g.chr <- gtf.ens[i,1] # char g.pos <- unlist(gtf.ens[i,4:5]) # int g.n <- sub(pattern=" gene_name ","", x=unlist( strsplit(gtf.ens[i,9], split=";", fixed=T ) )[2], fixed=T) # gtf.ens[,9] content has "" g.n <- gsub(pattern='[\"]','', x=g.n, perl=T) mat <- blast[blast[,1] %in% g.n, ] # filter by name-id mat2 <- mat[mat[,2] %in% g.chr,,drop=F] # filter by chr if(v){ cat("\n i",i," g.n",g.n," g.pos",g.pos," g.chr",g.chr,"\n"); print(mat); print(mat2) } nr <- nrow(mat2) if(nr==1){ # one hit b.pos <- unlist(mat2[,9:10]) # int # test if(sum(g.pos==b.pos)==2){ coll <- rbind(coll, c(g.n,g.chr,g.pos,b.pos,0)) s1 <- s1+1 }else{ coll <- rbind(coll, c(g.n,g.chr,g.pos,b.pos,1)) gtf.ens[i, 4:5] <- b.pos gtf.ens[(i+1), 4:5] <- b.pos gtf.ens[(i+2), 4:5] <- b.pos s2a <- s2a+1 } }else if(nr>1){ # more entries on one chr min.i <- which(mat2[,11] == min(mat2[,11])) # filter by p ranking very rare min.i.l <- length(min.i) if(min.i.l>1){ stop("\n equal ranks in blast table") }else{ b.pos <- mat2[min.i, 9:10] # int coll <- rbind(coll, c(g.n,g.chr,g.pos,b.pos,2)) gtf.ens[i, 4:5] <- b.pos gtf.ens[(i+1), 4:5] <- b.pos gtf.ens[(i+2), 4:5] <- b.pos s2b <- s2b+1 } }else{ # no blast hit - e.g. increase blast hits coll1 <- rbind(coll1, c(i,g.n,g.pos,g.chr)) cat("\n i",i," no entry in blast table") s3 <- s3+1 } #if(i==10){ break; return() } # gel dependent } dimnames(coll)[[2]] <- c("name-id","chr","start_g","end_g","start_b","end_b","flow") dimnames(coll1)[[2]] <- c("i","g.n","start_g","end_g","g.chr") write.table(gtf.ens, paste(fio,".gtf",sep=""), append=F, quote=F, sep="\t", row.names=F, col.names=F) write.table(coll, paste(fio,".log",sep=""), append=F, quote=F, sep="\t", row.names=F, col.names=T) write.table(coll1, paste(fio,".err",sep=""), append=F, quote=F, sep="\t", row.names=F, col.names=T) cat("\n blast positions: number of swaps",h,"/",blast.l) cat("\n gtf positions: no changes",s1,", changes-one hit",s2a,"-more hits",s2b,", no blast hit (no change)",s3,"/",length(gel),"\n") return() } #correct.pp.gtf.ens(fio="/home/korschi/1/mit.test.gtf", fin="/home/korschi/1/m_blast_mir_trna.txt", v=T) #correct.pp.gtf.ens(fio="/home/korschi/1/hsa_mir.gtf", fin="/home/korschi/1/m_blast_mir_trna.txt", v=F) #correct.pp.gtf.ens(fio="/home/korschi/1/mt.Homo_sapiens.miR.gtf", fin="/home/korschi/1/m2_blast_mir_trna.txt", v=F) # ... i 5044 no entry in blast table # blast positions: number of swaps 14559 / 29200 # gtf positions: no changes 2532 , changes-one hit 99 -more hits 0 , no blast hit (no change) 22 / 2653 #correct.pp.gtf.ens(fio="/home/korschi/1/mt.Homo_sapiens.miRtRNA.gtf", fin="/home/korschi/1/m2_blast_mir_trna.txt", v=F) # ... i 5044 no entry in blast table # blast positions: number of swaps 14559 / 29200 # gtf positions: no changes 2770 , changes-one hit 121 -more hits 0 , no blast hit (no change) 22 / 2913 check.pos.order.gtf.ens <- function(fin="/home/korschi/1/mit.Homo_sapiens.miRtRNA.gtf"){ # check the order start end in the ensemble like gtf file # start end should always be increasing gtf.ens <- read.table(fin, header=F, sep="\t", quote="'" ) gtf.ens.l <- nrow(gtf.ens) work <- NULL s <- 0 gel <- seq(1, gtf.ens.l, 3) # check 3 sub <- c(0,1,2) for(i in gel){ for(j in sub){ start <- gtf.ens[(i+j),4] end <- gtf.ens[(i+j),5] if(end<=start){ cat("\n block start",i," err row",(i+j)," start",start," end",end) work <- rbind(work, gtf.ens[(i+j), ]) s <- s+1 } } } cat("\n gtf length: ",gtf.ens.l,", div by 3:",length(gel),", hits: end<=start ",s,"\n") return(work) } #a <- check.pos.order.gtf.ens(fin="/home/korschi/1/mit.Homo_sapiens.miRtRNA.gtf") # ... # block start 9949 err row 9950 start 26286597 end 26286526 # gtf length: 9957 , div by 3: 3319 , hits: end<=start 356 #a <- read.table("/home/korschi/1/m_blast_mir_trna.txt", header=F, sep="\t", quote="'" ) import.ens.ids <- function(fin, cols=13){ # miR specific importer Ensemble gtf for original ENST ENSG ids # fin: gtf path/file, cols: number of columns in the last special column (;) # ENSG, ENST: number with 11 positions gtf <- read.table(fin, header=F, sep="\t") gtf.len <- nrow(gtf) # customize second matrix block second <- data.frame(matrix("", gtf.len, cols)) for(i in 1:gtf.len){ second[i,] <- unlist( strsplit(gtf[i,9], split=";", fixed=T ) ) # char df } # customize second <- second[, c(1,3,5,7)] # for(i in 1:gtf.len){ second[i, 1] <- gsub(pattern='gene_id ','', x=second[i, 1], perl=T) second[i, 2] <- gsub(pattern=' transcript_id ','', x=second[i, 2], perl=T) second[i, 3] <- gsub(pattern=' gene_name ','', x=second[i, 3], perl=T) second[i, 4] <- gsub(pattern=' gene_biotype ','', x=second[i, 4], perl=T) } return(second) } #a <- import.ens.ids(fin="/home/korschi/1/hu.mirna.list3.gtf", cols=13) gff.create.bed <- function(fin, fout){ # mirbase gff to bed format - region file for IGV # fin, mfout: gff path/file, bed path/file gff <- read.table(fin, header=F, sep="\t") # remove lines 'primary_transcript' logi1 <- grepl("primary_transcript", x=gff[,3]) gff <- gff[!logi1, ] gff.len <- nrow(gff) # customize second matrix block second <- data.frame(matrix("",gff.len,1)) for(i in 1:gff.len){ second[i,] <- unlist( strsplit(gff[i,9], split=";", fixed=T ))[3] # char df } # create a bed table bed.df <- data.frame(chr_name=rep("x",gff.len), start=rep("x",gff.len), end=rep("x",gff.len), name=rep("x",gff.len)) # fill for(i in 1:gff.len){ bed.df[i, 1] <- sub('chr', '', gff[i,1], fixed=T) bed.df[i, 2] <- gff[i,4] bed.df[i, 3] <- gff[i,5] bed.df[i, 4] <- sub('Name=', '', second[i, 1], fixed=T) } write.table(bed.df, fout, append=F, quote=F, sep="\t", row.names=F, col.names=T) cat("\n gff rows processed ",gff.len,"\n") return() } #gff.create.bed(fin="/home/korschi/1/hsa_mir.gff3", fout="/home/korschi/1/hsa_mir_gff3.bed") # gff rows processed 2880 gtf.create.bed <- function(fin, fout){ # tRNA gtf to bed format - region file for IGV # fin, mfout: gff path/file, bed path/file gtf <- read.table(fin, header=F, sep="\t") # remove lines 'primary_transcript' logi1 <- grepl("exon", x=gtf[,3]) gtf <- gtf[logi1, ] gtf.len <- nrow(gtf) # customize second matrix block second <- data.frame(matrix("",gtf.len,1)) for(i in 1:gtf.len){ second[i,] <- unlist( strsplit(gtf[i,9], split=";", fixed=T ))[4] # char df } # create a bed table bed.df <- data.frame(chr_name=rep("x",gtf.len), start=rep("x",gtf.len), end=rep("x",gtf.len), name=rep("x",gtf.len)) # fill for(i in 1:gtf.len){ bed.df[i, 1] <- sub('chr', '', gtf[i,1], fixed=T) bed.df[i, 2] <- gtf[i,4] bed.df[i, 3] <- gtf[i,5] bed.df[i, 4] <- sub('gene_name ', '', second[i, 1], fixed=T) } write.table(bed.df, fout, append=F, quote=F, sep="\t", row.names=F, col.names=T) cat("\n gtf rows processed ",gtf.len,"\n") return() } #gtf.create.bed(fin="/home/korschi/1/m2_tRNA.gtf", fout="/home/korschi/1/m2_tRNA_gtf.bed") # gtf rows processed 260 find.in.bed <- function(x, inF="", x.name.c, b.name.c){ # get ids and search in a BED file # return a matrix # x: vector of ids or df with name column, inF: path/file of BED # x.name.c,b.name.c: column number - x.name.c only useful if df cv <- is.vector(x) cdf <- is.data.frame(x) bed <- read.table(inF, header=F, sep="\t") if(cv & !cdf){ xl <- length(x) le <- NULL for(i in 1:xl){ le <- rbind(le, bed[bed[,b.name.c]%in%x[i],] ) } }else if(!cv & cdf){ ncx <- ncol(x) cat("\n x rows in",nrow(x)) le <- merge(x=x, y=bed, by.x=x.name.c, by.y=b.name.c) cat("\n x rows out",nrow(le)) le <- le[order(le[,2],decreasing=T),] # 2: first count col dimnames(le)[[2]][(ncx+1):(ncx+5)] <- c("chr","start","end","x","strand") } return(le) } #a4 <- find.in.bed(x=a3[1:10,1], inF="/home/korschi/1/sources/RNAcentral/r_star_pirbase_gold.bed", b.name.c=4) #a4 <- find.in.bed(x=a3[1:10,], inF="/home/korschi/1/sources/RNAcentral/r_star_pirbase_gold.bed", x.name.c=1, b.name.c=4)