# Question
Title: Retrieval of genomic position by using biomaRt package
I have a list of several protein names, their primary gene name and their amino acid sequence (as extracted from the 'SEQRES' section in their corresponding PDB files).

I'm looking for the genomic position (triad of translated nucleotides) of a specific amino acid in the sequence (i.e. residue number 64 in a protein sequence of 344 amino acids).

Is it possible to BLAT an amino acid sequence by using biomaRt in R? Do you suggest other methods to achieve this task?

# Answer
> 4 votes
BiomaRt cannot do that I'm afraid. There is an Ensembl REST API endpoint that will get you the genomic location of protein coordinates. It needs an Ensembl peptide ID an input though, so you could use the xref endpoint first to get that. There's a bit of sample code in R on those pages, you could use that to script together if you've got a list of these.

---
Tags: r, proteins, ngs, pdb, biomart
---


# Question
Title: GOseq analysis with evidence code filter
When using GOseq analysis on RNA-seq data, I often find many 'false positives'. What I mean with that is that some genes in certain categories are not really involved in the process, but are only involved due to 'Electronic Annotation' (IEA). Most of the genes are IEA. A list of all evidence codes is here.

Is it possible to (easily) do GOseq analysis without this evidence group?

I mean, I know how to do it the hard way (make a custom GO-genes annotation in e.g., biomaRt removing all IEA evidence code genes, and perform GOseq with custom genome). But is there also a easy way in GOseq available?

# Answer
A bit weird to answer my own question, but for who's interested I can show the 'hard' way using `biomaRt` and then use `goseq` with `gene2cat` defined.

So first you need biomaRt to get annotation of the genes. I have ensembl genes from mm10.

```
library(biomaRt)

table <- read.table(file, header = T, stringsAsFactors = F)

ensembl <- useMart("ENSEMBL_MART_ENSEMBL")

ensembl <- useDataset("mmusculus_gene_ensembl", mart = ensembl)

allGenes <- table$Geneid

GOs <- getBM(attributes=c('ensembl_gene_id', 'go_id', 'name_1006', 'go_linkage_type', 'namespace_1003'),
      filters = 'ensembl_gene_id',
      values = allGenes,
      mart = ensembl)

dim(GOs)
[1] 477082      5

table(GOs$go_linkage_type)
          EXP    IBA     IC    IDA    IEA    IEP    IGI    IMP    IPI    ISA    ISM    ISO    ISS    NAS     ND    RCA    TAS
 36931     36  27918    525  34038 210114    630   6891  30219   6205   4234     20  77365  24673    609  12821     99   3754

# I only want BP and no IEA evidence code
GOs_no_IEA <- GOs[GOs$go_linkage_type != "IEA" & GOs$namespace_1003 == "biological_process",]

dim(GOs_no_IEA)
[1] 111944      5

```

Then you can continue with goseq.

```
library(goseq)

All <- read.table(input_All, sep = "\t", header = T)

All_genes <- as.character(All$Geneid)

TopTable <- read.table(input, sep = "\t", header = T)

sig_genes <- as.character(TopTable$Ensembl)

# Make gene vector

gene.vector <- as.integer(All_genes %in% sig_genes)

names(gene.vector) <- All_genes

gene.lengths <- All$Length

pwf <- nullp(gene.vector, "mm10", "ensGene", bias.data = gene.lengths)
head(pwf)

GO.wall <- goseq(pwf, "mm10", "ensGene", gene2cat = GOs_no_IEA[, c(1:2)], use_genes_without_cat = F, test.cats = "GO:BP")

```
> 4 votes

---
Tags: rna-seq, go-enrichment, go
---


# Question
Title: How to split multifasta based on partial fasta header
I have download a batch of refseq fasta files and want to split them based on strain. this is complicated by accession numbers, project IDs etc. I downloaded the genomes via eDirect from NCBI. I know I can split the output but I want them together for some preliminary parsing (remove ambiguous and plasmid sequences etc.)

See example below.

```
>NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
>NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
>NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
>NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC

```

I want two files: Clostridium botulinum BKT028387.fas and Peptoanaerobacter stomatis OBRC8.fas

My actual dataset is much larger (~70,000 genomes) and not too sure about the headers in each case.

---

Im thinking the following could work but Im not sure how to execute it:

Using a grep command to extract the lines but split the output.

If I download the summary list from refseq and use the names as the query example:

```
grep -Fwf refseq.list sequences.fa

```

> Summary file
>
> 1. Peptoanaerobacter stomatis strain OBRC8, whole genome shotgun sequencing project 7,810,226 rc other DNA NZ\_AJXR00000000.1 GI:938474759
>
> 2. Clostridium botulinum BKT028387, whole genome shotgun sequencing project 4,134,593 rc other DNA NZ\_NOJY00000000.1 GI:1231817869
>
> But I could edit it down to
>
> Peptoanaerobacter stomatis strain OBRC8
>
> Clostridium botulinum BKT028387

Am I one the right track. How do I finish this?

# Answer
> 1 votes
The SeqBuddy --pull\_records command can grab sequences based on a regular expression against record IDs and/or metadata. Given your example:

```
$: seqbuddy.py sequences.fa --pull_records "Clostridium botulinum.*BKT028387" full > Clostridium_botulinum_BKT028387.fas
$: seqbuddy.py sequences.fa --pull_records "Peptoanaerobacter stomatis.*OBRC8" full > Peptoanaerobacter_stomatis_OBRC8.fas

```

# Answer
> 1 votes
I had a similar problem and found a solution using awk. Since this is the first time I use awk, the syntax might not be perfect, but you should get an idea of what is possible.

```
awk -F' ' '{if(/^>/){sub("strain ", ""); name=$2" "$3" "$4; print > name".fa"}else{print > name".fa"}}' yourfile.fa

```

*(more information below)*

This should transform a file like this:

yourfile.fa

```
>NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
ATGCATGC
ATGCATGC
>NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
ATGCATGC
>NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
ATGCATGC
ATGCATGC
>NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC
ATGCATGC

```

Into two files like this:

'Clostridium botulinum BKT028387.fa'

```
>NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
ATGCATGC
ATGCATGC
>NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
ATGCATGC

```

'Peptoanaerobacter stomatis OBRC8.fa'

```
>NZ_ALNK01000035.1 Peptoanaerobacter stomatis OBRC8 ctg120007290634 ATGC
ATGCATGC
ATGCATGC
>NZ_ALNK01000031.1 Peptoanaerobacter stomatis OBRC8 ctg120007290635 ATGC
ATGCATGC

```

**How it works**

* `-F' '` tells awk to split the line by space (accessible with `$1`, `$2`, ...)
* `if(/^>/)` looks for line starting with the `>` character, then `sub("strain ", "")` remove the word `strain` and `name=$2" "$3" "$4; print > name".fa"` print the line into a new file with a name based on what you asked for.
* else if the line does not start with `>`, `{print > name".fa"}` appends it to a file using the current value of `name`.

This solution will work regardless of the number of lines between the fasta headers. Note that I won't recommend using spaces in filenames, so I would replace `name=$2" "$3" "$4` with something like `name=$2"_"$3"_"$4`. I used this solution to split a large fasta file (\>5 GB) containing \>5k genomes with \>25k scaffolds and I found it fast and rather memory efficient.

# Answer
> 0 votes
I would convert the fasta file to a format with one line per record (I have already posted the FastaToTbl and TblToFasta scripts I use for this before), and then look for a specific pattern:

```
FastaToTbl file.fa | grep 'BKT028387' | TblToFasta > clostridium.fa
FastaToTbl file.fa | grep 'OBRC8' | TblToFasta > peptoanaerobacter.fa

```

I have downloaded a few sequences using the same command you described in your previous question, and it looks like the headers are not standardized. For example, while you have many that are in the format of `$species strain $strain`, like these:

```
Clostridioides difficile strain 7496-D/ST9
Clostridioides difficile strain 5616-DH/ST42,
Peptoanaerobacter stomatis strain OBRC8
Romboutsia maritimum strain CCRI-22766
Tepidibacter mesophilus strain JCM

```

Others don't have the word `strain` and just give the strain:

```
Peptostreptococcus anaerobius 653-L contig00088,
Peptoclostridium acidaminophilum DSM 3953

```

So, as a first step, you could deal with the ones that do have the word `strain` as the 4th word of the header line (the 1st is the seqID which I have not included in the examples above) and the actual strain as the 5th and get those out of the way:

```
FastaToTbl foo.fa |
    awk '($4=="strain"){gsub("/","_"); print $0 >> $2"_"$3"_"$5".tbl"}'

```

