M.tuberculosis: Annotation data for proteomics analyses (uniprot).

1 Annotation version: 20171205

1.1 Genome annotation input

1.1.1 Read a gff file

In contrast, it is possible to load most annotations of interest directly from the gff files used in the alignments. More in-depth information for the human transcriptome may be extracted from biomart.

## The old way of getting genome/annotation data
mtb_gff <- "reference/mycobacterium_tuberculosis_h37rv_2.gff.gz"

mtb_genome <- "reference/mtuberculosis_h37rv_genbank.fasta"
mtb_cds <- "reference/mtb_cds.fasta"

mtb_annotations <- sm(load_gff_annotations(mtb_gff, type="gene"))
rownames(mtb_annotations) <- mtb_annotations[["ID"]]

2 A quick pattern match

This little block is intended to seek out peptide sequences with the highly degenerate pattern: Y(E|D) with a varying number of amino acids between the Y and (E or D).

mtb_cds <- "reference/mtb_cds.fasta"

get_hits <- function(patterns, peptide_file) {
  pep_seq <- Biostrings::readAAStringSet(peptide_file)
  pep_lst <- as.data.frame(pep_seq)[["x"]]
  position_df <- data.frame(row.names=names(pep_seq))
  number_df <- data.frame(row.names=names(pep_seq))
  pct_df <- data.frame(row.names=names(pep_seq))
  nt_df <- data.frame(row.names=names(pep_seq))
  ct_df <- data.frame(row.names=names(pep_seq))
  for (pat in patterns) {
    column <- gregexpr(pattern=pat, text=pep_lst)
    names(column) <- names(pep_seq)
    for (r in 1:length(column)) {
      name <- names(column)[r]
      row <- column[[r]]
      len <- nchar(pep_lst[r])
      pct <- as.numeric(row) / len
      ## An important caveat here:  There may be more than one value.
      op <- options(warn=2)
      if (pct[1] < 0) {
        pct_df[r, pat] <- -1
      } else {
        pct_df[r, pat] <- toString(pct)
      }
      options(op)
      nt_str <- ""
      ct_str <- ""
      for (i in pct) {
        if (i < 0) {
          next
        }
        if (i < 0.1) {
          nt_str <- paste0(nt_str, ", ", i)
        }
        if (i > 0.9) {
          ct_str <- paste0(ct_str, ", ", i)
        }
      }
      nt_str <- gsub(pattern="^\\, ", replacement="", x=nt_str)
      ct_str <- gsub(pattern="^\\, ", replacement="", x=ct_str)
      nt_df[r, pat] <- nt_str
      ct_df[r, pat] <- ct_str
      position_df[r, pat] <- toString(as.numeric(row))
      number_df[r, pat] <- length(as.numeric(row))
    }
  }
  retlst <- list(
    "numbers" = number_df,
    "pct" = pct_df,
    "nterminals" = nt_df,
    "cterminals" = ct_df,
    "positions" = position_df)
  return(retlst)
}
patterns <- c("Y(E|D)", "Y.(E|D)", "Y..(E|D)", "Y...(E|D)", "Y....(E|D)", "Y.....(E|D)")
stuff <- get_hits(patterns, mtb_cds)
write.csv(file="excel/positions_by_pattern.csv", stuff[["positions"]])
write.csv(file="excel/number_by_pattern.csv", stuff[["numbers"]])
write.csv(file="excel/pct_by_pattern.csv", stuff[["pct"]])
write.csv(file="excel/nterminals.csv", stuff[["nterminals"]])
write.csv(file="excel/cterminals.csv", stuff[["cterminals"]])

2.0.1 Download from microbesonline

## First figure out the ID for the Mtb genome:
ids <- get_microbesonline_ids("37")
## Mycobacterium tuberculosis H37Rv is 83332
mtb_microbes <- load_microbesonline_annotations(ids=83332)
mtb_go <- load_microbesonline_go(id=83332)

mtb_uniprotws <- load_uniprotws_annotations(species="Mycobacterium tuberculosis")

## More than 1 species was returned, they follow; the first was chosen arbitrarily.

##   taxon ID
## 1   419947
## 2   443149
## 3   443150
## 4   652616
## 5   336982
## 6   478434
## 7     1773
##                                                        Species name
## 1            Mycobacterium tuberculosis (strain ATCC 25177 / H37Ra)
## 2                      Mycobacterium tuberculosis (strain CCDC5079)
## 3                      Mycobacterium tuberculosis (strain CCDC5180)
## 4 Mycobacterium tuberculosis (strain ATCC 35801 / TMC 107 / Erdman)
## 5                           Mycobacterium tuberculosis (strain F11)
## 6                Mycobacterium tuberculosis (strain KZN 1435 / MDR)
## 7                                        Mycobacterium tuberculosis

## Getting mapping data for A5U205 ... and P_GI

if (!isTRUE(get0("skip_load"))) {
  pander::pander(sessionInfo())
  message(paste0("This is hpgltools commit: ", get_git_commit()))
  this_save <- paste0(gsub(pattern="\\.Rmd", replace="", x=rmd_file), "-v", ver, ".rda.xz")
  message(paste0("Saving to ", this_save))
  tmp <- sm(saveme(filename=this_save))
}

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset 7de4503f6bb5724c28cce24af5dbee22bb1c0cae

## R> packrat::restore()

## This is hpgltools commit: Thu Apr 12 22:08:53 2018 -0400: 7de4503f6bb5724c28cce24af5dbee22bb1c0cae

## Saving to 01_annotation-v20171205.rda.xz