That will look for lines whose 4th field is `strain`, then it will replace any `/` in the line with `_` (there's a strain called `CF/ST37` and that cannot be in a filename) and print into a file whose name is the 2nd, 3rd and 5th fields of the line concatenated. On my test file, this produced the following files (as you can see, some are wrong since there is a comma after the strain, but that's because they're not standardized):

```
$ ls *tbl
Clostridioides_difficile_5479-NonSp_ST229,.tbl  Clostridioides_difficile_6589-DH_ST42.tbl         Clostridioides_difficile_7457-NonSp_novelST,.tbl
Clostridioides_difficile_5479-NonSp_ST229.tbl   Clostridioides_difficile_6592-G_ST8,.tbl          Clostridioides_difficile_7457-NonSp_novelST.tbl
Clostridioides_difficile_5480-E_ST6,.tbl        Clostridioides_difficile_6592-G_ST8.tbl           Clostridioides_difficile_7458-N_ST10,.tbl
Clostridioides_difficile_5480-E_ST6.tbl         Clostridioides_difficile_6593-Y_ST110,.tbl        Clostridioides_difficile_7458-N_ST10.tbl
Clostridioides_difficile_5482-L_ST55,.tbl       Clostridioides_difficile_6593-Y_ST110.tbl         Clostridioides_difficile_7459-G_ST8,.tbl
Clostridioides_difficile_5482-L_ST55.tbl        Clostridioides_difficile_6595-AL_ST58,.tbl        Clostridioides_difficile_7459-G_ST8.tbl
Clostridioides_difficile_5485-Y_ST2,.tbl        Clostridioides_difficile_6595-AL_ST58.tbl         Clostridioides_difficile_7462-DH_ST42,.tbl
Clostridioides_difficile_5485-Y_ST2.tbl         Clostridioides_difficile_6596-D_ST9,.tbl          Clostridioides_difficile_7462-DH_ST42.tbl
Clostridioides_difficile_5488-DH_ST42,.tbl      Clostridioides_difficile_6596-D_ST9.tbl           Clostridioides_difficile_7463-E_ST6,.tbl
Clostridioides_difficile_5488-DH_ST42.tbl       Clostridioides_difficile_6599-NonSp_ST103,.tbl    Clostridioides_difficile_7463-E_ST6.tbl
Clostridioides_difficile_5489-M_ST15,.tbl       Clostridioides_difficile_6599-NonSp_ST103.tbl     Clostridioides_difficile_7464-Y_ST110,.tbl
Clostridioides_difficile_5489-M_ST15.tbl        Clostridioides_difficile_6600-G_ST8,.tbl          Clostridioides_difficile_7464-Y_ST110.tbl
Clostridioides_difficile_5490-CF_ST37,.tbl      Clostridioides_difficile_6600-G_ST8.tbl           Clostridioides_difficile_7466-DH_ST42,.tbl
Clostridioides_difficile_5490-CF_ST37.tbl       Clostridioides_difficile_6606-BK_ST11,.tbl        Clostridioides_difficile_7466-DH_ST42.tbl
Clostridioides_difficile_5493-DH_ST42,.tbl      Clostridioides_difficile_6606-BK_ST11.tbl         Clostridioides_difficile_7467-DH_ST42,.tbl
Clostridioides_difficile_5493-DH_ST42.tbl       Clostridioides_difficile_6607-DH_ST42,.tbl        Clostridioides_difficile_7467-DH_ST42.tbl
Clostridioides_difficile_5495-AF_ST41,.tbl      Clostridioides_difficile_6607-DH_ST42.tbl         Clostridioides_difficile_7468-NonSp_ST97,.tbl
Clostridioides_difficile_5495-AF_ST41.tbl       Clostridioides_difficile_6610-N_ST10,.tbl         Clostridioides_difficile_7468-NonSp_ST97.tbl
Clostridioides_difficile_5497-Y_ST2,.tbl        Clostridioides_difficile_6610-N_ST10.tbl          Clostridioides_difficile_7471-BI_ST1,.tbl
Clostridioides_difficile_5497-Y_ST2.tbl         Clostridioides_difficile_6612-DH_ST28,.tbl        Clostridioides_difficile_7471-BI_ST1.tbl
Clostridioides_difficile_5500-DH_ST42,.tbl      Clostridioides_difficile_6612-DH_ST28.tbl         Clostridioides_difficile_7472-Y_ST2,.tbl
Clostridioides_difficile_5500-DH_ST42.tbl       Clostridioides_difficile_6613-NonSp_ST3,.tbl      Clostridioides_difficile_7472-Y_ST2.tbl
Clostridioides_difficile_5504-L_ST55,.tbl       Clostridioides_difficile_6613-NonSp_ST3.tbl       Clostridioides_difficile_7473-DH_ST42,.tbl
Clostridioides_difficile_5504-L_ST55.tbl        Clostridioides_difficile_6614-AL_ST34,.tbl        Clostridioides_difficile_7473-DH_ST42.tbl
Clostridioides_difficile_5505-Y_ST2,.tbl        Clostridioides_difficile_6614-AL_ST34.tbl         Clostridioides_difficile_7476-Y_ST2,.tbl
Clostridioides_difficile_5505-Y_ST2.tbl         Clostridioides_difficile_6615-NonSp_ST11,.tbl     Clostridioides_difficile_7476-Y_ST2.tbl
Clostridioides_difficile_5537-D_ST9,.tbl        Clostridioides_difficile_6615-NonSp_ST11.tbl      Clostridioides_difficile_7478-NonSp_novelST,.tbl
Clostridioides_difficile_5537-D_ST9.tbl         Clostridioides_difficile_6616-NonSp_novelST,.tbl  Clostridioides_difficile_7478-NonSp_novelST.tbl
Clostridioides_difficile_5538-BK_ST11,.tbl      Clostridioides_difficile_6616-NonSp_novelST.tbl   Clostridioides_difficile_7479-NonSp_ST3,.tbl
Clostridioides_difficile_5538-BK_ST11.tbl       Clostridioides_difficile_6617-A_ST43,.tbl         Clostridioides_difficile_7479-NonSp_ST3.tbl
Clostridioides_difficile_5539-Y_ST2,.tbl        Clostridioides_difficile_6617-A_ST43.tbl          Clostridioides_difficile_7480-M_ST15,.tbl
Clostridioides_difficile_5539-Y_ST2.tbl         Clostridioides_difficile_6618-NonSp_ST53,.tbl     Clostridioides_difficile_7480-M_ST15.tbl
Clostridioides_difficile_5541-CF_ST37,.tbl      Clostridioides_difficile_6618-NonSp_ST53.tbl      Clostridioides_difficile_7481-NonSp_ST59,.tbl
Clostridioides_difficile_5541-CF_ST37.tbl       Clostridioides_difficile_6619-CF_ST37,.tbl        Clostridioides_difficile_7481-NonSp_ST59.tbl
Clostridioides_difficile_5542-Y_ST2,.tbl        Clostridioides_difficile_6619-CF_ST37.tbl         Clostridioides_difficile_7482-DH_ST42,.tbl
Clostridioides_difficile_5542-Y_ST2.tbl         Clostridioides_difficile_6621-NonSp_ST10,.tbl     Clostridioides_difficile_7482-DH_ST42.tbl
Clostridioides_difficile_5547-DH_ST42,.tbl      Clostridioides_difficile_6621-NonSp_ST10.tbl      Clostridioides_difficile_7486-A_ST43,.tbl
Clostridioides_difficile_5547-DH_ST42.tbl       Clostridioides_difficile_6622-CF_ST37,.tbl        Clostridioides_difficile_7486-A_ST43.tbl
Clostridioides_difficile_5548-Y_ST2,.tbl        Clostridioides_difficile_6622-CF_ST37.tbl         Clostridioides_difficile_7488-NonSp_ST103,.tbl
Clostridioides_difficile_5548-Y_ST2.tbl         Clostridioides_difficile_6623-NonSp_ST54,.tbl     Clostridioides_difficile_7488-NonSp_ST103.tbl
Clostridioides_difficile_5549-NonSp_ST21,.tbl   Clostridioides_difficile_6623-NonSp_ST54.tbl      Clostridioides_difficile_7491-DH_ST42,.tbl
Clostridioides_difficile_5549-NonSp_ST21.tbl    Clostridioides_difficile_6625-Y_ST2,.tbl          Clostridioides_difficile_7491-DH_ST42.tbl
Clostridioides_difficile_5550-AH_ST67,.tbl      Clostridioides_difficile_6625-Y_ST2.tbl           Clostridioides_difficile_7492-CF_ST37,.tbl
Clostridioides_difficile_5550-AH_ST67.tbl       Clostridioides_difficile_6626-Y_ST110,.tbl        Clostridioides_difficile_7492-CF_ST37.tbl
Clostridioides_difficile_5551-CF_ST37,.tbl      Clostridioides_difficile_6626-Y_ST110.tbl         Clostridioides_difficile_7494-E_ST6,.tbl
Clostridioides_difficile_5551-CF_ST37.tbl       Clostridioides_difficile_6628-AL_ST58,.tbl        Clostridioides_difficile_7494-E_ST6.tbl
Clostridioides_difficile_5552-DH_ST42,.tbl      Clostridioides_difficile_6628-AL_ST58.tbl         Clostridioides_difficile_7496-D_ST9,.tbl
Clostridioides_difficile_5552-DH_ST42.tbl       Clostridioides_difficile_6629-NonSp_novelST,.tbl  Clostridioides_difficile_7496-D_ST9.tbl
Clostridioides_difficile_5553-L_ST55,.tbl       Clostridioides_difficile_6629-NonSp_novelST.tbl   Clostridioides_difficile_7497-NonSp_novelST,.tbl
Clostridioides_difficile_5553-L_ST55.tbl        Clostridioides_difficile_6631-DH_ST42,.tbl        Clostridioides_difficile_7497-NonSp_novelST.tbl
Clostridioides_difficile_5555-DH_ST42,.tbl      Clostridioides_difficile_6631-DH_ST42.tbl         Clostridioides_difficile_7498-M_ST15,.tbl
Clostridioides_difficile_5555-DH_ST42.tbl       Clostridioides_difficile_6632-NonSp_ST3,.tbl      Clostridioides_difficile_7498-M_ST15.tbl
Clostridioides_difficile_5557-DH_ST42,.tbl      Clostridioides_difficile_6632-NonSp_ST3.tbl       Clostridioides_difficile_7499-CF_ST37,.tbl
Clostridioides_difficile_5557-DH_ST42.tbl       Clostridioides_difficile_6634-DH_ST42,.tbl        Clostridioides_difficile_7499-CF_ST37.tbl
Clostridioides_difficile_5558-AH_ST67,.tbl      Clostridioides_difficile_6634-DH_ST42.tbl         Clostridioides_difficile_7501-BI_ST1,.tbl
Clostridioides_difficile_5558-AH_ST67.tbl       Clostridioides_difficile_6635-BI_ST1,.tbl         Clostridioides_difficile_7501-BI_ST1.tbl
Clostridioides_difficile_5559-DH_ST42,.tbl      Clostridioides_difficile_6635-BI_ST1.tbl          Clostridioides_difficile_7503-Y_ST2,.tbl
Clostridioides_difficile_5559-DH_ST42.tbl       Clostridioides_difficile_6637-AL_ST58,.tbl        Clostridioides_difficile_7503-Y_ST2.tbl
Clostridioides_difficile_5573-CF_ST37,.tbl      Clostridioides_difficile_6637-AL_ST58.tbl         Clostridioides_difficile_7505-Y_ST2,.tbl
Clostridioides_difficile_5573-CF_ST37.tbl       Clostridioides_difficile_7450-DH_ST42,.tbl        Clostridioides_difficile_7505-Y_ST2.tbl
Clostridioides_difficile_5616-DH_ST42,.tbl      Clostridioides_difficile_7450-DH_ST42.tbl         Clostridioides_difficile_NCKUH-21,.tbl
Clostridioides_difficile_5616-DH_ST42.tbl       Clostridioides_difficile_7451-NonSp_ST2,.tbl      Clostridioides_difficile_RC10,.tbl
Clostridioides_difficile_5754,.tbl              Clostridioides_difficile_7451-NonSp_ST2.tbl       Clostridioides_difficile_RC10.tbl
Clostridioides_difficile_5754.tbl               Clostridioides_difficile_7452-Y_ST110,.tbl        Clostridioides_difficile_RF17,.tbl
Clostridioides_difficile_6583-AH_ST67,.tbl      Clostridioides_difficile_7452-Y_ST110.tbl         Clostridioides_difficile_RF17.tbl
Clostridioides_difficile_6583-AH_ST67.tbl       Clostridioides_difficile_7453-Y_ST2,.tbl          Peptoanaerobacter_stomatis_OBRC8,.tbl
Clostridioides_difficile_6584-AQ_ST26,.tbl      Clostridioides_difficile_7453-Y_ST2.tbl           Peptoanaerobacter_stomatis_OBRC8.tbl
Clostridioides_difficile_6584-AQ_ST26.tbl       Clostridioides_difficile_7455-A_ST43,.tbl         Romboutsia_maritimum_CCRI-22766,.tbl
Clostridioides_difficile_6588-G_ST8,.tbl        Clostridioides_difficile_7455-A_ST43.tbl          Romboutsia_maritimum_CCRI-22766.tbl
Clostridioides_difficile_6588-G_ST8.tbl         Clostridioides_difficile_7456-DH_ST42,.tbl        Tepidibacter_mesophilus_JCM.tbl
Clostridioides_difficile_6589-DH_ST42,.tbl      Clostridioides_difficile_7456-DH_ST42.tbl

```

Finally, you can convert these back to fasta with:

```
for f in *tbl; do TblToFasta $f > ${f/.tbl/.fa} && rm $f; done

```

This will only get those sequences whose naming scheme is as described above. But that's a start. You can then see what else you have and extract the rest manually or by finding a similar approach that matches those.

---
Tags: fasta, data-management
---


# Question
Title: I need to use my code with fasta format. Can anybody help?
I need some help. I have written a program that simulates protein cleavage but I need to modify it so that it produces fasta formatted output.

here is all of my sequence in .fas file

```
>jgi|Cormi1|6|CCM_00006m.01
MSSFTRAIRPALQGGRRVLLAQRLSTPRIAHPIAATALRSLSTTAPRRSVDLSELQPTPI
THLSEVETAMADAVSKFATEVVLPKAREMDEAEAMDPAVVEQLFEQGLMGIEIPEEYGGA
GMNFTAAIIAIEELARADPSVSVLVDVHNTLCNTALMKWGSDALKKKYLPRLATNTVASF
CLSEPVSGSDAFALATRATETPDGYTITGSKMWITNSMEAGFFIVFANLDPAKGYKGITA
FVVDKDTPGFSIAKKEKKLGIRASSTCVLNFDDVAVPKENLLGKPGEGYKYAIGILNEGR
IGIAAQMTGLALGAFENAVRYVWNDRKQFGTLVGEFQGMQHQIAQSYTEIAAARALVYNA
ARKKEAGEDFVRDAAMAKLYASQVAGRVAGLAVEWMGGMGFVREGLAEKYWRDSKIGAIY
EGTSNIQLNTIAKLLQKEYTA*
>jgi|Cormi1|8|CCM_00008m.01
MLAGSNARILRYIVVAGIILIFFYHFSGSRYDTLPRPVPAPDHAQAASQPSPPNDEATTA
GRPVDSTQKSIPADRDRPNSGFPKSVSASTDAAASTLADIDPGPRANATFVSLARNADVW
DIARSIRQVEDRFNRRFNYDWVFLNDKPFDDNFRNITSALVSGTARFGQIPVEHWSVPDW
IDEDKATKVREDMKERKIIYGDSMSYRHMCRFESGFFFQHPIMKNYDFYWRVEPSIELFC
DIAQDPFRFMKENGKKYSFTISLLEYVDTIPTLWDSAVKFIENYPQHVAADNSMEFLSDD
GGQTYNMCHFWSNFEIGSLEWMRSDAYVDFFNSLDQDGGFFYERWGDAPVHSIAAGLLLN
KTELHFFDEIAYFHMPFTHCPNDVQSRIDLKCHCNPKDNFDWSRYSCTSRYFEMMGMQKP
DGVD*
>jgi|Cormi1|10|CCM_00010m.01
MSTSPSDGTVDVDLPDAPTAPFPLHSDTHTANHVENDEDDDPITASYPVFLNPALPLGRR
LLVLQHPNRTDDHGARAPPSELRLKAESGMVEVDMPLDYNEAYDRDKGLRWGRTLQTSMA
AKNGGSHGLAGGFGFGAVQQRGAGAGAGGARKRNEVDAEDLIDWNEAVRQDKVLRTQTLG
GQYPDADEVQYMVGVFQGNDLHLTPVSSLVHLRPQLHHLDAATQQERQAAAKEAGAAAAS
SSGGATAGPGATGNPAARAIHMSIKTTADGETVTTETIADRLRFVQGEAWRKLRYTDENE
EAAWEVYNESLFLSASEARAPPPGSEKKKKTGSNEEELPYRSLQGAVATFASKWESGQLL
EAVSGMKKPPPVVVVPAAVKKAKAAAEAEAKAKDVEEKRKAAEAQMAATATAAAAAEETR
RPKRAPRGGGGLGTAKRPSRAK*
>jgi|Cormi1|11|CCM_00011m.01
MASATTPILPAEGKRNVLITSALPYSNNLPHLGNVIGSVLSADLYARYCRARGYQTLYVC
GTDQYGTATETKALEEGTDPASLCAKYHALHKSIYDWFRIDFDIFDKTPTQHQTDIVQDI
FTKLWKNGYIYEDETIQPFCPVESHNTFLADRFVEGECSICGDKGARGDQCDACGSILDP
LKPEETGAVDGDAAIEVKATGWLVNPRCKVDGTAPEARKTKHLYLKLDDISEELKAWYKT
TADNGWSTNAKAITQSWIDKGLKPRAISRDLKWGVPIPAVEGLNVEDYARKVFYVWFDAC
IGYISITKTYTDGDNRDGKNWEKWWKNPENVELMQFQGKDNVGFHSIVFPSSLLGTKENW
TKVTRISATEYLNYEQGKFSKSRGIGVFGNNAKDTGIDPDIWRFYLLSRRPELADSEFKW
EEFIDVNNNELLKNVGNFIQRTLKFCVAKMDSTVPAATVSNDVIESHKQSANQILASYIE
HQEANKMRIGISDILTFSALGNKFLQDNGLNNKLLAEEPELCNTVINVAINHVHLLAGLL
YPYMPTTAESIFQQLGQEPNPFIPDVFGFTAVPVGHKIGTPGPLFSVIPAAKLDEWRDAY
GGAELKKQKAIEAEKAAAKKAKKLADKEKKKAKKAETADLKTEGASTAADISKLSLEEKK
*
>jgi|Cormi1|15|CCM_00015m.01
MAEEEVRSTVALVVGALLPALMTVFFNGLPGAVGGSNLPRWGVRLPILLTCAAQLWFMLV
SRVVPEPYLDEIFHIPQAQKYCEGRFFDWDDKITTPPGLYFFSVAVHRIASALRIPQLSV
CDPFSLRAVNFAGVLVVSYLALWCRQALEARQHEAISSPTPARVRAFSQYAVHTAVNIGL
FPLLFFFGGLYYTDVLSTGVVLAAFLNHLRRVGVPRSSVWSDLGTVALGVAALGMRQTNV
FWIVVFMGGLEAVHAVKTLRPARVDQPVMTTLAAQCRYFAWRYSLGDIHDLPVHKAYPDD
MIFTALSIGIAALCNPLRVVRQIWPYLAVLASFAAFVQWNGSVVLGDKSNHVATIHLAQM
LYIWPFFAFFSLPLLLPSGLSFLNMVVSFFRTHARASPTPDAAADGLLRATLRFFHSKVL
WPAYLAATAVLSLAVVRFNTIIHPFTLADNRHYMFYVFRYTIRRGGAVRLALALPYTVCR
WVVWDALGGHNSWSLGNDKPTAFVSHPFWTPPSAKVRRQTSTSYPRPPERASTAESRQAS
EAALKARLEQDSLYLSQQSVSTSTGLIFLLATALSLVTAPLVEPRYFILPWVVWRLMVPA
WRLQDHLRGGILEGSPAAKAVAHWSASYDLRLFLESAWHVVINLATGYVFLYKPYVWRAE
DGTVLEGGNLQRFMW*
>jgi|Cormi1|20|CCM_00020m.01
MTNTIVIVGQVLFPLLTRTTNCIANERRPPSSGGVSGLTSAYLLSKSKGNTVTLVAKHMP
GDYDIEYASPIAGANVLPMNPDADSRWERRTWPEFVWLCRNVPEAGFHFKNALVLRQHKD
QEADPSGGPIDPLFDRDPWYRHLTPDYRELAASELPPGQDSGCSFTSMCLNTPMYLFWLV
GQCLRQGVVFRRAVLHDLTEARRLAYSGRPADVLINATGLGARTLGGVRDAGMLPARGQI
VLVRNDGSSMLNTSAVDDGSGEILYTMQRAAGGGTILGGTYDLGNWESQPCPNVATRIMA
RAVRLDPALAGGKGVAGLDVIRHAVGLRPYRKGGVRVERETLRDEVDGASVHVVHNYGHA
GWGFQGSFGCAEKVVELVNEIRAANGEDLSKEPTLFNFEAAAASKL*
>jgi|Cormi1|21|CCM_00021m.01
MEATNSPARAAREPAVHVRALTFVVDVACGRHALSPLIPVALWLADAVLTCLIIWKIPYT
EIDWVAYMEQITQFVHGERDYTKIKGGTGPLVYPAAHVYTYTGLYYLTGEGENVLLAQQL
FAVLYMVTLGLVMLCYWKAKVPPYIFPLLILSKRLHSVFVLRCFNDCFAVFFLWLAIYSF
QRRLWALGSMAYTWGLGIKMSLLLVLPAVGAVLLHACGFSGALSLAWLMAQVQIVIALPF
IATNTNGYLARAFELSRQFQFEWTVNWRMLGEATFRRREWALLLLGLHAATLALFLATRW
MRPAGKTLRRMVPAWLRGRAALTPDEEIRVSRRVTPDFVLTTMLSANLVGLLFARSLHYQ
FYAYLAWSTPFLLWRALPGGAAVLTYPLWLLQEWAWNVFPSTSESSTAAVAVMATTVLVV
YFGTRNDGLTIKKPAAKKQ*
>jgi|Cormi1|32|CCM_00032m.01
MRLIAAQMLRGTAQRRGFATTPRRLDNYAFVGLGQMGYQMAKNLQSKLQPSDKISLYDIN
TDAMKRLQVEMKAASNGAAVELAASAHDAAKEADTVITVLPEPQHVQGVYKSILTPSLPR
RDRVFIDCSTIDPSTSRAVAASVAAAQQGAFVDAPMSGGVVGATAGTLTFMLGAPASLLP
RLEPALLTMGRKVLHCGDQGAGLSAKLANNYLLALNNIATAEAMNLGIRWGLDPKKLAGV
LNVSTGRCWPSEVNNPVKGVVETAPAGRDYSGGFGIALMKKDLRLAMVAAAEAGADMKLA
DAAFKVYEEASEREDCKGRDFSVVYRYLGGKE*
>jgi|Cormi1|35|CCM_00035m.01
MPGLPTSVDLDECISRLYKKELLAESVIEAICAKTKELLMRESNVVHVQAPVTVVGDIHG
QFYDLIEIFRIGGYCPDTNYLFLGDYVDRGMFSVETISLLVCLKLRYPERVHLIRGNHES
RGVTQSYGFYTECSRKYGNANVWHYFTDMFDFMTLSVVIDNQIFCVHGGLSPSIHSIDQI
KIIDRFREIPHEGPMADLVWSDPDPDRDEFSLSPRGAGYTFGAQVVKKFLAVNNMSHILR
AHQLCQEGFQVLYEDRLSTVWSAPNYCYRCGNMASVLEVNSDGERFFNVFAAAPENDQVK
DMQPGIDKAGDGNALPDYFL*
>jgi|Cormi1|37|CCM_00037m.01
MRYAFAYAALVSTASAHGLVTMIKGDNGVMMPGLGVTDGTPRDCSSNGCGSQADTSIIRN
RDMQQNPLGKTQGTGPVDAAVNVAAFMGTSKNAPKNNGASAGVGQEDDLSGLPAGAAKRH
EERKNMVRGLFGDLLGGAAGGAGAAGGAAGGAAGGLGGLLGGAGGAGGAAGGAAGGLGGL
LGGLGGGGGGTKSNLAAESMVADTAGMGSTQGMPTSDSTGTVSLVYRQINQDGAGPLTAA
VDGTSGGTDASAFQQAEVTKDVPGLGIGGLSLATNTDFPLQVKMPAGMTCEGTVGGVKNV
CIVRVRNSAAAGPFGGSAAFTQNATARKRAIAYRLRKRMEIAREAVPQIDYTKADAAEEL
EDEEDE*
>jgi|Cormi1|65|CCM_00065m.01
MRWHVYWGVVQLVELSNVDKRVAVACHVHTLALPDISLSHHPASSATVTPVIIVRLNRPE
SRHAYTDTMGASLVTAYGLLSADPRVRAIVLTSADASNKFFCAGMDFNVRAQRVEDPLQY
RDSGGTVSIAMHRCQKPIVAAVNGAAVGIGMTMTLPATVRVASTKAKCGFVFAQRGFCME
ACSSFFLPRLIGTAKAMHLVSTGAVYPATHRLFDDLFSEVVEPEDVLPTALKIAEGMAQN
CSIVSTQVMKDMMYSGPKSPEEAHLLESKIFYSLTCPDTDSAEGVQSFLQKRPPKFTGTM
ENQRPFGWPWWSNLDVRPKI*
>jgi|Cormi1|71|CCM_00071m.01
MSRIEKLSISGVRSFSPACREAIQFNTPLTLIVGYNGSGKTTIIECLKYATTGELPPNSK
GGAFIHDPKLCGEKEVLAQVKLQFRSINDRQHVATRSLQLTVKKTTRSQKTLDCSLVVLN
NGEKTTTSTRQAQLDEMIPERLGVSPAILDAVIFCHQDESLWPLSEPSVLKKRFDEIFEA
LKYTKAIENLKVLRKKQVEELGKLQNDEFHNKGNKDRGDRAEKRMTGLQAEIEAAREKCE
AITAHMQETQDKIKQKHEQANSFLQIVQNLGNKREQLDYRQDAVDELRQTIDELPEDDRS
LERQLVQYEESMRRQSEEAQRNKTQYGDLQQDLVKSRKSLSTKLAEQGKHQSDKEKHERQ
LKMRMEMVQEAAQHHGFRGYEDDLTDALIKQFNDKIQKLFAEKKRELERLQKENAAELDK
ASATISEFEGQKAARTQDRVSAKQRQGAIEKRTIVLQNEASSIDVDEGLKAILDGQLGEL
ESKFSQVHREFETAAWDNQISDENAKLWQLESESEKLGRELIECTRLASERAQLDFRRKE
LSDRKLKLTSLTETWRSKLDKVIGGDWSPEDLEAKFQALLAAQNKALQEARQQQEQTRQK
QARVDFKLKTVKETCEQKQTDQIECKKKVLQALQNVRDSAVIEDYTEEVSATEQQIETLR
DDLSLIDALKDYYTKCKRALDAKKTCLLCDRHFDDSQAASIDRLRVKIDRHLNSKTKIEA
EQDLSEASSTLQSLVAVRSQYDTYERLGKELPGLRDERKSIEADHESTERLLEEHDAVVA
AEEDKKRDMNDMSKTVSNITQTLREIKDSEGQVERIASQQQSVGTGRSADEIHDLQGTVS
EQMRASKNRVAKLSNDRQRARDQISTLELEKSELRNKISLTLGQLDKKQDLLGQIQSLKE
EQTHQREAIQRADEELEKVEPRITEARSIRDETLKRGRAKEQVVTESRDSVATSINEIKM
IESDIRDYEDRGGASNLASNQRAIATLEQNITTTETEISELTQRTNKLKQDIDNGDRTKK
NINDNLNYRKHLRMLETLRGEINELEDRNAHEDYERLQGEARALENQSNRLFAERGSVMG
SMKTKDEDLDRLLKEWEMDYKDAKQKYRESHIRVETTKAAIEDLAQCGAAVDKAVMQFHA
MKMAEVNRIAGELWQSTYQGTDIDTILIRSDNEGTTGRRNYNYRLCMVKQDAEMDMRGRC
SAGQKVLASIIVRLALAESFGVNCGLIALDEPTTNLDRDNIKSLAESLHNIIRTRQAQSN
FQLIVITHDEDFLRHMRCSDFCDSFFRVKRDDKQNSVISRESITKIF*
>jgi|Cormi1|75|CCM_00075m.01
MPSSTPLPPPKGFLSYPDVEEIISPRYEPRKDANPVLSGLPLVIASSLLTHSSFLAKIFY
NNAGLQKIKDLPFLDDVPYTFHPLVIPLGDNGPMVELGEAQLNPKNRGDTHHYTASDYHA
MYKSGQVTPLQVTETLLSITSKYAKEASIYQDAWADSHGADQLALDAAKASTERWAAGKA
LGVLDGVPIGVKDDLSVKGYIDHDGMAYKKGVPFFKSAEETVWPVQALQAAGAVVIGKNR
MHELGSDTNGLNPHQGTPTNHLNNQYYPGGSSSGAGSSLGAGVIPIAVGTDAGGSVRFPA
VYNGVYGLKPTHHRTLEMNSTMCVAGPMAANVADLTIAYRIMSQPDPTTSIQSRFAVSQP
PAAGTNKYIGVDRAWWGASDPRVIGACNKALEYFTSQGYELVDITIPYIDEAQTAHGGVT
IAEMTEAARRRGGDKTSWTTLIGPVNKLVLALGLETNAADLLKFNSMRTLIMRHMAWLFQ
KYPGLLILTPTAPFLGWEKKPAHAKYGVSDGDKTFLSMIYVFLANIAGLPALTAPIAYVD
PDQGEGKVCIGMMATGEWGAEEALLAWAADAETYLHEVYPEGRRRPKTWVDVLGLAAKQG
TV*
>jgi|Cormi1|77|CCM_00077m.01
MLSSTITRSRPTATMAKLAGGMLVAGLATKYTIGVAHAESPAAPPRSQVFGRGPALVSLV
LEDAEMVNHNTRRLRFRFQDKEAASGLPLTSALLTYSWPKDSWLPVIRPYTPITTADETG
HLDLLVKQYPNGKSSTHLHSLKPGHTLFFVAPLSGHQWQPNSYKHVTLIAGGAGITPIYQ
LAQGILNNPDDKTAVTVVFGVNSDRDVLFKKEFDEYEQRFPGRFRAVYTVSSPVPGSPHR
KGHVNEALLKEVTPAPGPGNMVFVSGPPPMETALVGSRREPGILEKYGYRKEQVYKF*
>jgi|Cormi1|78|CCM_00078m.01
MAADASFDNLARQFSHLPLIDKLPNCFPDINPVDIYRAHITSILSHITGADKSLIYPVLQ
WTATLDKGDLMLAIPALRIKGKPDDLGKQWVEKWPESPLVEKPVHNGAFISFFFKPEPMA
KTVIPMALQHKQDFGSNPANGLRDQNDPSKGRQKIVIEFSSPNIAKPFHAGHLRSTIIGG
FLANLYRSAGWDVATVNYLGDWGKQYGVLALGFEKYGSEEELVKDPINHLFTVYVKISKE
VAEEKEQADKLKAEGKEEEAQKILDNGVDEKARSYFKKMTEGDEQAIALWKRFRDYSIER
YKQTYARLNIHFDEYSGESKVSEDEMQQAAAKLKEMGLAHEDRGALLMDFTKSVPGKEGK
TLGKAILRKKDGTALYLTRDISELINRHERYKFDKMIYVIASEQDLHVKQFFKIVRMMGR
DDIADKCQHVNFGLVLGMSTRKGTVKFLNDILRDVGDHMHDVMKKNEDKYSQVEDPEQIA
DKLGISSVMVQDMNGKRINNYTFDLSRMTSFEGDTGPYLQYAHARLCSIYRKAELSDEAI
AAADLSLLTEPMAVNLLRVIAQWPDVVQNVLKTLEPTTILTYLFKLTHTLSTSYTHLRIV
GSEPEVMKARLALYNSTRVVLYNGMAILGLTPIERM*
>jgi|Cormi1|79|CCM_00079m.01
MRSRNASPHIGPWRNPARRVRHGTTRNYFDSGCSVWMSSPSRLSAPVPPASPLQKNKFSH
IAIYEQRAATGGVWNTSHVSQEPGFSIPRTTPTSEHEYVVVEDGHDAGARVELVSPVYDN
LETNIPHGLMSYTDLEFPAGTPLFPEHGTVLAYIQQYGRDVEHLVTFETQVRDVRKTTAE
DGRLGWTVETKDLKSGVVSAKSYDAVVAASGHYSDPFVPDIPGIAEFEAVQPGAILHSKF
YRRPEQFTGKKVVVVGNAASGIDISHQIATSAQLPVLISEKDVPGAPPAASNNWSRHVGQ
IVELLPVTRSVRFSSGRIEDGVDAVIFCTGYHYSFPFLGSLAPGITAPNGTYADHLWEHM
LYAPDPTLALLVIPKRIVPFPFAEAQMAVIARIWAGRLNVPNQSERDAWVRRCGAESAPG
ARHTFGYPQDVDYINRLHALSGQARPDAERGLENDGAGKRSPYWGEELSWVRSQIFAIKV
ASRALGERRKECKTLADLGFDFQAAK*
>jgi|Cormi1|87|CCM_00087m.01
MFLLATLFALSALPWTAPAMRSDHLAKLRQETVDMFYHGYNNYMLHAFPEDELRPVSCQP
LTRNPNKPADIGLNDVLGNYSLTLIDSLSTLAILAGGPANEKYTGSQALSDFQDSVAQFV
SYYGDGRQGSSGQGLRSRGFDLDSKVQVFETVIRGVGGLLSAHLFAVGDLPIRGYDPKPD
SKWASEDPLETAPISWHGGFRYDGQLLRLALDLSERLLPAFYTKTGIPYPRVNLRSGIPF
YINSPLHQDSEDTEPVVGNAEITETCSAGAGSLTLEFTVLSRLSGDPRFEQAAKRAFWEV
WERRSDIGLIGNGIDAERGVWIGPHSGIGAGMDSFFEYALKSHILLSGQETPNMTRSHRQ
STTGWLDPNSLHPPLPPEMHSSDAFLEAWHQAHASVKRYIYTDRNHFPYYSNNHRATGQP
YTMWIDSLGAFYPGLLALAGEVDEAVEANLVYTALWTRYGAIPERWSVRESNVEPGIAWW
PGRPEFIESTYHIYRATADPWYLRVGEMVLHDIQRRCWVACGWAGLENVLSGEKQDRMES
FFLGETTKYLYLLFDPEHPLNNLDAAYVFTTEGHPLIIPRRARGSSTPRRRASRTSRQEE
KDVAVYSYYDTRFTNSCPAARTPTEPLAGSNTAARRNLFDVSRFTNLYNTPNIHGPLESF
KTRDEHKGIITQYRATSNYSVFPWTLPPTMLPHNGTCPVPSNRILSWIEFPTVDPAPSLF
ATRGGGGDSNTPLRWYHNVGPTVERVDGLKLQLEQEFNEAAGTDVWMVTHVAGKTVGRHE
NVIFHADHVRHFRDDAFSVLRRKDMVEVVLLVQPEPTPQASRELPPAVMPLPVPGALDEV
SSLDAAPDNTSSTRRPSVRASGQPPATDPDSLFKSILRAVSSALEGSATASTPSPSPLRG
SGSGSGSRDDSSNSDDDVGVASIHSWLANTATGAGAFPLPSVRDALLRNSPDFSAADPAA
SLPWRHVYMAGSACDGPLPEDAPRRHQVIVMRRGGCPFGDKLAHVPNFVPDEASLQLVVV
VDEEGGGHVDDAPRPLLATEQRTPRGTPRLHGIPMVLLPASKGTYELFGEAAAVGLRRKY
VVKSQGYVIENAIVM*
>jgi|Cormi1|88|CCM_00088m.01
MSAGHRVARVHLKEKNTAQDASGRQPSSSSSATAQLASAKVVLGNNLRRFPDFPIPGINF
VDIMPLFANPEAHATLLNALELQITQTFVNTKPDVIVGLDARGFLFGPGLALRLGVPFAA
VRKQGKLPGPCVTAEYVKEYGKDFFQMQEDAIREGQKVLVVDDIIATGGSAKAASDLITQ
LKGHVMGYLFILEIPGLNGRDKLGAAPTTILLQDE*
>jgi|Cormi1|96|CCM_00096m.01
MHLVPKELDKLVITQLGLLAQRRLARGVKLNHAEATALIANNIHELIRDGNHTVADLMAL
GATMLGRRHVIPSVCHTLKEIQVEGTFPTGTYLVTVHNPISTNDGDLARALYGSFLPVPR
DADRLFPAALDPADFAPARQPGAVVAVKGDKITLNADRRRIALRVVSEGDRPIQVGSHYH
FIETNPQLVFDRARAYGFRLDIPAGTSYRFEPGDAKTVTLVEIAGHQVIRGGNGLATGGV
ERWRVDAIVERLQKAGYAHAPEPEVNGVVAGRRRAAAARVPRPYRMDRTAYTVMFGPTTG
DRVRLSSTDLWVKVERDYTTYGDECKFGGGKTLREGMGQASGRPDADSLDLVITNALIID
YTGIVKADIGVKNGIIVGIGKAGNPDVMEGVSPNMIVGSCTDVIAGENKIITAGAIDTHI
HLICPQQAQEALASGVTTFLGGGTGPSAGSNATTCTPGKHLMKQMLQACDALPVNVGITA
KGSDSAPEALREQVAAGACGLKIHEDWGATPSAIDTCLSVCDEMDVQCLIHTDTLNESGF
VESTIASFKNRTIHTYHTEGAGGGHAPDIISVVERANVLPSSTNPTRPFTRNTLDEHLDM
LMVCHHLSKDIPEDVAFAESRIRAETIAAEDVLHDTGAISMMSSDSQAMGRCGEVVLRTW
NTAHKNKLQRGCLPEDDGTGADNFRVKRYVSKYTINPALAQGFSHLVGSVEVGKLADLVV
WDPAWFGTKPHQVIKSGLIAWSQMGDPNASIPTVQPVIGRPMFASFVPAASVVFVSQASV
ASGAVASYGLRKRVEAVRGCRTVGKADMRFNDAMPKMRVDPESYTVEADGEVCTAEPAET
LPLTQSWYGILDNQNAPA*
>jgi|Cormi1|100|CCM_00100m.01
MAARLSAAALAPLRPISSSVCAARQLSQSLLQSQQTRGVAQTYLKKVAEGEERWKARAER
IQNGEQRHVWDILDERGYIKDVAGNTDRIKEVMRIKRIGAYVGIDPTADSLHVGHMLPLM
PLFWLWFHGYPAVTLIGGSTARIGDPTGRLQSREHITNSDISKNITKIHYQLSKLWHNVV
QMRIKYGYEDDWAAKRSLLNNNMWLQGLTLYDFIKRLARDTRIGPMLSRDTAKRKLTDGD
GMSLGEFMYPLLQGWDFWHMYNKIGVQMQIGGSDQYGNIVSGIESLKTIRETEEAPYARM
ESGWQHDPIGFTVPLLTDSAGVKFGKSAGNAIWLDEFKTSAFELYGYFMRRTDDEVEKLL
KLFTFLPMPKIQSIMEEHNADPAKRIAQHNLAFEFVSLIHGSQKALEEAQQHSFRFGGEL
PSIVKEPAASSGIINANTAPRSDIKLPRSVMKLSPAKLLFATGLAPSAAEGQRLVKMQGA
YVAAQPGQKRGLVPGNLNWTPMKMWYPEETSKFLLDDRMLILRKGKHNVRIVELVDDEEW
KASGEFYPGEPYTGEVRRRKQEIIAQAAEDGKTLTMTQVNKMLKKKGRQEGLRVANNPDI
ELPTKTEVRERVRLTHMHKKEERRTRGDEHEW*

```

My current code:

```
import re

myprotein = ""
my_protein = []
for protein in myprotein:
    myprotein = protein.rstrip('\n')
    my_protein.append(myprotein)
    my_pro = (''.join(my_protein))

   #cleaves sequence
peptides = re.sub(r'(?<=[^HRK][^P][^R])(?=[FLWY][^P])|(?<=[^HRK][^P][FLWY])(?=.[^P]?)','\n', my_pro)

print (peptides)

```

the output is like

```
MSS
F
TRAIRPAL
QGGRRVLL
AQRL
STPRIAHPIAATA
L
RS
LSTTAPRRSVD
L
SE
L
QPTPITH
L
SEVETAMADAVSK
F
ATEVVL
PKAREMDEAEAMDPAVVEQ
L
F
EQG
L
MGIEIPEE
Y
GGAGMN
F
TAAIIAIEE
L
ARADPSVSV
L
VDVHNTL
CNTA
L
MK
W
GSDA
L
KKKYLPRL
ATNTVAS
F
C
L
SEPVSGSDA
F
A
L
ATRATETPDG
Y
TITGSKM
WITNSMEAG
F
F
IV
F
AN
L
DPAKG
YKGITA
F
VVDKDTPGF
SIAKKEKK
LGIRASSTCV
L
N
F
DDVAVPKENL
L
GKPGEG
Y
K
Y
AIGI
L
NEGRIGIAAQMTG
L
A
L
GA
F
ENAVRY
VW
NDRKQFGT
L
VGE
F
QGMQHQIAQS
Y
TEIAAARA
LV
Y
NAARKKEAGED
F
VRDAAMAK
L
YASQVAGRVAG
L
AVE
W
MGGMG
F
VREGL
AEK
Y
WRDSKIGAI
Y
EGTSNIQ
L
NTIAK
L
LQKE
YTA

```

I need every sequence on a separate line, starting with the sequence header like:

```
>jgi|Cormi1|6|CCM_00006m.01 frag.1
MSS
>jgi|Cormi1|6|CCM_00006m.01 frag.2
F
>jgi|Cormi1|6|CCM_00006m.01 frag.3
TRAIRPAL

```

How could I make this work?

```
import sys
import re

sequence = open(input('Enter input filename: '),'r')
def digest(sequence):
    return re.sub(r'(?<=[^HRK][^P][^R])(?=[FLWY][^P])|(?<=[^HRK][^P][FLWY])(?=.[^P]?)','\n', sequence)

seq = ""
header = None
for line in sys.stdin:
    if not line.startswith('>'):
        seq = seq + line.rstrip()
    else:
        if header:
            for idx, fragment in enumerate(digest(seq).split('\n')):
                sys.stdout.write("%s frag.%d\n%s\n" % (header, idx + 1, fragment))
            seq = ""
        header = line.rstrip()

# final pass
for idx, fragment in enumerate(digest(seq).split('\n')):
    sys.stdout.write("%s frag.%d\n%s\n" % (header, idx + 1, fragment))

```

# Answer
One way is to store the header as a separate string from your protein residues. Use a list comprehension to build a list of formatted header-fragment pairs, and then print them to standard output.

```
#!/usr/bin/env python

import sys
import re

myprotein = ">jgi|Cormi1|6|CCM_00006m.01\nMSSFTRAIRPALQGGRRVLLAQRLSTPRIAHPIAATALRSLSTTAPRRSVDLSELQPTPI\nTHLSEVETAMADAVSKFATEVVLPKAREMDEAEAMDPAVVEQLFEQGLMGIEIPEEYGGA\nGMNFTAAIIAIEELARADPSVSVLVDVHNTLCNTALMKWGSDALKKKYLPRLATNTVASF\nCLSEPVSGSDAFALATRATETPDGYTITGSKMWITNSMEAGFFIVFANLDPAKGYKGITA\nFVVDKDTPGFSIAKKEKKLGIRASSTCVLNFDDVAVPKENLLGKPGEGYKYAIGILNEGR\nIGIAAQMTGLALGAFENAVRYVWNDRKQFGTLVGEFQGMQHQIAQSYTEIAAARALVYNA\n"

# parse record
my_protein = []
for line in myprotein.split('\n'):
    if not line.startswith('>'):
        my_protein.append(line)
    else:
        header = line

# build sequence
my_pro = ''.join(my_protein)

# cleave sequence
peptides = re.sub(r'(?<=[^HRK][^P][^R])(?=[FLWY][^P])|(?<=[^HRK][^P][FLWY])(?=.[^P]?)','\n', my_pro)

# build a formatted list
records = ["%s frag.%d\n%s\n" % (header, idx + 1, fragment) for idx, fragment in enumerate(peptides.split('\n'))]

# print list
for record in records:
    sys.stdout.write("%s" % (record))

```

This gives:

```
>jgi|Cormi1|6|CCM_00006m.01 frag.1
MSS
>jgi|Cormi1|6|CCM_00006m.01 frag.2
F
>jgi|Cormi1|6|CCM_00006m.01 frag.3
TRAIRPAL
...
>jgi|Cormi1|6|CCM_00006m.01 frag.89
LV
>jgi|Cormi1|6|CCM_00006m.01 frag.90
Y
>jgi|Cormi1|6|CCM_00006m.01 frag.91
NA

```

---

If you have multiple records in one file called `seqs.fa`, then you could use standard input to this script (*e.g.*, named `digest.py`) to process them all at once:

```
#!/usr/bin/env python

import sys
import re

def digest(sequence):
    return re.sub(r'(?<=[^HRK][^P][^R])(?=[FLWY][^P])|(?<=[^HRK][^P][FLWY])(?=.[^P]?)','\n', sequence)

seq = ""
header = None
for line in sys.stdin:
    if not line.startswith('>'):
        seq = seq + line.rstrip()
    else:
        if header:
            for idx, fragment in enumerate(digest(seq).split('\n')):
                sys.stdout.write("%s frag.%d\n%s\n" % (header, idx + 1, fragment))
            seq = ""
        header = line.rstrip()

# final pass
for idx, fragment in enumerate(digest(seq).split('\n')):
    sys.stdout.write("%s frag.%d\n%s\n" % (header, idx + 1, fragment))

```

Creating a function called `digest` lets you change the digestion rules in one place. It reduces code bloat and also reduces errors from changing things in one place and not another.

You could use this script like so:

```
$ ./digest.py < seqs.fa > seqs.digested.fa

```

The file `seqs.digested.fa` would look like:

```
$ cat seqs.digested.fa
>jgi|Cormi1|6|CCM_00006m.01 frag.1
MSS
>jgi|Cormi1|6|CCM_00006m.01 frag.2
F
>jgi|Cormi1|6|CCM_00006m.01 frag.3
TRAIRPAL
...
>jgi|Cormi1|100|CCM_00100m.01 frag.143
PTKTEVRERVRL
>jgi|Cormi1|100|CCM_00100m.01 frag.144
THMHKKEERRTRGDEHE
>jgi|Cormi1|100|CCM_00100m.01 frag.145
W*

```
> 5 votes

---
Tags: python, fasta
---


# Question
Title: Read sorted and indexed BAM files faster in C++?
I have some sorted and indexed alignment large BAM files. I'd like to read all the reads, and apply custom operation on each of them. The orders aren't important. I'm using `htslib` in C++ to do the reading. I have a single powerful machine.

The problem is ... reading a large BAM file is very slow even in C++.

**Q:** Is there any trick to speed up the reading? Reading BAM files multithreaded (I have many cores)?

# Answer
In this answer I want to show you some benchmarks that compares three different serial ways of reading data in C++, the third one being the fastest.

In the first method I use an `std::istreambuf_iterator`, in the second I read the file line by line into a `std::vector` and in the third one I use C flavoured operations.

```
/*
compile with g++ -std=c++11 -O2 reading_files.cpp -o reading_files
test the program with the following
dd if=/dev/random of=./dummy5GiB bs=1024 count=$[1024*1024*5]
 for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
*/

#include <fstream>
#include <iostream>
#include <vector>
#include <cstring> //for strcmp
#include <string>
#include <algorithm> //for min

using namespace std;

int main(int argc, char * argv[])
{
    if(argc != 3) {
        cerr << "usage: reading_files <path> <number 0 or 1 or 2>\n";
        return 1;
    }

    ifstream file(argv[1], ios::in|ios::binary);
    if(!file.is_open()) {
        cerr << argv[1] << " file not opened!\n";
        return 1;
    }
    cout << "---------------------------\n";
    if(!strcmp(argv[2],"0")) {
        cout << "method 0\n";
        vector<char> content((std::istreambuf_iterator<char>(file)), std::istreambuf_iterator<char>());
        cout << content.size() << " bytes read\n";
    } else if(!strcmp(argv[2],"1")) {
        cout << "method 1\n";
        string line;
        vector<string> content;
        while(getline(file,line)) {
            content.push_back(line);
        }
        cout << content.size() << " lines read\n";
    } else if(!strcmp(argv[2],"2")) {
        cout << "method 2\n";
        file.seekg(0, std::ios::end);
        unsigned long long file_size = file.tellg();
        cout << "file size is " << file_size << "\n";
        //maximum 5 GiB of memory, so to fit my RAM
        unsigned long long max_block_size = 1024ULL * 1024ULL * 1024ULL * 5;
        file.clear();
        file.seekg(0, std::ios::beg);

        unsigned long long block_size = min(max_block_size,file_size);
        char * buffer = (char*)malloc(block_size*sizeof(char));
        unsigned long long processed = 0;

        while(!file.eof()) {
            file.read(buffer, block_size);
            cout << "read a block of " << min(file_size - processed, block_size) << " bytes\n";
            //do work here with the block in memory (if the block size is > 0)
            processed += block_size;
        }

        free(buffer);
    }
    file.close();
    return 0;
}

```

In the third method I read the file in blocks of 5GiB, so to fit the free RAM of my machine.

Here is the benchmark of the program on a 5 GiB randomly generated file.

```
bash-3.2$ for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
---------------------------
method 0
5368709120 bytes read

real    0m35.194s
user    0m19.665s
sys 0m12.999s
---------------------------
method 1
20973135 lines read

real    0m40.478s
user    0m34.046s
sys 0m5.870s
---------------------------
method 2
file size is 5368709120
read a block of 5368709120 bytes
read a block of 0 bytes

real    0m4.413s
user    0m1.757s
sys 0m2.617s

```

Here is the benchmark on a real 32GiB .bam file using only the third method, here the splitting of the file read plays a fundamental role. I can't tell you how big because my hard disk memory gets saturated before the end of the execution and so I have to kill the program.

```
time ./reading_files ./HG00252.mapped.ILLUMINA.bwa.GBR.low_coverage.20130415.bam 2
---------------------------
method 2
file size is 34316054058
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 2103799338 bytes

real    0m50.656s
user    0m9.135s
sys 0m35.728s

```

The key point is to minimise the file operations reading large block of data at once (as long as they fit the memory) and operate on them.

Using `mmap` could lead to an improvement, but I don't think it would be as large as using the third method instead of the first two.
> 4 votes

# Answer
If you need to iterate through the entire BAM file, decompression is usually the main bottleneck. Three things you can do:

* Install libdeflate (https://github.com/ebiggers/libdeflate ) and then the development version of htslib, and make sure the latter is compiled to use libdeflate. This roughly doubles decompression efficiency in my experience.
* Call htslib's bgzf\_mt() function right after opening the file for reading, to request multiple decompressor threads. I generally go with 2-6, depending on what I'm actually doing with the reads.
* Have a separate read/decompress-ahead thread perform the bam\_read1() calls. Ideally, when your main thread is processing read #X, background threads are already working on decompressing reads e.g. (X+100) through (X+200).

Also, if your workload is "embarrassingly parallel", you may also benefit from launching multiple processes simultaneously, configuring each to handle a similarly-sized chunk of the BAM.
> 3 votes

# Answer
In both deepTools and MethylDackel we implement multithreading by having individual threads process separate genomic regions. At some point you hit an IO limit, but unless your processing is trivial then IO won't typically be limiting with single-threaded applications. Since you already have sorted/index BAM files, you can take advantage of the same strategy.
> 2 votes

---
Tags: sam, c++
---


# Question
Title: Find specific functional genes from a list?
I have compiled a list of genes that are distinct from 2 different conditions. From this list (of about 600 human genes), I want to find out if any of the genes in this list have acetyl transferase function i.e. are acetylases.

Is there an online platform where I could say put in this list of genes and the output I can select genes which have this property/function?

I.e some sort of gene annotation for functions for a list of human genes?

# Answer
BioMart. Filter by your list of genes and a GO term linked to acetyl transferase activity.
> 2 votes

---
Tags: sequence-annotation, functional-annotation
---


# Question
Title: How does htslib/samtools access optional BAM fields?
I am confused regarding how various software libraries deal with optional fields in a BAM:

Based upon the BAM specification, there are 11 mandatory fields to a BAM:

```
QNAME, FLAG, RNAME, POS, MAPQ, CIGAR, RNEXT, PNEXT, TLEN, SEQ, QUAL

```

The BAM files of today normally have many optional fields as well, which must be of the format `TAG:TYPE:VALUE`, whereby `TAG` is a two-character string that matches `/[A-Za-z][A-Za-z0-9]/` and `TYPE` is a single case-sensitive letter which defines the format of `VALUE`. Each `TAG` can only appear once in one alignment line.

Certain BAMs may have no optional fields, while others have several.

My understanding is that optional tags may be in a different order from row to row, and then sometimes they simply do not exist.

My questions:

1. How does htslib/samtools (e.g. `samtools view`) access these optional fields, which differ in number from BAM to BAM (possible row to row)?

In `htslib/sam.h`, a BAM alignment struct is defined as follows:

```
/*! @typedef
 @abstract Structure for one alignment.
 @field  core       core information about the alignment
 @field  l_data     current length of bam1_t::data
 @field  m_data     maximum length of bam1_t::data
 @field  data       all variable-length data, concatenated; structure: qname-cigar-seq-qual-aux
typedef struct {
    bam1_core_t core;
    int l_data;
    uint32_t m_data;
    uint8_t *data;
#ifndef BAM_NO_ID
    uint64_t id;
#endif
} bam1_t;

```

and followed by several methods which access mandatory BAM fields, e.g.

```
/*! @function
 @abstract  Get the name of the query
 @param  b  pointer to an alignment
 @return    pointer to the name string, null terminated
 */
#define bam_get_qname(b) ((char*)(b)->data)

```
2. How does htslib access each of the optional fields, given this `structure`? Does it keep the order of the optional fields, or is this information unavailable?
3. When a user runs htslib or even `samtools view`, how does the BAM row "end", i.e. what data structure exists to show that the BAM row has no more optional fields, and that the algorithm should "move on" to the next BAM row.

Any clarification for how the source code works in this case is much appreciated.

# Answer
> 4 votes
> How does htslib access each of the optional fields, given this structure?

htslib uses a macro, `bam_get_aux`, for this purpose. You can check the source code of this macro — it essentially just advances the data pointer of the current record to the start of the optional data, which is after all other data. But you shouldn’t care about this, and just use the relevant exported functions for reading optional (“auxiliary”) data:

* `bam_aux_get` returns a pointer to an optional record, given by its two-character tag name.
* `bam_aux2‹C›` converts the pointer returned by the previous function to the actual type of the optional tag; here, `‹C›` is a type identifier: `i` for an integer, `Z` for a zero-terminated character string, etc. Refer to the documentation and the `htslib/sam.h` file for a complete list.

> Does it keep the order of the optional fields, or is this information unavailable?

As long as you don’t modify optional fields, the order is kept, but unspecified. *Do not rely on the order*!

> When a user runs htslib or even `samtools view`, how does the BAM row "end"

BAM records store a header that contains the record size.

# Answer
> 3 votes
1. In a BAM file, all of the auxiliary tags (the optional fields you refer to) are held what is effectively a `void *` field, one after the other. Each entry within this field is structured as a two `char`s storing the tag name (e.g., `NM`, followed by a `char` storing the type (e.g., `Z` or `i`), followed by the data in whatever format is appropriate (`char`, `float`, etc). Strings in the aux tags are `NULL` terminated, so their lengths aren't explicitly stored.
2. To access a given field, htslib has to iterate over all of them, in whatever order they happen to be stored, until it finds the correct one. The output order with `samtools view` is the order in which they're stored.
3. The size of `*data` is stored as `l_data`, so once it hits that then it knows that the record is finished.

Regarding point 1, this is a standard C trick of casting between data types. So, while `*data` is actually a `uint8_t *` (aka, a byte array), the first two bytes of each tag are `char`s and whatever is after them (well, after the byte specifying the type) is then cast to the appropriate type (be it `float` or `char*` or what not). This is all probably clearer if you look at the code.

---
Tags: sam, samtools, algorithms, htslib
---


# Question
Title: Cufflinks Error: sort order of reads in BAMs must be the same
I am running Cufflinks for transcriptome assembly using the .bam file generated by Hisat2. I tried both bam and sorted bam files

```
cufflinks --no-update-check -o ./ -p 15 accepted_hits.bam

```

it gives the following error:

```
Error: sort order of reads in BAMs must be the same

```

UPDATE: It is working fine now. It seems there were some issues with the reference file. The BAM files generated from Tophat2 also gave the same error. I downloaded the updated reference file and it worked great.

# Answer
> 6 votes
Just add `--dta-cufflinks` to your Hisat2 command, so that the output alignment file provides the attributes needed by Cufflinks (XS flags). This should do the job.

From manual:

> Report alignments tailored specifically for Cufflinks. In addition to what HISAT2 does with the above option (--dta), With this option, HISAT2 looks for novel splice sites with three signals (GT/AG, GC/AG, AT/AC), but all user-provided splice sites are used irrespective of their signals. HISAT2 produces an optional field, XS:A:\[+-\], for every spliced alignment.

---
Tags: rna-seq, assembly
---


title: Influence of pH on protein peptide mass

text:
# Question
Title: Influence of pH on protein peptide mass
I'm using the *Peptides* package in R to calculate the mass of peptides with the mw() method. The thing is that I want to calculate them for different pH values, and thus I wrote this little function to take the pH value into account according this image:

My code is this:

```
phmw = function(pept_seq , unimod_location, pH ){

  monoisotopic_hydrogen_mass = 1.00794

  if(pH < 7){

    # mol_weight = mol_weight_of_peptide + mol_weight_of_Hydrogen
    mol_weight = mw(pept_seq, monoisotopic = T) + monoisotopic_hydrogen_mass

  }else if(pH > 7){

    # mol_weight = mol_weight_of_peptide - mol_weight_of_Hydrogen
    mol_weight = mw(pept_seq, monoisotopic = T) - monoisotopic_hydrogen_mass

  }else if(pH == 7){

    # mol_weight = mol_weight_of_peptide
    mol_weight = mw(pept_seq, monoisotopic = T)
  }

  return(mol_weight)

}

```

Do you think that this approach is correct or not?

# Answer
> 6 votes
Your approach in general seems to be OK, but I have a few suggestions:

1. You are basing all your calculations on whether the pH is above or below 7. However the N and C termini have pKa values that determine whether these groups will be ionized (8 and 3, respectively). You should check whether the pH is above/below these values instead of 7.
2. Different amino acid side chains can also be ionized and should be accounted for when calculating the MW. Amino acid pKas can be found here
3. Code wise: you pass a variable unimod\_location in your function but it's never used. If it's not being used then take it out.

# Answer
> 3 votes
Amino-acid content (eg Lys, Asp) can have drastic effect on observance in a particular MS-mode.

You should consider including `Peptides::pI()` or `Peptides::charge()` to determine the likelihood of (de-)protonation in your calculations.

---
Tags: r, proteins
---


# Question
Title: How to convert the given mathematical computation (on biological problem) to mathematical fomula, equation?
I have crossposted this question in maths StackExchange. The problem is dominantly mathematical (this question) but the application of the problem is mainly biological. Hoping that people in this forum have faced similar problems, I am posting this in bioinformatics forum to get some ideas and probably solutions to it.

I am working on a problem that allows haplotype phasing. I have developed this computation method (expressed in detail below) and have developed a python code (not shown here) to solve the issue.

Now, I am having a hard time trying to translate this computation in mathematical and/or statistical language. I have revisited maths and I am able to comprehend the process at several places, but I cannot translate those ideas onto my computation (or algorithm).

---

Here is my computation (step by step) with a very workable example. I have tried to make it as comprehensive as possible, let me know if there is any confustion:

Fig 01: Example of the haplotype file

```
contig   pos     ref   all-alleles   ms01e_PI   ms01e_PG_al   ms02g_PI   ms02g_PG_al   ms03g_PI   ms03g_PG_al   ms04h_PI   ms04h_PG_al   ms05h_PI   ms05h_PG_al   ms06h_PI   ms06h_PG_al
2      15881764   .      .           4           C|T           6           C|T           7           T|T           7           T|T           7           C|T           7           C|T
2      15881767   .      .           4           C|C           6           T|C           7           C|C           7           C|C           7           T|C           7           C|C
2      15881989   .      .           4           C|C           6           A|C           7           C|C           7           C|C           7           A|T           7           A|C
2      15882091   .      .           4           G|T           6           G|T           7           T|A           7           A|A           7           A|T           7           A|C
2      15882451   .      .           4           C|T           4           T|C           7           T|T           7           T|T           7           C|T           7           C|A
2      15882454   .      .           4           C|T           4           T|C           7           T|T           7           T|T           7           C|T           7           C|T
2      15882493   .      .           4           C|T           4           T|C           7           T|T           7           T|T           7           C|T           7           C|T
2      15882505   .      .           4           A|T           4           T|A           7           T|T           7           T|T           7           A|C           7           A|T

```

This (fig 01) is my main data file. The main idea is to take a pool of several samples which have short phased haplotype blocks represented as `PI (phased index)` and `PG_al (phased genotype)`. For any sample, the site where two consecutive haplotype blocks are not joined represents the break point (and has different `PI` values). In the above example only sample `ms02g` has break point and needs to be phased. All other samples have haplotype block that bridges these two consecutive haplotype blocks and contains the data required to extend the phase state of `sample ms02g`.

Fig 02: Representing a break-point in the sample: ms02g

```
contig    pos    ms02g_PI    ms02g_PG_al
2     15881764     6         C|T
2     15881767     6         T|C
2     15881989     6         A|C
2     15882091     6         G|T
                   ×——————————×—————> Break Point
2     15882451     4         T|C
2     15882454     4         T|C
2     15882493     4         T|C
2     15882505     4         T|A

```

So, in the above haplotype file there is a breakpoint in `sample ms02g` at the position *15882091-15882451*. For sample ms02g(PI-6) the haplotypes are `C-T-A-G` and `T-C-C-T`. Similarly, at PI-4 the haplotypes are `T-T-T-T` and `C-C-C-A`. Since the haplotype is broken in two levels, we don’t know which phase from level-6 goes with which phase of level-4. But, all other samples have full haplotype intact that bridges this position. We can therefore use this information from other samples to join the two consecutive haplotype in sample ms02g.

Using markov-chain transition to extend the haplotype block:

Since, all other samples are completely phased bridging that breakpoint, I can run a markov-chain transition probabilities to solve the phase state in the `sample ms02g` . To the human eye/mind you can clearly see and say that left part of ms02g (PI-6, i.e C-T-A-G is more likely to go with right block of PI-4 C-C-C-A), thereby creating the extended haplotype block as `C-T-A-G-C-C-C-A` and `T-C-C-T-T-T-T-T`.

Below, I show how I can apply the first order markov transition matrix to compute the likelyhood estimates, calculate the log2Odds and then assign and extend the haplotype in proper configuration. And, I feed this logic to the computer using python.

Calculation of likelyhood estimates using markov transition:

### Step 01: prepare required haplotype configuration

* The top `PI-6 is Block 01` and the bottom `PI-4 is Block 02`.
* The phased haplotype in the left within each block is `Hap-A` and on the right is `Hap-B`.

Fig 03: representation of the haplotype breakpoint and block assignment

```
ms02g_PI    ms02g_PG_al
  6         C|T \
  6         T|C | Block 1
  6         A|C |
  6         G|T /
  ×——————————×—————> Break Point
  4         T|C \
  4         T|C |
  4         T|C | Block 2
  4         T|A /
           ↓   ↓
       Hap-A   Hap-B

```

So, the two consecutive blocks can be extended in one of the two possible haplotype configurations:

```
Parallel Configuration:
 Block01-HapA with Block02-HapA, so B01-HapB with B02-HapB
Vs. Alternate Configuration:
 Block01-HapA with Block02-HapB, so B01-HapB with B02-HapA

```

### Step 02: Compute transition matrix and estimate likelihood for each configuration.

Fig 04: Representation of the allele transition matrix (from alleles of former block-01 to alleles of later block-02).

```
Possible
transitions     ms02g_PG_al
 │              ┌┬┬┬  C|T \
 │              ││││  T|C | Block 1
 │              ││││  A|C |
 │              ││││  G|T /
 └────────────> ││││   ×—————> Break Point
                │││└> T|C \
                ││└─> T|C | Block 2
                │└──> T|C |
                └───> T|A /
                     ↓   ↓
                 Hap-A   Hap-B

**note: In this example I am showing transition only from first
  nucleotide of Block01 to first nucleotide of Block02. Actually we
  prepare transition from all nucleotides of B01 to all nucleotide of
  B02 in both parallel and alternate configuration.

```
* I count the number of transitions from each nucleotide of PI-6 to each nucleotide of PI-4 for each haplotype configuration across all the samples and convert them to transition probabilities.
* And multiply the transition probabilities from the first nucleotide in PI-6 to all nucleotides of PI-4. Then similarly multiply the transition probability from 2nd nucleotide of PI-6 to all nucleotides in PI-4, and so on.
* When transition probabilities are calculated for all possible combination (from each position of PI-6 to each position of PI-4), then I compute the cumulative transition probabilities for each possible haplotype configuration.

Fig 05 : Representation of nucleotide counts (emission counts) at positions 15881764 and 15882451.

```
pos\allele     A    T    G    C
15881764       0    8    0    4
15882451       1    7    0    4

```

Fig 06 : Representation of transition matrix counts (from pos 15881764) to (pos 15882451).

This transition matrix is computed from nucleotides (A,T,G,C) at block01 to nucleotides (A,T,G,C) at block02 for all the positions. Transition counts are then converted to transition probabilities.

```
from     to
          A    T    G    C
  A       0    0    0    0
  T       1    6.5  0    0.5
  G       0    0    0    0
  C       0    0.5  0    3.5

```

Note: if the PI matches between two blocks the transition are counted as 1, else 0.5. Sample ms02g itself is also taken as an observation.

### Step 03: Compute the maximul likelihood for each configuration

Fig 06 : Likelihood estimate for parallel configuration using transition counts (probabilities)

```
Parallel configuration:
Block-1-Hap-A (C-T-A-G) with Block-2-Hap-A (T-T-T-T)
  CtT × CtT × CtT × CtT = (0.5/4)*(0.5/4)*(0.5/4)*(0.5/4) = 0.000244
+ TtT × TtT × TtT × TtT = (0.5/2)*(0.5/2)*(0.5/2)*(0.5/2) = 0.003906
+ AtT × AtT × AtT × AtT = (0.5/3)*(0.5/3)*(0.5/3)*(0.5/3) = 0.0007716
+ GtT × GtT × GtT × GtT = (0.5/2)*(0.5/2)*(0.5/2)*(0.5/2) = 0.003906
——————— ————————— ——————— ————————— Max Sum (likelihoods) = 0.008828
                                    Average (likelihoods) = 0.002207

Block-1-Hap-B (T-C-C-T) with Block-2-Hap-B (C-C-C-A)
  TtC × TtC × TtC × TtA = (0.5/8)*(0.5/8)*(0.5/8)*(0.5/8) = 0.00001526
+ CtC × CtC × CtC × CtA = (2.5/10)*(2.5/10)*(2.5/10)*(2.5/10) = 0.003906
+ CtC × CtC × CtC × CtA = (1.5/8)*(1.5/8)*(1.5/8)*(1.5/8) = 0.001236
+ TtC × TtC × TtC × TtA = (0.5/4)*(0.5/4)*(0.5/4)*(0.5/4) = 0.000244
——————— ————————— ——————— ————————— Max Sum (likelihoods) = 0.0054016
                                    Average (likelihoods) = 0.0013504

```

note: \- "AtC" -\> represent "A" to "C" transition \- "+" represents the summation of the likelyhoods

Fig 07 : Likelihood estimate for alternate configuration using transition counts (probabilities)

```
Alternate configuration:
Block-1-Hap-A (C-T-A-G) with Block-2-Hap-B (C-C-C-A)
  CtC × CtC × CtC × CtA = (3.5/4)*(3.5/4)*(3.5/4)*(3.5/4) = 0.5861
+ TtC × TtC × TtC × TtA = (1.5/2)*(1.5/2)*(1.5/2)*(1.5/2) = 0.3164
+ AtC × AtC × AtC × AtA = (2.5/3)*(2.5/3)*(2.5/3)*(2.5/3) = 0.4823
+ GtC × GtC × GtC × GtA = (1.5/2)*(1.5/2)*(1.5/2)*(1.5/2) = 0.3164
——————— ————————— ——————— ————————— Max Sum (likelyhoods) = 1.7012
                                    Average (likelihoods) = 0.425311

Block-1-Hap-B (T-C-C-T) with Block-2-Hap-A (T-T-T-T)
  TtC × TtC × TtC × TtA = (6.5/8)*(7.5/8)*(7.5/8)*(6.5/8) = 0.5802
+ CtC × CtC × CtC × CtA = (6.5/10)*(7.5/10)*(7.5/10)*(6.5/10) = 0.237
+ CtC × CtC × CtC × CtA = (5.5/8)*(6.5/8)*(6.5/8)*(6.5/8) = 0.36875
+ TtC × TtC × TtC × TtA = (3.5/4)*(3.5/4)*(3.5/4)*(2.5/4) = 0.4187
——————— ————————— ——————— ————————— Max Sum (likelyhoods) = 1.60465
                                    Average (likelihoods) = 0.4011625

```

note:

```
  - the sum of the likelihoods > 1, in the above example.
  - So, we can rather use the product of the likelyhoods.

```

Fig 08 : Likelhood estimate of Parallel vs. Alternate configuration.

```
Likelihood of Parallel vs. Alternate configuration

= likelihood of Parallel config / likelihood of Alternate config

= (0.002207 + 0.0013504)/ (0.425311 + 0.4011625)
= 0.0043043  (i.e 1/232)
Therefore, haplotype is 232 times more likely to be phased
in "alternate configuration".

log2Odds = log2(0.0043043) = -7.860006
So, with our default cutoff threshold of "|5|", we will be extending
the haplotype in "alternate configuration"

```

### Final Output data:

```
contig   pos   ref   all-alleles   ms02g_PI   ms02g_PG_al
2   15881764   .    .                  6      C|T
2   15881767   .    .                  6      T|C
2   15881989   .    .                  6      A|C
2   15882091   .    .                  6      G|T
2   15882451   .    .                  6      C|T
2   15882454   .    .                  6      C|T
2   15882493   .    .                  6      C|T
2   15882505   .    .                  6      A|T

```

So, how do I translate my overall computation process in mathematical and statistical language.

I am looking to represent computation of:

* emission counts (and probabilities)
* transition counts (and probabilities)
* summation of the transition counts for each configuration
* computation of likelihood ratio and it's conversion to log2Odds

I am trying to represent my computation (algorithm) in the way it's shown in these links, but have a hard time doing it. I've read, this, and this, and this.

FYI - I am biologist but have managed to learn python and write a python program to work it out. The code isn't relevant in this question though.

Now, it's my time to work out some maths and represent my algorithm.

# Answer
I am not sure I understood all, but what you do is like doing a Markov chain but instead of saying that the previous position(s) decides the fate of the next one(s) is the block what determines which other block follows it.

Mathematically I think that could be explained in terms of linkage disequilibrium (and wikipedia page), but instead of using a single position you use a block. I am not familiarly enough with the field to see if a similar (or this method) has been previously described, but I think you could find inspiration about how to describe it.

A related field of work that also deals with this kind of data (and problems) is structural chromatin, where they work with TADs (their blocks of Topologically Associating Domains) and study which ones are near (spatially instead of probabilistic) with others regions, that field could be also interesting to look for formulas and how do they present this kind of analysis.
> 1 votes

---
Tags: python, bioconductor, hidden-markov-models, haplotypes
---


# Question
Title: How should the popular press compare similarity of genomes?
Note this is a question from a lay reader.

I've read in the popular press that ~1-4% of the genome of non-African *Homo sapiens* is inherited from Neanderthals, or that Melanesians derive a similar amount from the Denisovans.

This seems like a lot, until you read (elsewhere!) that Neanderthals are 99.7% identical to humans, so most of the fraction inherited is identical in any case.

Or taking it another way, reading that two individuals will typically have 99.9% identical DNA seems pretty similar, until you realise you are only 0.3% 'different' from a Neanderthal, or 15% 'different' from your dog.

It seems to me this use of percentages is confusing, misleading, and open to abuse. When comparing species (for example) I would think that they should speak of differences in terms of base pairs that we could not otherwise have had.

How should the lay (or popular science) press best communicate similarities or differences in genomes when talking about individuals, races, or species?

# Answer
1–4% is from an evolution point of view. 99.7% is from a sequence point of view. Because they are measuring different things, they can be compatible with each other. The correct interpretation is: 1–4% of a non-African genome is inherited from Neanderthals and the sequence of this 1–4% differs from modern human sequence by 0.3%. PS: so, for each non-African, Neanderthals contributed 1–4% * 0.3% of base pairs not present in modern human. This sounds tiny, but even a single SNP may lead to severe disease. This tiny fraction has known functional impact to diabetes and high-altitude adaption.

---

Others pointed out two related questions:

These questions are really about the ambiguity in the definition of "gene". Importantly, traditionally in popular science and population genetics, a gene could mean an allele. This makes differences.

For example, by common uses of "gene" in biology, the answer to "how many genes do we share with our mother" would be nearly 100% – every human has nearly the same set of genes. However, clearly this is not what the asker wants. By gene, he/she really means allele. The expected answer is 50% – we inherit 50% of alleles from our mother.

In the end, what is a "gene" exactly? I don't know. Biologists can often figure out the meaning of word "gene" in the right context, but non-biologists will sometimes have troubles with the ambiguity.
> 7 votes

# Answer
I don’t think it’s misleading. Your question seems to be based on the mistaken assumption that the percentages you mention — 1%–4% Neanderthal origin, 99.7% similarity between individuals, etc — should somehow add up to 100%.

As you note yourself, that’s obviously not the case. There was a similar question to the effect, “how can we share 99% of our genes with chimpanzees if we only share 50% of our genes with our mother and 50% with our father?” and the answer is fundamentally the same. Namely, that these numbers refer to different things. And the DNA or gene content we share with these individuals and groups of individuals isn’t additive, or mutually exclusive. Think of it as a Venn diagram:

(Ratios obviously exaggerated, since showing a diagram where mother and father share 99% wouldn’t be very readable.)
> 5 votes

# Answer
I don't think there's any way to summarize data like this in a way that is both comprehensive in the "covering all relevant facts and concepts" meaning and in the "able to be understood" meaning. Relatedness is not a simple quantity that can be summarized in one statistic. There's what percentage of total genome is shared. There's percentage of *coding* genome shared (I would expect the non-coding portion of the genome to be disproportionately represented in similarities between species). There's in-group versus out-group variation (i.e. how does the variation between Neanderthals and modern humans compare to the variation within modern humans and within Neanderthals). Then there's more abstract concepts such as consanguinity index and metrics of genetic drift (i.e. species that separated more recently have less drift between them).

Ideally, a journalist would pick the statistic that is most relevant to the issue being discussed in the article, give an overview of how the statistic is calculated and what it represents, and explain why they think it is the most relevant statistic, and perhaps explain what some other statistics are and whey the journalist chose not to focus on those. But this would be a lot to expect from a journalist, and for a reader to go through.

Fact is, 87.6% of percentages are largely or completely meaningless. If someone claims that teenagers get into accidents 20% more often than the general public, does that mean that the number of accidents per teenagers is 20% more than accidents per general person? Accidents per mile driven is 20% more? Accidents per hour? Is it comparing all teenagers to the general public? Teenagers with driver's licenses to the general public with driver's licenses? Is number of accidents, or dollar amount of damages? Then you have things like Simpson's Paradox and selection bias. A percentage is supposed to be one number divided by another, and giving just that number is implicitly asserting that it summarizes the data. But generally speaking, knowing those two original numbers, or at the very least where those numbers come from, is an important part of interpreting the data.
> 1 votes

---
Tags: genome
---


# Question
Title: counting individuals with the present gene and absent gene for each population
I have a folder containing 500 files from different individuals:

```
$ ls
HI.3577.005.Index_25.QPG4_KS2W-35_all_PA
HI.3577.005.Index_27.QPH4_SD1W-13_all_PA
HI.3577.005.Index_9.QPA2_KS2W-8_all_PA
HI.3577.007.Index_1.Q004B_SD1A-1_all_PA
HI.3577.007.Index_21.Q076C_SD2W-23_all_PA
HI.3577.007.Index_23.Q079_KS1A-15_all_PA

```

Each of these individuals belongs to a population (in the example above individuals belong to KS2W, SD1W, KS2W, SD1A, SD2W, and KS1A respectively). Each of my files looks as follows (this sample belongs to from "MK4" population):

```
$ head HI.3579.001.Index_9.Q089_MK4-4R_all_PA
Ha8_00040788 1
Ha4_00024045 1
Ha4_00025366 1
Ha16_00022130 0
Ha16_00023451 0
Ha8_00040789 1
Ha4_00025367 1
Ha4_00024046 0
Ha16_00022131 1
Ha16_00023452 1

```

or another sample from MK4 population

```
$ head HI.221.009.Index_9.P089_MK4-4R_all_PA
Ha8_00040788 1
Ha4_00024045 1
Ha4_00025366 1
Ha16_00022130 1
Ha16_00023451 1
Ha8_00040789 1
Ha4_00025367 1
Ha4_00024046 0
Ha16_00022131 1
Ha16_00023452 1

```

the first column is genes and the second columns refers to the presence of that gene (coded as 1) or absence of gene (coded as 0).

For all individuals belonging to each of the populations, I would like to count the total number individuals that have the gene and individuals that do not have the gene.

For example, for the 2 examples above, as they both belong to MK4 population my desired output is going to looke like this:

```
$ head MK4
Ha8_00040788   2  0
Ha4_00024045   2  0
Ha4_00025366   2  0
Ha16_00022130  1  1
Ha16_00023451  1  1
Ha8_00040789   2  0
Ha4_00025367   2  0
Ha4_00024046   0  2

```

The first column is the gene, the second column is the number of individuals in the MK4 population that have the gene and the third column is the number of individuals that do not have that gene.

I have made a file containing population IDs

```
$ head pop.txt
KS2W
SD1W
SD1A
KS1A

```

Can anyone suggest me how I could do this? thanks

# Answer
> 1 votes
I can show you how I would approach it in `R`. You'll have to adjust and customize it to your own data of course. I only had two files to play with, from your example. Especially for the `merge` step you should build a `for` loop around it so that it will merge every next file automatically, since you have 500 files.

Here an example for the KS2W files:

```
myFiles <- dir(pattern="KS2W")

myFiles
[1] "HI.3577.005.Index_25.QPG4_KS2W-35_all_PA.txt"
[2] "HI.3577.005.Index_9.QPA2_KS2W-8_all_PA.txt"

fileNames <- list()

for(i in 1:2){
    fileNames[[i]] <- read.table(myFiles[i], sep = " ")
}

merged <- merge(fileNames[[1]], fileNames[[2]], by.x = "V1", by.y = "V1", all = F, sort = F)

present <- rowSums(merged[,2:3])

zeros <- rowSums(merged[,2:3] == 0)

output <- cbind(present, zeros)

rownames(output) <- merged[,1]

output
              present zeros
Ha8_00040788        2     0
Ha4_00024045        2     0
Ha4_00025366        2     0
Ha16_00022130       1     1
Ha16_00023451       1     1
Ha8_00040789        2     0
Ha4_00025367        2     0
Ha4_00024046        0     2
Ha16_00022131       2     0
Ha16_00023452       2     0

```

To add more files with `merge`, you could add something like this:

```
merged <- merge(fileNames[[1]], fileNames[[2]], by.x = "V1", by.y = "V1", all = F, sort = F)

for(i in 3:30){
  merged <- merge(merged, fileNames[[i]], by.x = "V1", by.y = "V1", all = F, sort = F)
}

```

You'll have to adjust all the numbers and coordinates to fit your data.

# Answer
> 1 votes
The following python script takes a population name as argument, searches for all files for that population (assuming some regularity in the file names), and records the counts of "0" and "1" for all the genes encountered:

```
#!/usr/bin/env python3

from sys import argv
from glob import glob
from collections import defaultdict

filenames = glob("HI*Index_*_{pop}-*_all_PA".format(pop=argv[1]))

gene_counts = defaultdict(lambda : defaultdict(int))

for filename in filenames:
    with open(filename, "r") as genes_file:
        for line in genes_file:
            (gene, num) = line.strip().split()
            gene_counts[gene][num] += 1

for (gene, counts) in gene_counts.items():
    print(gene, counts["1"], counts["0"], sep="\t")

```

With the two example files you give being present in the current directory:

```
$ ./count_gene.py MK4
Ha8_00040788    2   0
Ha4_00024045    2   0
Ha4_00025366    2   0
Ha16_00022130   1   1
Ha16_00023451   1   1
Ha8_00040789    2   0
Ha4_00025367    2   0
Ha4_00024046    0   2
Ha16_00022131   2   0
Ha16_00023452   2   0

```

I haven't tested with python 2. If you use python 2, you should at least add `from __future__ import print_function` as first import, or adapt the `print` part.

Then, to run on all populations:

```
for pop in $(cat pop.txt)
do
    ./count_gene.py ${pop} > ${pop}
done

```

---
Tags: python, perl
---


# Question
Title: What is a poly(T) primer anchor sequence?
I am reading Tang et al. (2009).

In the method description it is stated that:

> A single cell is manually picked under a microscope and lysed. Then mRNAs are reverse-transcribed into cDNAs using a poly(T) primer **with anchor sequence** (UP1) and unused primers are digested.

What is an anchor sequence in this case?

1. Tang, F. et al. mRNA-Seq whole-transcriptome analysis of a single cell. Nat. Methods 6, 377–382 (2009).

# Answer
As shown in figure 1 of that paper, the anchor sequence is additional sequence that is added to the end of the polyT primer that provides a known and non-target binding site for additional PCR reactions. I gave a presentation about a month ago about the process that Oxford Nanopore Technology use in their strand-switching cDNA kit, which may be useful for understanding the process used by Tang *et al.*. Here are my strand-switch/PCR images from that presentation:

1. PolyT primer addition, with mismatched anchor sequence

2. Strand switch primer addition, with overhanging anchor sequence

3. PCR amplification, using primers that are complementary to the anchor sequences

In the case of the Tang *et al.* paper, they use a polyA-tailing method for strand-switching, but the process is essentially the same.
> 4 votes

---
Tags: scrnaseq, primer
---


# Question
Title: Question regarding the function pairwiseAlignment in R
I'm a bit confused and would appreciate your assistance, please.

I'm trying to use the `pairwiseAlignment` function from the `Biostrings` package to align two protein sequences.

I'm not quite sure where to put my reference sequence, the sequence I want to compare to. Is it under "subject" or "pattern"?

Also, is the score I get the edit distance? Is it the alignment score? Or is it the the alignment length? When I use local alignment, will the score be higher the more the sequences resemble each other? And when I use global alignment, will the score be lower the more the sequences resemble each other?

# Answer
For semantic purposes, it's better to use the reference sequence as the pattern. Note that only the pattern argument can be used as a multi-element vector, i.e. you can match multiple patterns against a single subject with this function, but not a single pattern against multiple subjects. Knowing that might alter what is defined as a "pattern" or a "subject".

The scores are associated with the alignment method that is used. Alignment scores are higher for more similar matches and lower for poor matches, but the weightings associated with matches and mismatches dictate what that actually means. It might be that a score of -2 happens to be the optimal alignment score, because all other possible alignments are more negative.

It's possible to use `pairwiseAlignment` to generate the levenshtein distance (most commonly-used edit distance) between two sequences by setting the match score to 0, the mismatch score to -1, gap opening penalty to 0, gap extension penalty to 1, and then negating the resultant `pairwiseAlignment` score.
> 2 votes

# Answer
According to the help page, it seems more logic to put the reference sequence in `subject`, and your query sequence in `pattern`, but it won't really matter for the outcome (except that the pattern sequence will be on top of the alignment).

For your second question:

> the scores for the optimal pairwise alignments is returned
> 2 votes

---
Tags: r, sequence-alignment
---


title: SNP data Download Help

text:
# Question
Title: SNP data Download Help
I am new to the field of Bioinformatics. Single Nucleotide Polymorphism or SNP data are available on the the Internet to download.SNP related data are valuable for variety of reasons ranging from disease related charatertcics and precision medicine to name the few.

Can somebody explain me the steps of downloading this kind of data from NCBI or other resource. For example if I want to download the the SNP in bunch of humans on a particular Gene from the specific start to end. What are the sequence of steps I would do online to download such data?

**Update**

I want SNP data for example for a particular gene say X from positions in range say base pairs 100,000 to 200,000 for individuals say 100 in txt or csv format or standard format that is available?

Like if I choose Gene X and position ranging from 100 to 109 base pairs in 5 persons I want following kind of data in txt or csv format with say two SNV at position 101 and 109

1. A **A** A A A G G G G **G**
2. A **G** A A A G G G G **A**
3. A **A** A A A G G G G **A**
4. A **A** A A A G G G G **A**
5. A **G** A A A G G G G **G**

regards

# Answer
> 2 votes
I would recommend to do this in R, using biomaRt.

For example if you want all the SNPs from IDH1 (ENSG00000138413) gene:

```
library(biomaRt)

ensembl<- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl",mart=ensembl)

IDH1 <- getBM(attributes=c('ensembl_gene_id','variation_name'), filters = 'ensembl_gene_id', values ="ENSG00000138413", mart = ensembl)

IDH1

```

You can add much more information, take a look at the help page of biomaRt.

---
Tags: snp, sequence-analysis, data-download
---


title: Simulate fast5 MinION run

text:
# Question
Title: Simulate fast5 MinION run
I'm planning on using Oxford Nanopore MinION sequencer to test a script that I'm working on. However, I would not like to run multiple MinION runs to validate the tool. I'd like to simulate various minION runs and get raw reads and all pore events in a fast5 file. Do you know of any Oxford Nanopore (MinION) sequencing simulator that allows to get fast5 as output? I looked into NanoSim but it seems like there is no option to get fast5 as an output.

# Answer
> 1 votes
Using a simulator for example MinION output doesn't make sense when there are an abundance of freely-available FASTQ and FAST5 sequences for general purpose use (and especially for the purposes of testing out algorithms). Simulation of sequences can only ever express the components of DNA that are explicitly modelled, which is never going to fully represent reality.

As long as you're doing your own private research, you can use the FAST5 sequences from Clive Brown's own genome:

https://github.com/nanoporetech/ONT-HG1

Bear in mind the license attached to that data:

https://github.com/nanoporetech/ONT-HG1/blob/master/LICENSE.md

> License to use for research use only (i.e. non-commercial) governed by the law of England and Wales.
>
> You agree to this License by downloading the data on this site whether in whole or in part.
>
> You cannot re-distribute modified copies of these data, or derivative works, without my express written permission. Where such permission is granted, you must cite the original work and its Copyrights and you must include this license with any work that includes any of my genome.

For a more open license, you can use the nanopore-wgs reads here (derived from human cell line NA12878):

https://github.com/nanopore-wgs-consortium/NA12878

> We encourage the reuse of this data in your own analysis and publications which is released under the Creative Commons CC-BY license. Therefore we would be grateful if you would cite the reference below if you do.

If you don't want human data, there's mouse parasite data here:

http://www.ebi.ac.uk/ena/data/view/PRJEB20824

And *E. coli* data here:

https://www.ebi.ac.uk/ena/data/view/PRJEB22772

... and if you just want 8h of sensor noise as a FAST5 file, I've got that here in one of my DropBox folders (as a ~70MB file):

https://www.dropbox.com/s/20j1anlz0motrs4/MINICOL233\_20170118\_FN\__MN15625\_sequencing\_throughput\_ONLL3127\_84599\_ch320\_read10302\_strand.fast5?dl=0

Reminder: if you do ever present the results of analyses based on a data source, it's a good idea to acknowledge the source (even if not expected). If nothing else, it will let the creator know that someone cares about their data, and perhaps encourage them to release more in the future.

# Answer
> 0 votes
I wrote a tool that will rerun an experiment. So all you would need to do is just download an ONT dataset from ENA or something and you can just rerun it as many times as you want.

https://github.com/mbhall88/taeper

---
Tags: minion
---


title: Converting a BAM file into VCF

text:
# Question
Title: Converting a BAM file into VCF
I have NGS illumina RNA-seq reads from *M. musculus* (mm10). I am trying to find variants along the strand portion of the reads in the refseq (mm10).

I mapped a pair of sequence files and generated a BAM file. Now I need to generate a VCF file from this BAM file. I’ve checked Galaxy for direct tools, but I couldn’t identify one. Any suggestions on how to turn a BAM file into VCF? Thank you in advance.

# Answer
It's not really possible to **convert** `bam` to `vcf`. `bam` is a mapping file, it does not contain the information about variants, this information needs to be inferred in process called variant calling. I find important to mention that it's not just a different format of the same thing.

How to **call variants** (a `vcf` file) from mapped reads (a `bam` file) is very broad question, it depends on the sequencing technology (Illumina, PacBio, ...), the type of sequencing (whole genome sequencing, RAD, RNA-seq, ...), the sampling (calling of individuals vs population variability), the quality of a reference (human vs denovo assembly made by a postdoc in your lab)... For a lot of cases the answer is simply GATK, partially because it has well written manual, but mainly just because many people use it.

Some more detailed cases are well answered in more specific questions :
> 19 votes

# Answer
You can use Freebayes, provided you have your BAM file and the reference genome.

Example:

```
freebayes -f genome.fa aln.bam > var.vcf

```

I'd suggest you go through the whole github documentation, and choose the option(s) more suitable for your needs (ploidy, read support, etc.)
> 5 votes

---
Tags: sam, vcf, variant-calling, variation
---


# Question
Title: tools to reconcile experimental transcripts with reference annotation
Looking for tools to reconcile alignment file of experimental transcripts mapped to genome (SAM/BAM) with the reference transcriptome annotation (GTF) from Ensembl (organism: *D. melanogaster*).

The aim would be to check which transcripts reported in the alignment file are also reported in the reference annotation, generating a new annotation (e.g., GTF) including the new transcripts as seen in the alignment file (experimental transcripts to genome).

Any idea?

-- Edit:

Up to now, I tried an in-house script that does the following:

1. generate a BED file from the BAM file with the exon coordinates (start/end) per transcript - using BEDOPS v2p4p27
2. convert the reference GTF to BED to retrieve the coordinates per annotated feature
3. compare the experimental vs the reference BED intervals
4. infer the suite of annotated exons per experimental transcript, thus, ideally, the isoform it corresponds to

However, with this simple approach, I only get around 20% of the experimental transcripts that can be identified in all their length using reference annotation. This might be due either to a too simple approach (e.g. if the start/end of the exon in the experimental transcript is shifted by 1 base as compared to the reference, it won't be detected), or to a real biological signal (since I am also using data from long-read sequencing technologies).

# Answer
I've never tried this myself, so I don't know how easy this is...

One option would be to start with GMAP, which is meant to align whole transcripts against the genome. The really nice thing about this is that it can directly produce GFF3 files. You can then use that with your Ensembl GTF with `cuffcompare` or whatever the equivalent is in stringTie. You should then be able to get the information you want from the transfrag codes.
> 6 votes

# Answer
As per my answer to @\_julien\_roux on twitter:

Trying to find novel transcripts within the context of an existing annotation is much less straightforward. You probably need to do a "genome-guided assembly" with Trinity and PASA: http://pasapipeline.github.io/#A\_ComprehensiveTranscriptome

We did something similar in much simpler organisms in our recent paper: http://dx.doi.org/10.1186/s12864-017-3505-0

Your biggest problems are going to be alternative spliced gene models, the large amount of fragmentation inherent in de novo transcriptome assembly and largely underreported non-coding RNAs in most genomes. It was easier for us as the dictyostelids have few and small introns - fly is going to be a different ball-game. Be prepared to spend a lot of time optimising parameters before finding a result you're happy with (maybe) ;)

You may think the fly genome pretty good and complete, but be prepared to find lots of errors, corrections and novelty. We were very surprised at the amount of new annotation we were able to find in *D. discoideum* \- an organism where every single gene had been manually curated!

Good luck!
> 5 votes

# Answer
Thanks to @Devon\_Ryan's suggestion, I had a go at using GMAP with our own *de-novo* assembly. Here's what I started with:

1. Genome assembled from nanopore long reads, corrected using Illumina cDNA reads
2. genome-guided Trinity transcriptome assembly, with reads grouped based on mapping to the assembled genome, then assembled by Trinity

The transcriptome assembly was filtered to only include transcripts that had reasonable coverage from cDNA reads, and further filtered to cluster similar transcripts and only use the longest isoform. This process is detailed in our paper. If I were to do that again, I'd probably create a supertranscript containing all exons, rather than the longest isoform.

I first created an index for the genome using `gmap_build`:

```
gmap_build -k 15 -j 8 -d 'MINBRA' \
  Nb_ONTCFED_65bpTrim_t1.contigs.hpcleaned.fasta -D MINBRA_GMAP

```

I then used gmap to map the transcriptome to the genome, generating gene models as a GFF3 file:

```
gmap -t 10 -d MINBRA -D ./MINBRA_GMAP \
  -A locusTagged_CDLI_CD98LMOHC50_TBNOCFED.tran.fasta \
  -f 2 > Nb_TrinTran_vs_Canu.gff3

```

Here's an example from the annotation file that was produced by GMAP:

```
...
tig00022147 MINBRA  gene    413327  414778  .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.path1;Name=MINBRA_TGG_483_c1_g1_i1
tig00022147 MINBRA  mRNA    413327  414778  .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.path1;coverage=100.0;identity=99.9;matches=1452;mismatches=0;indels=2;unknowns=0
tig00022147 MINBRA  exon    413327  414778  99  +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna1.exon1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna1;Target=MINBRA_TGG_483_c1_g1_i1 1454 1 .
tig00022147 MINBRA  CDS 413417  413557  100 +   0   ID=MINBRA_TGG_483_c1_g1_i1.mrna1.cds1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna1;Target=MINBRA_TGG_483_c1_g1_i1 231 91 .
###
tig00001947 MINBRA  gene    37526   38961   .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.path2;Name=MINBRA_TGG_483_c1_g1_i1
tig00001947 MINBRA  mRNA    37526   38961   .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna2;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.path2;coverage=100.0;identity=98.7;matches=1435;mismatches=1;indels=18;unknowns=0
tig00001947 MINBRA  exon    37526   38961   98  +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna2.exon1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna2;Target=MINBRA_TGG_483_c1_g1_i1 1454 1 .
tig00001947 MINBRA  CDS 37614   37750   97  +   0   ID=MINBRA_TGG_483_c1_g1_i1.mrna2.cds1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna2;Target=MINBRA_TGG_483_c1_g1_i1 231 91 .
###
tig00022147 MINBRA  gene    464550  465456  .   +   .   ID=MINBRA_TGG_489_c13_g1_i2.path1;Name=MINBRA_TGG_489_c13_g1_i2
tig00022147 MINBRA  mRNA    464550  465456  .   +   .   ID=MINBRA_TGG_489_c13_g1_i2.mrna1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.path1;coverage=100.0;identity=99.3;matches=907;mismatches=0;indels=6;unknowns=0
tig00022147 MINBRA  exon    464550  465456  99  +   .   ID=MINBRA_TGG_489_c13_g1_i2.mrna1.exon1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.mrna1;Target=MINBRA_TGG_489_c13_g1_i2 913 1 .
tig00022147 MINBRA  CDS 464909  465053  98  +   0   ID=MINBRA_TGG_489_c13_g1_i2.mrna1.cds1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.mrna1;Target=MINBRA_TGG_489_c13_g1_i2 507 361 .
###
...

```

GMAP told me on standard output what transcripts were not mappable to the genome (136 of 36830, or about 0.4%).

I've now used STAR to map these cDNA reads to the transcripts (using the GFF3 annotation file for initial evidence) in order to confirm that the mapping carried out by GMAP is reasonable. I ended up throwing out another 7.4% of transcripts. First, I generated a genome index using the generated gff3 file:

```
~/install/star/STAR-2.5.4b/bin/Linux_x86_64/STAR --runThreadN 10 \
   --runMode genomeGenerate --genomeDir ./star_genome \
   --genomeFastaFiles Nb_ONTCFED_65bpTrim_t1.contigs.hpcleaned.fasta \
   --sjdbGTFfile Nb_TrinTran_vs_Canu.gff3 \
   --sjdbGTFtagExonParentTranscript Parent

```

Followed by mapping my RNASeq experiments using STAR. STAR can sort the genome mapping, but not the transcriptome mapping, so an additional sort/index was done on that:

```
for x in $(seq 1 9); do
  ~/install/star/STAR-2.5.4b/bin/Linux_x86_64/STAR --runThreadN 10 \
     --readFilesCommand zcat --genomeDir ./star_genome \
     --readFilesIn *-1563-0${x}-*.fastq.gz \
     --outFileNamePrefix mapped/1563-0${x}_vs_BNOCFED \
     --outSAMtype BAM SortedByCoordinate --quantMode TranscriptomeSAM
  samtools sort mapped/1563-0${x}.toTranscriptome.*.bam \
    > mapped/1563-0${x}.toTranscriptome.sortedByCoord.out.bam
  samtools index mapped/1563-0${x}.toTranscriptome.sortedByCoord.out.bam
done

```

BAM hits were collated into a single file:

```
(for x in $(seq 1 9);
    do samtools idxstats 1563-0${x}.toTranscriptome.sortedByCoord.out.bam | \
      perl -pe "s/^/1563-0${x}\t/";
 done) | gzip  > readCounts_all_vs_transcriptome.tsv.gz

```

Which was then processed using R to identify unmapped transcript annotations:

```
#!/usr/bin/Rscript
data.df <-
    read.delim("readCounts_all_vs_transcriptome.tsv.gz", header=FALSE,
               col.names=c("Sample","Transcript",
                           "Length","Mapped","Unmapped"));
gene.counts <- tapply(data.df$Mapped + data.df$Unmapped,
                      list(data.df$Transcript, data.df$Sample), sum);
sum(rowSums(gene.counts)==0) / nrow(gene.counts);
TranscriptsToRemove <- rownames(gene.counts)[rowSums(gene.counts) == 0];
cat(TranscriptsToRemove, file="transcripts_to_remove.txt", sep="\n");

```

The unmapped annotations were then excluded from the gff3 file:

```
grep -v -F -f mapped/transcripts_to_remove.txt Nb_TrinTran_vs_Canu.gff3 \
  > mappable_Nb_TrinTran_vs_Canu.gff3

```
> 1 votes

---
Tags: sequence-alignment, transcriptome, software-recommendation
---


# Question
Title: Subset data frame by similar elements from rows and column between two data frames
My first data frame whose first column i intend to use which is PATIENT\_ID

```
PATIENT_ID  SEX RACE    FAB AGE BM_BLAST_PERCENTAGE WBC
TCGA-AB-2805    Male    WHITE   M0  77  67  92
TCGA-AB-2822    Male    WHITE   M0  65  99  2.9
TCGA-AB-2831    Male    WHITE   M0  59  61  1.4
TCGA-AB-2890    Male    NH/C    M0  69  60  6.1
TCGA-AB-2907    Male    WHITE   M0  69  82  1.6
TCGA-AB-2913    Male    WHITE   M0  61  95  61.6
TCGA-AB-2922    Male    WHITE   M0  83  75  77.9
TCGA-AB-2936    Female  NH/C    M0  60  50  2.1

```

My second data frame where i have the respective patient ID as individual column with expression values

```
Hugo_Symbol TCGA-AB-2803    TCGA-AB-2805    TCGA-AB-2806    TCGA-AB-2807    TCGA-AB-2808
A1BG    0   0   0   0   0
A1CF    0   0   0   0   0
A2BP1   71.4095 48.5475 50.8339 45.7498 47.1328
A2LD1   2.0576  3.3467  5.7288  0.993   5.4988
A2ML1   467.0782    16.0643 7.0019  343.5948    255.9152
A2M 11.3169 0.6693  17.823  116.1867    182.2467
A4GALT  0   2.6774  1.2731  4.9652  0.7855
A4GNT   0   0   0   0   0
AAA1    514.4033    532.7979    504.1375    697.1202    443.8335
AAAS    0   0   0   0   0
AACSL   190.3292    74.2972 84.6595 559.0864    317.3606
AACS    0   0   0   1.9861  0
AADACL2 0   0   0   0   0
AADACL3 0   0   0   0   0
AADACL4 1.0288  0   0.6365  0   0

```

What I tried so far transposing the second data frame so that i can get the individual column as my first column ,and subset by using dplyr library but so far im not able to ..any help or suggestion how to do it would be highly appreciable

I would like to have my expected data frame as such

```
PATIENT_ID  Gene B  GeneA   GeneC   GeneA   GeneA
TCGA-AB-2803-03 0   0   0   92.5926 876.5432
TCGA-AB-2805-03 0.6693  1.3387  0   92.6305 421.6867
TCGA-AB-2806-03 0   14.1884 16.3654 119.599 625.0796
TCGA-AB-2807-03 2.006   9.3247  13.5154 178.4211    537.2393
TCGA-AB-2808-03 0   42.3016 28.3975 170.2042    653.5742
TCGA-AB-2810-03 0   1.5668  8.9665  152.0747    1087.5648
TCGA-AB-2811-03 0   0   6.9556  30.3103 1414.6602
TCGA-AB-2812-03 0   14.3289 11.9466 120.3591    962.5575
TCGA-AB-2813-03 0   3.4938  1.3813  51.1883 594.7593
TCGA-AB-2814-03 0   56.8659 118.7051    87.1817 472.6912
TCGA-AB-2815-03 0   2.0692  3.6986  104.6575    502.5234

```

# Answer
If you're still looking for an answer on this, I would use the `gather` and `spread` functions in the `tidyr` package (part of the `tidyverse` suite of packages - which also includes `dplyr` and `ggplot2`:

```
library(tidyverse)
df2_long <- gather(PATIENT_ID, value, -Hugo_Symbol)
df1 %>%
  left_join(df2_long) %>%   # merges the 2 data frames together
  spread(Hugo_Symbol, value)

```
> 1 votes

# Answer
You can take the PATIENT\_ID from the first dataframe and use them to select from your second dataframe. Than transpose, so something like this:

```
ID_Names <- First_df$PATIENT_ID

rownames(Second_df) <- Second_df$Hugo_Symbol

Selected_df <- Second_df[, ID_Names]

Expected_df <- t(Selected_df)

```
> 4 votes

---
Tags: r
---


# Question
Title: Is it OK to use Blast+ to query NCBI's 16s rRNA database with 16s DNA sequences?
I have 16S DNA sequences from diversity sequencing, and want to query these against NCBI's 16S rRNA database. If I ultimately want to get taxids for the species that my sequences have, does it make a difference that my sequences are 16S DNA and not RNA? I am using Blast+.

I imagine this is how analysis of 16S sequences normally goes (get 16S DNA sequences, query them against a 16S DB), but does anyone know if Blast+ or the 16S rRNA DB itself allow for DNA queries?

In other words, if I were to blast a sequence of 16S DNA and its exact transcript RNA, would I get the same results?

# Answer
> 2 votes
Your confusion is understandable, since it isn't explained in the NCBI description of the database, but most of the sequences in the database are derived from DNA that encodes the 16S rRNA *gene*, not from the rRNA itself. In some cases, sequences are derived from cDNA that is synthesized using the actual rRNA as a template, but it's still DNA sequencing. Direct RNA sequencing is extremely uncommon. Depending on the cDNA sequencing protocol (first strand vs. second strand sequencing), the transcript sequence might be the reverse-complement of the 16S rDNA, but the only difference that you would see is that the numbering of the bases would be reversed. If you are sequencing rDNA from organelles, there may also be a few bases that are different, due to RNA editing.

---
Tags: ngs, blast
---