# Question
Title: How to indicate the END of a haplotype block in VCF?
In VCF I know how to indicate that two genotypes are in the same phase by using consecutive "0|1" and "1|0" genotype fields, for example. However, how do I deal with the case that the first two genotypes are in phase and the second 2 are in phase, but there is no phase between the 2nd and 3rd?

According to the format I think I need to write:

```
chr1 100 . A T . . . GT 1|0
chr1 200 . A T . . . GT 0|1
chr1 300 . A T . . . GT 1|0
chr1 400 . A T . . . GT 0|1

```

but that implies that the SNP at 200 and the SNP at 300 are in phase, when in fact they are not.

How to indicate the END of a haplotype block in VCF?

# Answer
According to the VCF specification, haplotype blocks are defined by the "Phase Set" \[PS\] tag:

> PS : phase set. A phase set is defined as a set of phased genotypes to which this genotype belongs. Phased genotypes for an individual that are on the same chromosome and have the same PS value are in the same phased set. A phase set specifies multi-marker haplotypes for the phased genotypes in the set. All phased genotypes that do not contain a PS subfield are assumed to belong to the same phased set. If the genotype in the GT field is unphased, the corresponding PS field is ignored. The recommended convention is to use the position of the first variant in the set as the PS identifier (although this is not required). (Non-negative 32-bit Integer)

A good example of this can be found on a 10X Genomics support page about phasing:

```
CHROM   Pos REF ALT GT  PS
chr1    1000    A   C   0|1 1000
chr1    1010    T   G   1|0 1000
chr1    2000    C   T   0|1 2000
chr1    2005    T   G   0/1 2000
chr1    2008    G   C   0|1 2000

```

> In this example we have two phase blocks, denoted by PS=1000 and PS=2000. PS=1000 spans position 1000-1010, and PS=2000 spans position 2000-2008. In PS=1000, haplotype 1 contains the REF A allele at position 1000, and the ALT G allele as position 1010, while haplotype 2 contains the ALT C allele at position 1000 and the REF T allele at position 1010.

Note that this phase set identifier needs to be a single numeric value, which means that sub-blocks and overlapping blocks cannot be defined by this field.
> 3 votes

# Answer
Two options that immediately come to mind, depending on what you hope to accomplish down the pipeline: you might simply mark the unphased sites with a "/" instead of "|" to indicate they are not phased; or, if you have multiple haplotype blocks that are all independent of each other, you could use the chromosome field to identify each unique haplotype (e.g. chr1h1, chr1h2, chr2h1, chr2h2, etc).
> 1 votes

# Answer
There is a new idea of reporting the "phased haplotype/genotype block" with unique block ID.

phaser uses "PI" (unique phase-index) and "PG" (phased-genotype).

GATK is also using this idea of "unique-block" and "phased-genotype"; See the "HP" tag.

Personally, I feel the representation of "phaser" is way better and comprehensive.

I would represent your data as:

```
chr1 100 . A T . . . GT:PI:PG    1/0:3:1|0
chr1 200 . A T . . . GT:PI:PG    0/1:3:0|1
chr1 300 . A T . . . GT:PI:PG    1/0:4:1|0
chr1 400 . A T . . . GT:PI:PG    0/1:4:1|0

```

where first two lines have same "PI" (i.e 3) stating they are phased with each other. Same with the later two lines which have another "PI" (i.e 4).
> 1 votes

---
Tags: file-formats, vcf, haplotypes
---


# Question
Title: How to calculate statistical significance of sequence motifs
I would like to show that certain types of RBP motifs are enriched in RNA editing islands (i.e. clusters of RNA editing). However, I am unsure about how to think about sequence motifs with respect to their occurrence in other genomic features.

I understand how to find the probability of a motif at a location. i.e.

P(sequence is at position i) = P(A)^\[A\] * P(C)^\[C\] * P(G)^\[G\] * P(T)^\[T\]

where (.) is the base pair and \[.\] is the number of bps.

Instead what I would like to find is:

P(sequence S is contained in a feature type T) = ???

where feature type T is an gene, intron, editing island, etc. I think I should incorporate length since I will mainly be comparing genes or introns vs. editing islands. Also, I am not sure what to do about the editing islands being located within genes. How can I keep from counting the same motif twice?

Any ideas would be greatly appreciated. Thank you for your time.

# Answer
> 4 votes
As far as I remember the exact probability computation is an open problem. The reason is that potential motifs can overlap, which makes probability computations for an arbitrary string non-trivial, and depends on the motif.

For example if he have a binary string of four digits, the probability of "01" will be 11/16, while the probability of "11" will be 8/16=1/2.

The simplest approximation is to assume that probability of the motif on every position of the sequence is equal and independent. In such case (this comes from this course):

(Probability of text of length N having at least t occurrences of a k-mer Pattern, A is the number of letters in the alphabet, n = N – t * k)

This approximation is very rough, there are better here.

# Answer
> 1 votes
As Iakov Davydov points out above mathematical calculation is very difficult. Particularly if you do not have uniform nucleotide, dinucleotide etc. frequencies. Thus I'd recommend using a method based on randomisation.

You'd want to randomise only within genes since the likelihood of getting an editing island is likely to be dependent on expression level of the RNA.

You might like to look at what people analysing CLIP data do for these sorts of analyses. See for example Wang et al. Or I have code for this sort of thing in my `iCLIPlib`. See the module `kmer.py`. (WARNING: this code is still in development, and is not "released" yet, so I offer no promises about bugs or how easy it would be for your to get working). However, these approaches generally only *rank* the significance of kmers rather than giving absolute significance values.

Genome Annotation Tester (GAT) might be of use here. You would define genes as your workspace, islands as your annotations, and motif matches as your segments of interest.

---
Tags: statistics, motifs
---


# Question
Title: How to QC overamplified shotgun library?
I made some NEB ultra II whole-genome shotgun libraries with around ~200ng of template going into the indexing PCR step. The PCR was run for 13 cycles and I ended up having high library concentrations between 60-90 ng/µL. Should I be worried about complexity or PCR duplicates? How would all of you QC to make sure the library is complex enough and not overamplified too badly?

# Answer
FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) is a good start: it will show you observed/expected ratios for the most overrepresented sequences, and running it on a few datasets you know to be high quality should give you intuition for what looks good and what doesn't.

My labmate also produced something called preseq (http://smithlabresearch.org/software/preseq/) that uses the library complexity of an initial sequencing experiment to predict how many more unique reads you would get with additional sequencing of the same library. It could give you some intuition for whether or not it is worth sequencing your library deeper to get the information you need.
> 1 votes

---
Tags: ngs, illumina, genome-sequencing
---


# Question
Title: How to interpret the SMRT Link base modification algorithm output?
I have got PacBio data from a Sequel machine at a coverage \> 50X. I have run it in SMRT Link with their Base Modification Detection algorithm. I am now trying to figure out how to interpret the output and where to find the information about the modified bases. One of the outputs is a GFF3 file. The top of this file looks like this:

```
> ##gff-version 3

```

...

```
> I kinModCall  modified_base   1168    1168    31  -   .   coverage=48;context=TTAAACTAATTTCCAAAAATTTTGAACAGTTTTTAATACAA;IPDRatio=2.96
> I kinModCall  modified_base   1497    1497    32  -   .   coverage=50;context=CTTGAAAACACGAATATCTGTTTAGACAATAGTTTTGCCAA;IPDRatio=3.19
> I kinModCall  modified_base   1509    1509    33  -   .   coverage=48;context=CATGAAAAATTTCTTGAAAACACGAATATCTGTTTAGACAA;IPDRatio=2.84
> I kinModCall  modified_base   1529    1529    39  -   .   coverage=46;context=AAGAAAAGACCAAGAAAAACCATGAAAAATTTCTTGAAAAC;IPDRatio=3.84
> I kinModCall  modified_base   1904    1904    36  -   .   coverage=54;context=TAAAACCTGACAGAAAGAATTGGTCCCGTCATCTACTTTCT;IPDRatio=2.86
> I kinModCall  modified_base   2152    2152    31  -   .   coverage=59;context=TTAAAAATTAATATAAAATCAGGGAAATTTTTTTAAATTTT;IPDRatio=2.98
> I kinModCall  modified_base   2226    2226    36  -   .   coverage=63;context=TGGGCGGGGCGTCTCTGCCGGGAAAATTTTGCGTTTAAAAA;IPDRatio=2.59
> I kinModCall  modified_base   2228    2228    34  -   .   coverage=52;context=CGTGGGCGGGGCGTCTCTGCCGGGAAAATTTTGCGTTTAAA;IPDRatio=2.98
> I kinModCall  modified_base   2237    2237    33  -   .   coverage=60;context=GCACGGTTTCGTGGGCGGGGCGTCTCTGCCGGGAAAATTTT;IPDRatio=3.07
> I kinModCall  modified_base   2240    2240    44  -   .   coverage=63;context=GCGGCACGGTTTCGTGGGCGGGGCGTCTCTGCCGGGAAAAT;IPDRatio=2.89
> I kinModCall  modified_base   2311    2311    31  -   .   coverage=63;context=AAAATCGACGAAAATCGGTATAAAATCAAATAAAAATAGAA;IPDRatio=4.35
> I kinModCall  modified_base   2313    2313    39  -   .   coverage=62;context=GAAAAATCGACGAAAATCGGTATAAAATCAAATAAAAATAG;IPDRatio=4.00
> I kinModCall  modified_base   2314    2314    34  -   .   coverage=63;context=AGAAAAATCGACGAAAATCGGTATAAAATCAAATAAAAATA;IPDRatio=2.63
> I kinModCall  modified_base   2333    2333    35  -   .   coverage=62;context=AAAAAAAAGAGAAAAAATGAGAAAAATCGACGAAAATCGGT;IPDRatio=3.77
> I kinModCall  modified_base   2334    2334    35  -   .   coverage=63;context=CAAAAAAAAGAGAAAAAATGAGAAAAATCGACGAAAATCGG;IPDRatio=2.77

```

I expected something along the lines of 'chromosome I position 1497 methylcytosine', etc.

So it seems that the algorithm has detected modified bases at the following positions but it cannot call them, maybe because of a lack of coverage?

Update: the screenshot below shows there is no CSV file in the Data tab of the completed analysis.

# Answer
It sounds like you're looking for the `modifications.csv` file, rather than the `modifications.gff` file. The `.gff` file is for viewing in `SMRT View` (or other viewers). You can see schema for both files at PacBio's Methylome Analysis Technical Note. Specifically, `modified_base` means the algorithm could not determine the exact base and modification.

I also found a video tutorial for performing the analysis in SMRTLink on YouTube

**EDIT**: Per PacBio's Methylome Analysis Technical Note documentation, the `.csv` "...is also produced when motif analysis is not active, such as with the RS\_Modification\_Detection analysis protocol." Not sure why it's missing, perhaps because the algorithm couldn't clearly identify any of the modifications. You'll probably need to contact tech support on this issue.

In either case, the `.gff` you posted includes the significant hits (p \< 0.01), but the algorithm was unable to classify them, likely because of insufficient coverage. PacBio's Methylome Analysis Technical Note documentation gives specific coverage recommendations for each event type, but 100x is the lower bound. Lots of good information there.

Probably not what you wanted to hear, but I hope that gets you going.
> 2 votes

---
Tags: long-reads, pacbio
---


# Question
Title: Creating frequency plot from weblogo with frequency unit on the Y axis
Hello I am trying to create frequency plot for a motif with frequency unit on the Y axis from 0 -100 nt as the unit as represented in this paper.(Seiler et al., 2018)

# Answer
Within weblogo, add the `--units probability` argument and edit the figure with Illustrator (PDF/SVG) or Inkscape (SVG) to modify the y-axis, changing fractional probability (`[0,1]`) to a percentage (`[0,100]`) by multiplying by 100. Refer to the manual for more detail.
> 4 votes

---
Tags: sequence-alignment, motifs
---


# Question
Title: How can I reproduce a manual NCBI search with Biopython Entrez module?
I'm trying to make a Biopython script that reproduces a manual search on the NCBI website.

My manual search gives me the following URL:

https://www.ncbi.nlm.nih.gov/nuccore?term=ANOS1\[gene%20name\]%20AND%20refseq\[filter\]

And the "Search details" box indicates:`ANOS1[gene name] AND refseq[filter]`.

I get 338 results that I "Send to" "File" as "GI List", save as `ANOS1_orthologs_manual.txt`.

**How can I do the same using Biopython?**

# Answer
This can be done by using the "Search details" as a search term in `Entrez.esearch`:

```
import os
OPJ = os.path.join
base_dir = os.getcwd()

from Bio import Entrez
Entrez.email = "my_email_here"
query = "ANOS1[gene name] AND refseq[filter]"
handle = Entrez.esearch(term=query, db="nuccore", retmax=400)
ids = Entrez.read(handle)["IdList"]
with open(OPJ(base_dir, "ANOS1_orthologs_max400_ids.txt"), "w") as ids_file:
    ids_file.write("\n".join(ids))
    ids_file.write("\n")

```
> 6 votes

---
Tags: python, public-databases
---


# Question
Title: How to retrieve the yeast locus tag from uniprot via sparql?
Let's say I want to retrieve the yeast-specific identifier for a certain protein, in the example below it will be P00330 which I would like to link to YOL086C.

When I go to uniprot's sparql UI and type

```
PREFIX up:<http://purl.uniprot.org/core/>
PREFIX rdfs:<http://www.w3.org/2000/01/rdf-schema#>

SELECT *
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein rdfs:seeAlso ?db .
    ?db up:database <http://purl.uniprot.org/database/SGD> .

}

```

it returns

```
protein db
http://purl.uniprot.org/uniprot/P00330  http://purl.uniprot.org/sgd/S000005446

```

How would I now go on and find the `YOL086C`? When I add the line

```
?db ?p ?o

```

to the query above, I don't see anything that would directly lead me to the identifier.

I also tried to access another database

```
?db2 up:database <http://purl.uniprot.org/database/EuPathDB>

```

but this returns all kind of databases while I am interested only in the `FungiDB`. How could I select this one?

How could I link the two IDs to each other?

# Answer
> 4 votes
Ok, figured it out by looking at the scheme uniprot provides:

```
PREFIX up:<http://purl.uniprot.org/core/>

SELECT ?protein ?yeast_id
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein up:encodedBy ?gene .
    ?gene up:locusName ?yeast_id
}

```

So one starts at the accession number in the middle and then uses the path to the actual `locusName`.

Also the path via the database works fine, I just did not use the `seeAlso` correctly (`protein` was only bound to `db` but not to `db2`, that's why it returned all entries):

```
PREFIX up:<http://purl.uniprot.org/core/>
PREFIX rdfs:<http://www.w3.org/2000/01/rdf-schema#>

SELECT ?protein ?db ?db2
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein rdfs:seeAlso ?db .
    ?protein rdfs:seeAlso ?db2 .
    ?db up:database <http://purl.uniprot.org/database/SGD> .
    ?db2 up:database <http://purl.uniprot.org/database/EuPathDB>

}

```

This will return the desired output (one would still have to process the `db2` entry, of course):

```
protein db  db2
http://purl.uniprot.org/uniprot/P00330  http://purl.uniprot.org/sgd/S000005446  http://purl.uniprot.org/eupathdb/FungiDB:YOL086C

```

---
Tags: uniprot
---


# Question
Title: seqret error when trying to transform GFF3 and FASTA files into EMBL format
I am using the tool seqret from emboss to transform an annotation file in GFF3 format and a FASTA file into an EMBL file because Wormbase does not supply an EMBL file with annotation and sequence together. I am running it on Mac OS X. I got the gff3 and fasta files from Wormbase. I installed emboss using conda from the bioconda channel. I use the following command:

`seqret -sequence c_elegans.PRJNA13758.WS263.genomic.fa -feature -fformat gff -fopenfile c_elegans.PRJNA13758.WS263.annotations.gff3 -osformat embl -auto -outseq c_elegans.PRJNA13758.WS263.embl`

It fails with the following error:

```
/Users/me/anaconda/bin/seqret`: line 7: 54774 Killed: 9               _seqret "$@"

```

`/Users/me/anaconda/bin/seqret` is a short file that indeed contains `_seqret "$@"` on line 7:

```
#!/bin/sh
BIN_DIR=$(dirname $(which seqret))
export EMBOSS_ACDROOT=$BIN_DIR/../share/EMBOSS/acd/
export EMBOSS_DATA=$BIN_DIR/../share/EMBOSS/data/
export PLPLOT_LIB=$BIN_DIR/../share/EMBOSS/
_seqret "$@"

```

I tried to erase it but that did not help.

I also tried EMBLmyGFF3 in a Python 2 environment:

`EMBLmyGFF3 c_elegans.PRJNA13758.WS263.annotations.gff3 c_elegans.PRJNA13758.WS263.genomic.fa --topology linear --molecule_type 'genomic DNA' --transl_table 1 --species 'Caenorhabditis elegans' --locus_tag CELE --project_id PRJNA13758 -o c_elegans.PRJNA13758.WS263.annotations.embl`

but I just had this message on display for hours:

```
#############################################################################
# NBIS 2016 - Sweden                                                        #
# Authors: Martin Norling, Niclas Jareborg, Jacques Dainat                  #
# Please visit https://github.com/NBISweden/EMBLmyGFF3 for more information #
#############################################################################

```

Update 11th March:

I started EMBLmyGFF3 yesterday at 5 pm and this morning it started outputting messages such as:

```
08:22:08 ERROR feature: >>inverted_repeat<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 ERROR feature: >>assembly_component<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 WARNING feature: Unknown qualifier 'genbank' - skipped
08:22:08 ERROR feature: >>tandem_duplication<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 WARNING feature: Unknown qualifier 'public_name' - skipped
08:22:08 WARNING feature: Unknown qualifier 'polymorphism' - skipped
08:22:08 WARNING feature: Unknown qualifier 'other_name' - skipped
08:22:08 WARNING feature: Unknown qualifier 'consequence' - skipped
08:22:08 WARNING feature: Unknown qualifier 'variation' - skipped
08:22:08 ERROR feature: >>biological_region<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 WARNING feature: Unknown qualifier 'balancer_type' - skipped
08:22:08 WARNING feature: Unknown qualifier 'balancer' - skipped
08:22:08 ERROR feature: >>protein_match<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 ERROR feature: >>transcript_region<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 ERROR feature: >>transcribed_fragment<< is not a valid EMBL feature type. You can ignore this message if you don't need it.

```

...

```
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:45:36 ERROR feature: >>piRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:46:00 ERROR feature: >>pseudogenic_tRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:47:59 ERROR feature: >>antisense_RNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:47:59 ERROR feature: >>lincRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:48:03 ERROR feature: >>pre_miRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:48:03 ERROR feature: >>miRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:50:37 ERROR feature: >>snRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
09:04:00 ERROR feature: >>miRNA_primary_transcript<< is not a valid EMBL feature type. You can ignore this message if you don't need it.

```

For some of these features, it seems OK to not report them. But for others, such as piRNA, I need to find a solution. *C. elegans* has got 1000s of piRNAs, it would be nice to have them in the EMBL file. I need to find a list of valid EMBL feature types.

Update + ~ 24h computation:

...

```
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
12:21:19 ERROR feature: >>scRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
Conversion done

```

My EMBL file is only 1.19 GB compared to the GFF3 3.34 GB, so I probably lost a lot of annotation that is not supported by EMBL.

# Answer
That error is telling you that the script you were running was killed by signal `9`, also known as `SIGKILL`, while executing the command `_seqret "$@"`. This is the way the operating system will kill processes that are asking for more memory than it can provide.

Note that removing `_secret "$@"` is not the solution since that's what's actually doing the job. `$@` is a bash array that contains all the arguments passed to a bash script. The `/Users/me/anaconda/bin/seqret` file is just a simple script that sets up a few variables and then launches `_seqret`, passing it the arguments you gave to `/Users/me/anaconda/bin/seqret` as `$@`. So deleting that lines just makes the script useless.

Based on your question, I am guessing you're running macOS on a laptop, so it isn't unreasonable that something that needs to read an entire genome and what looks like a full set of annotations will run out of memory on a laptop. Of course, this is a pretty small genome, but the annotations are still 396MB compressed, so depending on your machine's specs and what kind of internal processing the command is doing, it isn't too surprising that it fails.

So, either run this on a more powerful machine, or try splitting the annotations file into smaller files and running it sequentially on each of those.
> 2 votes

# Answer
I'm one of the EMBLmyGFF3 developers. If the compressed annotation is 396MB as mentioned by @terdon it is not a small annotation... In the worse case we have tested we never waited more than few minutes. We haven't done any particular speed test. It would be valuable if you could tell us how long it has taken in your case (If you didn't kill the process before the end...).

At least we could try to implement some nice features, like to see that the tool is working , or/and see which percentage of the file has been already processed.

If it is really slow we could think to try to improve its speed.
> 3 votes

---
Tags: file-formats, format-conversion, sequence-annotation, emboss
---


title: Impute phenotype under some constraints

text:
# Question
Title: Impute phenotype under some constraints
Errors happen frequently in the lab and I got one sample that we mixed one important information (the region of the sample) between two samples. Now we don't know which sample is from which region, but we know all the other variables (Date of the sample, patient, and other variables).

So far, my solution has been to make that values as `NA` to be able to analyze it, but I would like to know if there is any solution. If I can convey to test or other analysis downstream that sample A and sample B are from Location1 or Location2 and if A is from Location1, B is from Location2 and otherwise round.

There are several packages for imputing quantitative data, but I couldn't find any to store this kind of information or to impute qualitative data. I can recall a recent paper with single-cell sequencing that estimated phenotype and some other variables but I can't find it neither I think it is possible with RNAseq.

Although maybe I could build a classifier to actually do that (but that would be too long and I'm not sure if it is worth it)

**Notes**

1. My samples are from more than two locations.
2. My samples consist of 16S microbiome DNA-seq and RNA-seq

# Answer
**Question 1: Is there a tool for imputing phenotypes while also using genetic data?**

Yes, try Andy Dahl's Phenix from the Marchini lab.

Disclosure: I did my PhD in the Marchini lab, although I was not directly involved in this research.

**Question 2: Is there a way to impute categorical data using methods intended for continuous data?**

If the variable of interest has only two categories, then you should be able to simply recode the categories as arbitrary (differing) numbers and get sensible results. But do make sure you understand the tool you are using to see if this will work.

For more than two categories, if your imputation approach allows the simultaneous imputation of multiple variables, you could try recoding your categories as contrasts or one-hot. After recoding, you could impute the variables resulting from recoding together.
> 2 votes

---
Tags: data-management, imputation
---


# Question
Title: Convert bam file to highly compressible bam
I have a large collection of bam files and I want to post-process each of them into another bam where I can make queries about:

* the reads position and pair-endness,
* insert sizes,
* MAPQ and other flags, etc. of the reads

but where I don't need to make use of:

* the sequence string of the reads,
* the Q-scores qualities or,
* the read ids, apart from keeping pair-endness in the read id,

So I would like to know what's a good tool/parameters to process each bam and turn them into a small size bam that still retains the information needed above. CRAM format is not an option at this point as several of the tools I use are not CRAM-compatible yet.

Any ideas?

# Answer
You can just set the fields you don't need to `*`:

```
samtools view -h foo.bam | awk 'BEGIN{FS="\t"; OFS="\t"}{if($1~/^@/) {print $0} else {print "*", $2, $3, $4, $5, $6, $7, $8, $9, "*", "*"}}' | samtools view -bo smaller.bam -

```

This will set each read's name to `*`, but you can still see where its mate maps with the PNEXT and RNEXT fields.

The resulting BAM file looks like:

```
*   83  chr1    3003308 37  51M =   3003157 -202    *   *
*   99  chr1    3003364 33  51M =   3003547 234 *   *
*   99  chr1    3003410 37  50M =   3003636 269 *   *
*   99  chr1    3003438 255 50M =   3003760 373 *   *

```
> 12 votes

# Answer
I wrote http://lindenb.github.io/jvarkit/Biostar173114.html for biostars "What is the best way to make a bam file smaller by removing unwanted information?"

```
Usage: biostar173114 [options] Files

    -keepAtt, --keepAttributes
      keep Attributes
      Default: false
    -keepCigar, --keepCigar
      keep cigar : don't remove hard clip
      Default: false
    -keepName, --keepName, --name
      keep Read Name, do not try to create a shorter name
      Default: false
    -keepQuals, --keepQualities
      keep base qualities
      Default: false
    -keepRG, --keepReadGroup
      if attributes are removed, keep the RG
      Default: false
    -keepSeq, --keepSequence
      keep read sequence
      Default: false
    -mate, --mate
      keep Mate/Paired Information

```

e.g:

```
$ java -jar dist/biostar173114.jar src/test/resources/S1.bam --keepSequence | head -n 20
@HD VN:1.5  SO:coordinate
(...)
@SQ SN:RF11 LN:666
@RG ID:S1   SM:S1   LB:L1   CN:Nantes
R0  128 RF01    1   60  70M *   0   0   GGCTATTAAAGCTATACAATGGGGCCGTATAATCTAATCTTGTCAGAATATTTATCATTTATATATAACT  *
R1  64  RF01    8   60  70M *   0   0   AAAGCTATACAATGGGGAAGTATAATCTAATCTTGTCAGAATATTTATCATTTATATATAACTCACAATG  *
R2  64  RF01    11  60  70M *   0   0   GCTATACAATGGGGAAGTATAATCTAATCTTGTCAGAATATTTATCATTTATATATAACTCACAATCCGC  *
R3  64  RF01    12  60  70M *   0   0   CTATACAATGGGGAAGTATAATCTAATCTTGTCAGAATATTTATCATTTATATATAACTCACAATCCGCA  *
R4  128 RF01    27  60  70M *   0   0   GTATCATCTAATCTTGTCATAATATTTATCATATATATATAACTCACAATCCGCAGTTCAAATTCCAATA  *
R5  64  RF01    44  60  70M *   0   0   CAGAATATTTATCATTTATATATAACTCAGAATCCGCAGTTCAAATTCCAATATACTATTCTTCCAATAG  *
R6  128 RF01    44  60  70M *   0   0   CAGAATATTTATCATTTATATATAACTCACAATCCGCAGTTCAAATTCCAATATACTATTCTTCCAATAG  *

```
> 3 votes

---
Tags: sam, samtools, format-conversion, sambamba
---


# Question
Title: Salmon output: transcripts quantified with zero reads support
I quantified some samples using Salmon. According to the documentation of output format, the last column of the 'Quantification File' represents the number of reads that supported the given transcript, namely:

> NumReads — This is salmon’s estimate of the number of reads mapping to each transcript that was quantified. It is an “estimate” insofar as it is the expected number of reads that have originated from each transcript given the structure of the uniquely mapping and multi-mapping reads and the relative abundance estimates for each transcript.

In some case, I had in the output some transcripts that had `NumReads=0.000000`. Does this mean that I should discard them from the list of transcripts obtained in that sample? How comes that those transcripts end up in the output file in first place, since they don't have any reads support and the quantification was done using the reads themselves?

# Answer
> 4 votes
All transcripts should always be returned in the output, whether there's evidence supporting their expression or not. This is for the sake of simple convenience, since typically one would be merging these values across samples into a matrix. Missing values will then tend to become either `NA` or cause an error, when what would be appropriate would be to have values of 0. This also facilitates calculations such as "what's the median number of reads per transcript?" or "what percent of transcripts have an expected reads count of at least X".

---
Tags: rna-seq, quantification, salmon
---


# Question
Title: Set directory where MinKNOW writes FAST5 files
I would like to set the folder in which MinKNOW writes the raw data. How can I do this? I do not currently know where MinKNOW will output my data. What is the default directory where MinKNOW outputs the data? My OS is Mac OS X.

# Answer
> 3 votes
I found a page explaining where this is set up in a `user_conf` file here:

https://community.nanoporetech.com/support/faq/test1/minknow/minknow/how-do-i-change-the-directory-for-the-reads-folder?search\_term=user\_conf

Windows

Open Notepad with administrative rights (search for notepad, right click and then select run as administrator). Select open file and navigate to C:\ProgramFiles\OxfordNanopore\MinKNOW\conf and then open user\_conf to view the script. Under output\_dirs change the two instances of "C:\data" (default read location) to the directory required. Please note the "\\" is important.

Mac OS

Go to /Applications/MinKNOW.app/Contents/Resources/conf/user\_conf and open user\_conf to view the script. Under output\_dirs replace the default path "/Library/MinKNOW/data" with your path of choice.

Ubuntu

Go to /opt/ONT/MinKNOW/conf and open user\_conf to view the script. Under output\_dirs replace the default path "/var/lib/MinKNOW/data/reads" with your path of choice.

---
Tags: nanopore, ngs, minion, minknow
---


title: simple data subsetting on r

text:
# Question
Title: simple data subsetting on r
I have a big df like following (pseudo) df and I want to retrieve values in different col. according to the rows in the first col.

```
A B C D
a 1 1 1
a 2 3 1
a 3 3 3
a 4 5 6
b 5 6 5
b 7 6 5
c 7 6 5
c 5 4 3
c 5 5 4

```

I want to retrieve the values in col. B, C, and D according to the rows in col A and perform further analysis of them. For example, multiply the values in col B and C for all the rows of "a" in col A. Can someone suggest a commend for such a naive question. Thanks.

# Answer
Here's an approach that uses native R, using the vectorized (fast) family of `apply()` functions.

```
> df2 <- df[df$A == 'a', which(names(df) %in% c("B","C"))]
> apply(df2, 1, prod)

```

Use apply/lapply/mapply/etc. where you can in place of alternatives, where you don't know their run characteristics.

Here's a more concrete example:

```
> df
  A B C D
1 a 1 1 1
2 a 2 3 1
3 a 3 3 3
4 a 4 5 6
5 b 5 6 5
...
> df2 <- df[df$A == 'a', which(names(df) %in% c("B","C"))]
> df2
  B C
1 1 1
2 2 3
3 3 3
4 4 5
> apply(df2, 1, prod)
 1  2  3  4
 1  6  9 20

```

**Edit**: If you need to do this on each unique value of `df$A`, you can nest `apply` calls to make this scale:

```
> perA.list <- apply(data.frame(unique(df$A)), 1, function(x) apply(df[df$A == x, which(names(df) %in% c("B","C"))], 1, prod))
> perA.list
[[1]]
 1  2  3  4
 1  6  9 20

[[2]]
 5  6
30 35

[[3]]
 7  8  9
35 20 25

```

Each index in this list correspondents to the index of the value from `unique(df$A)`, i.e. `a`, `b`, and `c`.

If you don't want to do this for all unique values in `df$A`, but for some specific element in `df$A`, then you would just need to write a function. 95% of the code is in this answer to write such a function, so I'll leave that as a homework exercise for you to implement.
> 2 votes

---
Tags: r, subset
---


title: Command not found error on HISAT2

text:
# Question
Title: Command not found error on HISAT2
I have been attempting to map reads to a reference genome using HISAT2 using the Pertea, et al 2016 Nature Protocols paper. I am able to successfully generate an index (using the hisat2-build command), but when it comes to the actual mapping, I run into problems. The code I'm using (from the paper):

```
$hisat2 -p 8 --dta -x daphnia_tran -1 C2*R1*.fastq -2 C2*R2*.fastq -S C2paired.sam

```

And I get the error message:

```
./HISAT2_Basic_14.sh: line 26: -p: command not found

```

When I remove the -p argument, the error message is similar, but says "--dta: command not found"

Any ideas on what the issue is?

Thank you!

# Answer
> 3 votes
tldr: Remove the `$` from the command.

I imagine you're literally typing `$hisat2`, where you mean to instead type `hisat2`. The `$` is just meant to show the end of the command prompt. If you actually type `$hisat2`, you're not executing the hisat2 program, but whatever the `$hisat2` variable is set to. Since it's likely not set to anything, it's just ignored. For example:

```
ryan@maximus:~$ $hisat2 ls
bin              convertChromosomesChain.py  ...

```

Since there's no `$hisat2` variable, I'm actually just running `ls`. Likewise, your command is just reducing to `-p 8 --dta -x daphnia_tran -1 C2*R1*.fastq -2 C2*R2*.fastq -S C2paired.sam`, which is the source of the error.

---
Tags: sequence-alignment, transcriptome, genomics
---


# Question
Title: While opening a Bam file with the SeqAn library I get 'seqan::FileOpenError'
When I try to open a simple bam file with the SeqAn library using the following code:

```
CharString bamFileInName = getAbsolutePath("/home/omair/TFG/Files/9827_2#49.bam")
BamFileIn bamFileIn(toCString(bamFileInName));

```

I get the following error:

```
terminate called after throwing an instance of 'seqan::FileOpenError'
  what():  Could not open file /home/omair/Programs/seqan//home/omair/TFG/Files/9827_2#49.bam: iostream error

Process finished with exit code 134 (interrupted by signal 6: SIGABRT)

```

The file exists and the path is correct.

# Answer
Remove the `getAbsolutePath` call — your path is already absolute. `getAbsolutePath` garbles it, as you can see in the error message.
> 3 votes

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


title: Changing the legend names

text:
# Question
Title: Changing the legend names
```
ggplot(dataset, aes( x = Norm, fill = Condition, y=log2(measure))) + geom_boxplot() +
        scale_fill_manual(values=c("brown1","darkolivegreen4","burlywood3")) + theme_classic() +
        labs(title="Comparison",x="Norm", y = "counts") +
        scale_fill_continuous(labels = paste("condition1", "condition2", "condition3"))

```

I am trying to change the labels of the legend. Currently, it shows the labels that are present in my dataset. However, I want to change it to "condition1", "condition2", "condition3". Above is the code I am trying. It gives me following error.

```
Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale. Error: Discrete value supplied to continuous scale In addition: Warning message: Removed 18 rows containing non-finite values (stat_boxplot).

```

Could someone help me here?

# Answer
> 5 votes
You don't have any continuous fill scales, since you're filling according to `Condition`, which is a factor. This is leading to the `Error: Discrete value supplied to continuous scale` error. You just need to combine your previous `scale_fill_manual(values=...)` command with the `labels=` part of your incorrect `scale_fill_continuous()` command:

```
ggplot(dataset, aes( x = Norm, fill = Condition, y=log2(measure))) +
    geom_boxplot() +
    scale_fill_manual(values=c("brown1","darkolivegreen4","burlywood3"),
                      labels = paste("condition1", "condition2", "condition3")) +
    theme_classic() +
    labs(title="Comparison",x="Norm", y = "counts")

```

Something along those lines should work.

---
Tags: r, ggplot2
---


title: Syntenic gene browser

text:
# Question
Title: Syntenic gene browser
I found that there is a Syntenic gene browser at GEvo. Do anyone know where could I find a similar browser for a local installation?

# Answer
The UCSC genome browser supports local mirrors. You could make a local mirror that has the "TransMap Alignments" cross-species gene track. From the track description page:

> These tracks contain cDNA and gene alignments produced by the TransMap cross-species alignment algorithm from other vertebrate species in the UCSC Genome Browser. For closer evolutionary distances, the alignments are created using syntenically filtered LASTZ or BLASTZ alignment chains, resulting in a prediction of the orthologous genes in human. For more distant organisms, reciprocal best alignments are used.

They provide a helpful tutorial on setting up a local installation here.
> 1 votes

---
Tags: sequence-annotation, genome-browser
---


# Question
Title: Writing a perl script to holding information for two genes
Basically I have a perl script in which I have an array (where each element of the array references a hash) and need to be able to print the array with a dumper function. Thus I need to be able to split the $line on white spaces and save into variables and then construct anonymous hash and push into @genes:

```
#!/usr/bin/perl
use strict;
use warnings;
use Data::Dumper; # print complex data structure
my @genes; # declare the array
while(<>) { # this means that as long as lines come from the pipe we keep going
  my $line = $_; # a line that come from the pipe (we go line by line)
  next unless $line =~ /^\d+/; # skip lines except those reporting genes
  <split the $line on white spaces and save into variables>
  <anonymous hash and push into @genes>
}
print Dumper(\@genes);
exit;

```

I don't even know how to approach this, this is my first time perl scripting and I am really confused. Any insight would be extremely helpful.

So far I have:

```
my $gene = {'id' => $id, 'start' => $start, 'end' => $end, 'frame' = $frame, 'score' => $score};
my @line_split = split(' ', $line);
my $id = $line_split [0], $start = $line_split [1], $end = $line_split [2], $frame = $line_split [3], $score = $line_split [4]

```

# Answer
Perhaps this small test script will help demonstrate some of the principles:

```
#!/usr/bin/env perl

use strict;
use warnings;
use Data::Dumper;

my @arr;

# define the current line
my $line = "foo\tbar\tbaz\n";

# remove the trailing newline character
chomp $line;

# put the line's elements into an anonymous array
my ($a, $b, $c) = split('\t', $line);

# insert key-value pairs into anonymous hash
my $hash = {a => $a, b => $b, c => $c};

# add the hash table to the array
push @arr, $hash;

# print the array
print Dumper \@arr;

```

Run this script and see what it prints out.
> 6 votes

# Answer
The following Perl documentation pages should be informative:

* split \- for splitting a scalar at all matches of a defined pattern
* perlreftut \- discusses the approach of anonymous variables and how to combine them
> 2 votes

---
Tags: gene, perl, homework
---


# Question
Title: How to efficiently compute the exact percentage of non-unique k-mers in a genome for given k?
I'm looking for some software that can "efficiently" (time and memory) compute the exact percentage of non-unique k-mers in a genome for given k. I don't need the k-mers or the abundances itself, I just need the percentage.

Alternatively, the result could be:

1. the number of different k-mers observed in the genome
2. the number of unique k-mers, i.e. k-mers that occur only once in the genome, **or** the number of non-unique k-mers, i.e., k-mers that occur more than once in the genome.

which can be easily used to compute the requested percentage.

**Requested features:**

* don't count k-mers with ambiguous nucleotides
* count a k-mer and its reverse complement together
* k upto 50 or 100
* for large genomes (upto some GB)
* efficiency

For the study in mind, I like to vary k. Hence, I will run the software several times with different k.

Any help is welcome.

# Answer
You can do all of that with khmer. For example, `abundance-dist-single.py` produces a file with columns: `k-mer abundance`, `k-mer count`, `cumulative count`, and fraction of total distinct k-mers`. So for question 1 you would sum column 2. For question 2 you would just get the`k-mer count`associated with a`k-mer abundance\` of 1.

That package also provides a python API (see the read the docs link above) if you need to customize things more.
> 5 votes

---
Tags: genome, software-recommendation, k-mer
---


# Question
Title: Counting repeated kmers sequences that match at least x % of reads sequence
Working on a fastQ file, I would like to get the occurrences of repeated sequences for all possible kmers of a given length that cover at least 90% of the read's length for the whole data set.

example :

> for a length 6 with the kmer "ATTGGG" and a data set containing reads of 300 bases each

I want to obtain the occurrences of the following kmers :

```
ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGG......  (90 % of the read's length without  gaps)
ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGGA.....  (91%)
ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGGAT....  (92%)

```

this for all possible combinations of kmers of length = 6.

I have done some research and tested many softwares such as (KMC2,Jellyfish,Khmer,Dsk,Scturtle) but none matched what I would like to do, because I want to obtain the occurrences of the repeated kmer without any possible gaps in-between with a certain % of read coverage.

I've done a naive algorithm that works but the optimization isn't the best obviously. Is there a tool that does the trick ?

# Answer
The following javascript does what you want. You need node.js to run it. It should be easy to translate the code to Python. I have not carefully tested it. Use with caution.

EDIT (response to new comments): the program has been changed to compute the sum of lengths. Note that only a tandem repeat of length k\*2 or longer is counted. For example in sequence `ATTGGGATTGGGATTcGGGATTGGG`, the script returns 15, because the second is not long enough. Counting length smaller than k\*2 will be slower and more complex, but is often not the right thing to do. For example, in `ATTGGGATTGGGcccccccgaaatcgatagcatcgaGGGATTGcgatc`, we would not count the second `GGGATTG` as a tandem repeat.

On performance, for an input string of length $l$, the time complexity of this script is $O(l)$ associated with a small constant. It will be efficient and adequate to practical uses.

If you use TRF, you should run it on your sequences and then test if any tandem repeats it finds match `GGGAAT` – the procedure will be similar to my script, except that TRF will also find impure repeats.

```
// given the repeat unit length _k_, find all tandem repeats of 2*k or longer
function trf_k(k, str)
{
    var streak = 0, a = [];
    for (var i = k; i <= str.length; ++i) {
        if (i < str.length && str[i] == str[i - k]) {
            ++streak;
        } else {
            if (streak >= k)
                a.push([i - streak - k, i]);
            streak = 0;
        }
    }
    return a;
}
// find the sum of tandem repeat length with _kmer_ being the repeat unit
function trf_kmer(kmer, str)
{
    var a = trf_k(kmer.length, str), sum = 0;
    for (var i = 0; i < a.length; ++i)
        if (str.substring(a[i][0], a[i][1]).indexOf(kmer) >= 0)
            sum += a[i][1] - a[i][0];
    return sum;
}

var str = 'xATTGGGATTGGGATTGxGATTGGGATTGGGATTGGGATTGxx';
console.log(trf_kmer('GGGATT', str));

```
> 6 votes

---
Tags: fastq, k-mer, reads, repeat-elements
---


title: Intuitive explanation of GSVA analysis

text:
# Question
Title: Intuitive explanation of GSVA analysis
I'm trying to understand the way the GSVA analysis is working behind the scenes.And I was wondering if there is any way to understand it more intuitively the whole process.

So at first according to paper it *starts by evaluating whether a gene i is highly or lowly expressed in sample j in the context of the sample population distribution.* They use these kernel estimations of the cumulative density functions to transform the initial values so not to be affected by the problematic intensities.

After this "transformation" and a following normalization, GSVA calculates the enrichment scores *using the Kolmogorov-Smirnov (KS) like random walk statistic*.

As I know, Kolmogorov-Smirnov checks for differences in distributions. Which distributions does it check? Gene-set's against all the others genes? And what is the role of the random walk?

So is there any intuitive way to understand this kind of Kolmogorov-Smirnov (KS) like random walk statistic? How does it actually work? Which one is the null and which the alternative hypothesis in that case?

# Answer
> 4 votes
The most intuitive explanation is also explained in the background section:

> Conceptually, this methodology can be understood as a change in coordinate systems for gene expression data, from genes to gene sets.

Which I think it is explained in Figure 1 of the paper, which is also one of the most informative I found about the methods of an algorithm1.  I'll will use it to explain what's going on.

1. Fit a distribution function per gene

> For each gene expression profile ... a non-parametric kernel estimation ... function is performed.

I can't give more detail because I haven't understood it fully

2. Then they are normalized to make the ranks symmetric around zero

It is an inline formula between equations 2 and 3.
3. The Kolmogorov-Smirnov like random walk statistic is applied to those kernel estimations normalized distributions.

Which two distributions does it compare? The genes in a gene set and all the others (but I'm not totally sure)

The role of random walk is explained a bit later from your quote:

> \[the random walk\] produces a distribution over the genes to assess if the genes in the gene set are more likely to be found at either tail of the rank distribution

Which is later used to calculate the enrichment score by looking the maximum deviation from zero, or by summing the largest and the lowest deviations.

To address some other questions:

> GSVA calculates sample-wise gene set enrichment scores as a function of genes inside and outside the gene set, analogously to a competitive gene set test

---

1: Try finding which fittings and approximations are done inside limma

---
Tags: statistics
---


# Question
Title: Can I pilon-polish long reads with Illumina short reads to improve structural variant detection?
I have pacbio Sequel data at 50 x coverage for a strain of animal. I would like to find structural variants compared to the reference genome sequence. At the moment, I align my reads to the reference genome using NGM-LR and call the variants using Sniffles. The resulting VCF file contains variants that are well supported in a visual inspection in IGV as well as some more questionable variants (especially larger events). I have used the reads corrected by canu for alignment. I think it splits my long reads into very long reads and not-so-long reads and uses the shorter reads to correct the longer reads. However, in IGV, my reads still have a considerable number of mistakes compared to the reference. I was wondering if it is OK to use pilon to polish my long reads with some Illumina data from the same sample? The main purpose of pilon is to correct assemblies but I was wondering if it can be used to correct subreads from PacBio sequencing?

# Answer
I'd recommend using Canu for correction rather than Pilon, because it has a component that is specifically designed for read-level correction. The newest version of Canu (v1.7) will track *all* reads through the correction and assembly process:

> Ensure that every raw read is either corrected or used as evidence for correcting some other raw read. This serves to rescue short plasmids in high coverage datasets, and it is no longer necessary to set corOutCoverage to achieve the same result.

More details here.
> 2 votes

---
Tags: variant-calling, long-reads, structural-variation, pacbio, read-correction
---


title: Run gffread in multi-thread mode

text:
# Question
Title: Run gffread in multi-thread mode
Is there any possibility to run gffread in multi-thread mode? The answer seems to be 'no' from the manual (or `gffread -h`), as no multi-thread option is mentioned.

I'm mostly using this utility to extract transcript sequences (FASTA) from annotation files (GTF). In single-thread mode, runtime is acceptable in most cases, but would like to speed up a bit in some cases (especially for big genomes).

Note that I already tried, as mentioned in the manual, to index the genome first (which, again, cuts the runtime in most but not all cases):

> Note that the retrieval of the transcript sequences this way is going to be much faster if a fasta index file (genome.fa.fai in this example) is found in the same directory with the genomic fasta file.
>
> Such an index file can be created with the samtools utility prior to running gffread, like this: `samtools faidx genome.fa`
>
> Then in subsequent runs using the -g option gffread will find this FASTA index and use it to speed up the extraction of transcript sequences.

I'm also open to different alternatives than `gffread` to convert GTF -\> FASTA.

# Answer
> 3 votes
Since there doesn't seem to be an easy way to run this in parallel, you could instead break the job into sections. For example, separate each chromosome into its own gtf file, extract the sequences using that file and then cat them all together.

The commands would be something like this (using the human hg38 genome and the gff GENCODE annotations downloaded from here):

Split the file into one file per chromosome (or whatever else you have as the first field of your gff/gtf file):

```
tmp=$(mktemp -d); ## create a temp directory
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' gencode.v27.annotation.gtf

```

Extract the sequences described in each file, running each command in the background (`&`):

```
for f in $tmp/*; do
    gffread -w "$f.fa" -g hg38.fa "$f" &
done

```

Concatenate them into one file:

```
cat "$tmp/*fa" > all.fa

```

You can combine all this into a single command, using `wait` to make sure you don't concatenate until all of the sub-commands have finished:

```
tmp=$(mktemp -d);
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' gencode.v27.annotation.gtf &&
    for f in "$tmp"/*; do
        gffread-0.9.12.Linux_x86_64/gffread -w $f.fa -g hg38.fa "$f" &
    done
    wait
    cat "$tmp"/*fa > all.fa

```

Or, just save it as a script called pgffread.sh (or whatever):

```
#!/bin/sh
gtf="$1"
genome="$2"

tmp=$(mktemp -d);
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' "$gtf" &&
    for f in "$tmp"/*; do
      gffread -w "$f".fa -g "$genome" "$f" &
    done
wait
cat "$tmp"/*fa
rm -rf "$tmp"

```

And then run it giving the genome and annotation files as parameters:

```
pgffread.sh gencode.v27.annotation.gtf hg38.fa > all.fa

```

On the system I ran it on, which has 48 cores, the time difference was significant (although it's so fast anyway, I really don't know if it is worth the effort unless you do this very often):

```
$ time gffread -w out.fa -g hg38.fa gencode.v27.annotation.gtf
real    1m52.072s
user    1m3.160s
sys     0m47.257s

$ time ./pgffread.sh gencode.v27.annotation.gtf hg38.fa > all.fa
real    0m17.533s
user    1m5.228s
sys     0m43.039s

```

---
Tags: rna-seq, fasta, file-formats, gtf
---


# Question
Title: Can I use DRIMseq to model interaction effect of variables?
I am trying to use DRIMseq for DTU with 2 treatments on two different strains of mice. In my previous analysis with genes only, I had found a strong interaction effect of strain and treatments. Is it possible to model that with DRIMseq?

Also when I try to add strain term and use it as the batch effect as described in the tutorial, I get the following error:

```
! Using a subset of 0.1 genes to estimate common precision !
Error in optimHess(par = par, fn = dm_lik_regG, gr = dm_score_regG, x = x,  :
  non-finite value supplied by optim
In addition: There were 50 or more warnings (use warnings() to see the first 50)

```

Please see the sessionInfo() below.

```
R version 3.4.0 (2017-04-21)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_GB.UTF-8     LC_MONETARY=en_US.UTF-8
 [6] LC_MESSAGES=en_GB.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] DRIMSeq_1.6.0                            pheatmap_1.0.8                           ggplot2_2.2.1
 [4] edgeR_3.16.5                             limma_3.30.13                            TxDb.Mmusculus.UCSC.mm10.knownGene_3.4.0
 [7] GenomicFeatures_1.26.4                   AnnotationDbi_1.36.2                     Biobase_2.34.0
[10] GenomicRanges_1.26.4                     GenomeInfoDb_1.10.3                      IRanges_2.8.2
[13] S4Vectors_0.12.2                         BiocGenerics_0.20.0                      tximport_1.2.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.12               RColorBrewer_1.1-2         plyr_1.8.4                 pillar_1.2.1               compiler_3.4.0
 [6] XVector_0.14.1             bitops_1.0-6               tools_3.4.0                zlibbioc_1.20.0            biomaRt_2.30.0
[11] digest_0.6.12              bit_1.1-12                 gtable_0.2.0               RSQLite_2.0                memoise_1.1.0
[16] tibble_1.4.2               lattice_0.20-35            pkgconfig_2.0.1            rlang_0.2.0                Matrix_1.2-10
[21] DBI_0.7                    stringr_1.3.0              rtracklayer_1.34.2         Biostrings_2.42.1          locfit_1.5-9.1
[26] bit64_0.9-7                grid_3.4.0                 XML_3.98-1.10              BiocParallel_1.8.2         magrittr_1.5
[31] reshape2_1.4.3             blob_1.1.0                 MASS_7.3-47                scales_0.5.0               Rsamtools_1.26.2
[36] GenomicAlignments_1.10.1   SummarizedExperiment_1.4.0 colorspace_1.3-2           labeling_0.3               stringi_1.1.7
[41] lazyeval_0.2.1             munsell_0.4.3              RCurl_1.95-4.10

```

# Answer
One way to see the interaction effect would be including it in the design matrix as explained in the limma tutorial. Something on the line of:

```
design_full <- model.matrix(~ treatment*strain, data = samples(d))

```

Or you could build it yourself in a matrix like

```
 (Intercept) strainA_treatmentA strainA_treatmentB  strainB_treatmentA strainB_treatmentB
       1                      0                  1                  0                   0
       1                      1                  0                  0                   0

```

About the error, it seems like the data is not in the right format, you have either `NA` or `Inf` somewhere in the matrix. Check it and rethink the processing you made.
> 2 votes

---
Tags: differential-expression
---


# Question
Title: How are positions visually indicated in the UCSC Genome Browser?
I'm a little confused about how the position indices work in the browser. Based on the picture below, is base number 755440 an `A` or a `G`?

The general question is how do the indexing guides work? Is the base right under the `|` the one corresponding to the index or is it the one just before the `|`?

# Answer
It's an A. The number refers to the position just before the `|`. This is clearer on different zoom levels, but I double checked by going to the very beginning of a chromosome:

As you can see in the image above, the `|` is always after the base it refers to. I also confirmed by checking a few random positions using another genome browser which places the positions unambiguously above each base.

For more details, see this UCSC blog post.
> 5 votes

---
Tags: ucsc, genome-browser
---


# Question
Title: Pepsin digest (cleavage) does not work using RE?
Aim is to code the theoretical peptic cleavage of protein sequences in Python. The cleavage (cutting) rule for pepsin is: `1234^567`. This is a position, `^` stands for the cleavage point.

There are rules that it mustn't have at all and rules where it must have at least one. The rule for mustn't is in position 2, it mustn't be R, K or H. In position 3 and 7, there can't be a P and in position 4, there mustn't be an R.

If it follows "the must not" rules, then it will check "the must" rules so that if any of one them is met, it will cleave. So if there is a P in position 1 or 2, or if there's an L, W, Y or F in position 4 or 5 then, it will be cleaved (use newline).

This is the code I have so far, but it doesn't work yet:

```
import re
print('INSERT YOUR PEPTIDE')
X= input();
r=re.findall('(?<=[P])(?<=[P,^R,^K,^H])(?<=([^P])(?<=[L,W,F,Y,^R])(?<=[L,W,Y,F])/D(?<=[^P])',X)
print(r)

```

My input is:

```
TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI

```

How can I make this work?

# Answer
I am really not sure if I've understood your rules correctly, but I think you want this:

* If all of these are true:

+ position 2: not R, K or H
    + position 3: not P
    + position 4: not R
    + position 7: not P
* And, if at least one of these is true:

+ position 1: P
    + position 1: P
    + position 4: L or W or Y or F
    + position 5: L or W or Y or F

Then cleave.

If so, then you probably don't want to do it using regular expressions. First because it will be hard, and second because it will be absolutely awful to maintain if you need to come back to it. Instead, try something like this:

```
#!/usr/bin/env python3
# X = input();
X = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
seq = list(X)
for i, res in enumerate(seq):
    if i > len(seq) - 7:
        break
    if (
            seq[i+1] in ('R','K','H') or
            seq[i+2] == 'P' or
            seq[i+3] == 'R' or
            seq[i+6] == 'P'
    ):
        continue
    if res == 'P':
        start = i
    elif seq[i+1] == 'P':
        start = i+1
    else:
        continue
    if (
            seq[i+4] in ('L', 'W', 'Y', 'F') or
            seq[i+5] in ('L', 'W', 'Y', 'F')
    ):
        print("%s\n%s\n%s" %
            ("".join(seq[0:start]),
             "".join(seq[start:start+7]),
             "".join(seq[start+7:])
            )
        )

```

This will produce the following output when run:

```
$ foo.py
TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMV
PKCFNEL
LMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI

```

---

If you insist on using regular expressions, you could separate it into two expressions and use something like this:

```
import re
X = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
re1='P[^RKH][^P][^R][LWYF].[^P]'
re2='P[^P][^R][LWYF].[^P]'
r=re.findall('(?:P[^RKH][^P][^R][LWYF].[^P])|(?:P[^P][^R][LWYF].[^P])',X)
print(r)

```

Before you do, however, please consider the wise words of XKCD 208:
> 2 votes

# Answer
Ok I think I almost solve this problem so I still using RE

```
import re
myprotein = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
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)

```

and it produce to

```
TTG
W
M
F
KG
LY
Q
L
KSDEARMHN
LQRW
YV
F
KVIE
F
W
GINRKDQHDSRY
WF
SH
F
CW
K
L
EHGMVPKC
FNE
L
LMPPMSTHRSMF
NA
Y
SMS
F
CVQVKE
WRGSISRMCDHHQQQCQ
Y
W
RL
MAGGT
LI

```

but the true answer should be like (checking from https://web.expasy.org/peptide\_cutter/ )

```
TTG
W
M
F
KG
LY
Q
L
KSDEARMHN
LQRW
YV
F
KVIE
F
W
GINRKDQHDSRY
WF
SH
F
CW
K
L
EHGMVPKC
FNE
L
LMPPMSTHRSMF
NA
Y
SMS
F
CVQVKE
WRGSISRMCDHHQQQCQ
Y
W
RL
MAGGT
L
I

```

so how can I make this work?
> 2 votes

---
Tags: python
---


title: Digisome generation

text:
# Question
Title: Digisome generation
I scoured and searched the web for some implementation of digisome creation, but I had a hard time finding any good source code.

By digisomes, what I mean is images with chromosomal masks/borders that have been simulated. I am seeking this because I am trying to do a pixel level prediction on chromosomes using a deep network. And, I need ground truths with the chromosomes being illuminated. So ideally, I would like a binary image with just the chromosomes being illuminated.

How exactly are digisomes generated? Is a single actual chromosome extracted from some metaphase image, then morphed in interesting ways, and finally, all the different morphs are stitched together in one big collage format so that it there are numerous chromosomes all seemingly different?

Also, if someone could point me to some good source code, that'd be helpful.

# Answer
> 2 votes
I found a very nice github repo that simulates chromosomes and their labels for overlapping chromosomes - and with a little bit of tweaking, it can be amended for a variety of different simulations.

This was also used in a Kaggle post; so right up the alley for deep learning.

---
Tags: database
---


# Question
Title: Hierarchial PCA Clustering with duplicated row names
I have a matrix object in R of 1500 rows and 20 columns. Based on my own algorithm of analysis, I know there are roughly 7 clusters. In my matrix object, each row is one data point, and is labeled to which cluster it belongs by its row name. Thus, there are duplicated row names. I wanted to see if PCA would give me the approximate same clustering of my data.

```
Segment1 Segment2 Segment3 Segment4 Segment5 Segment6 Segment7 Segment8 Segment9 Segment10 Segment11 Segment12 Segment13 Segment14 Segment15 Segment16 Segment17 Segment18 Segment19 Segment20
B  38.40000 41.75200 44.38400 44.18400 45.37600 37.49600 41.36800 33.93600 38.00800  42.51200 46.49600  40.48000  45.40800  46.32800  43.78400  39.88800  38.84000  40.56800  42.03200 38.89185
C  45.53846 50.08462 39.91538 36.95385 34.96154 39.74615 38.01538 35.75385 35.54615  36.69231 35.20769  38.05385  39.29231  37.96923  37.30000  36.86923  39.19231  38.81538  43.69231 38.06400
C  46.05176 41.69412 38.80000 37.75529 39.67529 39.07765 39.17647 38.24941 39.58588  38.63529 38.30588  41.87765  38.97412  40.13647  42.27294  38.24471  35.41647  40.80000  38.07059 42.11294
B  44.20000 43.42857 44.80000 35.20000 35.91429 37.82857 51.45714 44.68571 46.68571  48.74286 41.25091  39.45455  38.17091  40.70182  40.39273  41.28727  40.63636  41.50909  41.68364 41.29455
E1 45.06909 41.09818 40.02909 42.50182 42.34909 39.84727 41.42909 40.47273 40.28000  40.51636 41.25091  39.45455  38.17091  40.70182  40.39273  41.28727  40.63636  41.50909  41.68364 41.29455
E3 40.87407 39.27704 44.13630 43.25037 35.86667 37.30667 38.76148 40.74667 38.93333  43.16148 37.47259  37.73630  38.34370  39.00148  36.96889  37.76593  39.14667  37.92593  37.62963 38.89185

```

(The first column are the row names, either being A,B,C,D,E1,E2 or E3 according to my clustering)

I used various packages and functions like **pcomp()**, **princomp()** from the **stats** package, and **PCA()** from the **FactoMineR** package. If I wanted to do hierarchical clustering with the **HCPC()** function from the same package, it was not possible because of the duplicated row names.

Using:

```
fviz_pca_ind(pca1, geom.ind = "point", axes = c(1,3) , col.ind = new$Class, palette = c("#66C2A5" ,"#FC8D62" ,"#8DA0CB", "#E78AC3" ,"#A6D854", "#FFD92F", "#E5C494"), addEllipses = T, legend.title="Classes")

```

I got the individual datapoints on my two PC axes with ellipses around them.

But for clustering:

```
PCA_FR <- PCA(mat_data, ncp = 4, graph = FALSE)
Hier_PCA<- HCPC(PCA_FR, graph = F)

```

I couldn't compute Hier\_PCA because of duplicated row names.

```
> Hier_PCA<- HCPC(PCA_FR, graph = F)
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
In addition: Warning messages:
1: In data.row.names(row.names, rowsi, i) :
some row.names duplicated:  3,4,7,8,10,11,12,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,

```

How can I best cluster my data using the first three PCs, and visualize the 1500 datapoints, but labeled only in 7 ways?

# Answer
> 2 votes
Try to make your rownames unique first.

So if you have this matrix with duplicate rownames:

```
m <- replicate(20, rnorm(6))
rownames(m) <- c("B", "C", "C", "B", "E1", "E3")
m
         [,1]        [,2]      [,3]      [,4]       [,5]        [,6]       [,7]
B   2.4071997  1.30040528 1.0765500 0.2005440  1.1732384  2.34462132  0.3702172
C  -0.2683534 -0.74607893 0.1611073 0.3145887  0.8208741  1.26575765  0.4205310
C   0.1997264  0.03318396 0.1967151 0.7785350  0.9382146 -0.32475195 -0.9832998
B  -1.2726538  0.22519992 1.1556121 0.7161238  0.5601368  1.02900724 -0.2726709
E1  0.8062722 -1.80899128 0.3936155 0.5724118 -0.8163454  0.04306317 -0.4074828
E3  0.3320706  0.21218846 0.9210126 0.9151837  0.4588621 -1.27784583  1.3292399
         [,8]       [,9]       [,10]      [,11]       [,12]       [,13]
B   0.2305559  0.2964808 -0.01035426 -2.6062776 -0.03143381 -1.48828554
C   0.1622548 -0.6244459 -0.66571655 -0.3529994  1.07540375 -1.58949697
C   0.5984752 -0.6423845 -1.00001700  0.2212160  0.05499773  1.36182121
B  -0.5594747 -0.3570748  1.44098577 -0.4939490 -0.93386684  0.54919399
E1  1.6378047  0.8702920  3.06641490  0.1948452  0.36339348 -1.04694826
E3 -2.2422997  1.0682864  0.14326146 -0.4513848 -0.25819231  0.01364221
        [,14]      [,15]      [,16]      [,17]       [,18]      [,19]
B  -2.0410367  2.4858879  1.0564085 -1.7353230 -0.63024160 -1.2642884
C  -1.2216387  0.5533529 -0.2757584 -0.2119105 -0.68899706  1.0587382
C   0.6401224 -1.5998479  1.4072007 -1.5062900 -0.60161624  0.9781313
B   0.4309208  1.4310142 -0.1702066 -0.7195977  0.37466185  0.2605195
E1  1.2538412  0.3413762 -1.9346107  0.7314531 -0.03105008  0.7341104
E3 -0.3146694 -1.3028086  0.5043235 -0.2568268  0.08336845  1.3097557
        [,20]
B  -1.1245496
C  -0.4552617
C  -1.5498408
B   0.1261818
E1 -0.1180164
E3  2.1035368

# Make them unique with make.names function

rownames(m) <- make.names(rownames(m), unique=T)

m
          [,1]        [,2]      [,3]      [,4]       [,5]        [,6]
B    2.4071997  1.30040528 1.0765500 0.2005440  1.1732384  2.34462132
C   -0.2683534 -0.74607893 0.1611073 0.3145887  0.8208741  1.26575765
C.1  0.1997264  0.03318396 0.1967151 0.7785350  0.9382146 -0.32475195
B.1 -1.2726538  0.22519992 1.1556121 0.7161238  0.5601368  1.02900724
E1   0.8062722 -1.80899128 0.3936155 0.5724118 -0.8163454  0.04306317
E3   0.3320706  0.21218846 0.9210126 0.9151837  0.4588621 -1.27784583
          [,7]       [,8]       [,9]       [,10]      [,11]       [,12]
B    0.3702172  0.2305559  0.2964808 -0.01035426 -2.6062776 -0.03143381
C    0.4205310  0.1622548 -0.6244459 -0.66571655 -0.3529994  1.07540375
C.1 -0.9832998  0.5984752 -0.6423845 -1.00001700  0.2212160  0.05499773
B.1 -0.2726709 -0.5594747 -0.3570748  1.44098577 -0.4939490 -0.93386684
E1  -0.4074828  1.6378047  0.8702920  3.06641490  0.1948452  0.36339348
E3   1.3292399 -2.2422997  1.0682864  0.14326146 -0.4513848 -0.25819231
          [,13]      [,14]      [,15]      [,16]      [,17]       [,18]
B   -1.48828554 -2.0410367  2.4858879  1.0564085 -1.7353230 -0.63024160
C   -1.58949697 -1.2216387  0.5533529 -0.2757584 -0.2119105 -0.68899706
C.1  1.36182121  0.6401224 -1.5998479  1.4072007 -1.5062900 -0.60161624
B.1  0.54919399  0.4309208  1.4310142 -0.1702066 -0.7195977  0.37466185
E1  -1.04694826  1.2538412  0.3413762 -1.9346107  0.7314531 -0.03105008
E3   0.01364221 -0.3146694 -1.3028086  0.5043235 -0.2568268  0.08336845
         [,19]      [,20]
B   -1.2642884 -1.1245496
C    1.0587382 -0.4552617
C.1  0.9781313 -1.5498408
B.1  0.2605195  0.1261818
E1   0.7341104 -0.1180164
E3   1.3097557  2.1035368

```

---
Tags: r, pca
---


# Question
Title: Extracting strand specific reads from MinION cDNA-PCR protocol
I recently performed my first MinION run, and now I'm trying to analyze the data. Being pretty new to the bioinformatics field, I was hoping some of you could help me out.

As a bit of background, I'm trying to define transcript isoforms, and down the line, I'm hoping to quantify isoform abundances. As my first attempt, I used a 1D kit that first performs RT starting at the polyA tail and then template switching to enrich for full length transcripts. To initiate RT, a poly-dT primer is added and after template switching a second primer is added. Using those added primer sequences, I then performed PCR to amplify the full-length cDNAs. This seemed to work well, so I went ahead and sequenced the library.

Now, on to my problem. I have aligned my fastq file to my genome with both bwa-mem and minimap2. In both cases, I get alignment of reads to the genome where I expect transcription to be occurring, but half of the reads are in the sense direction and half are in the antisense direction. This is presumably due to the fact that the cDNA is double stranded and either strand can be recruited to pore to be sequenced. However, because I'm interested in mapping transcript isoforms which could include anti-sense transcripts, this isn't exactly what I hoped for.

I'd like to figure out a way to extract only the reads that correspond to the actual mRNA transcript's sequence and not the reverse complement. Does anyone have ideas as far as programs that would help with this?
I'm thinking that I should be able to extract out only the molecules that have the polyA or have the strand switching primer sequence or something like that, but I'm not sure how to properly implement this.

# Answer
> 3 votes
This has been set up as a bioinformatics protocol on protocols.io.

I've been using LAST to identify the adapter orientation relative to the genome, and then using that information to split BAM files up and recombine them to create two strand-specific files that are displayable in a genome browser.

I have written my own script to process LAST results into a CSV format, which makes it easier to do line-by-line data filtering.

*Note:* these scripts have been slightly modified from scripts that I have run. Consider them demonstrative: pay attention to the words rather than the script.

# Read correction

I prefer starting off my data analysis with a read correction with Canu (ideally v1.7 now that it's out, because that attempts correction of *all* reads, but here I use canu v1.6). I use `minimap` as the mapper to speed this up. The `genomeSize` parameter should be approximately a tenth to a fortieth of the number of bases in your dataset to make sure that no sequences are excluded (bigger is better, as long as Canu doesn't freak out about memory consumption):

```
~/install/canu/canu-1.6/Linux-amd64/bin/canu overlapper=minimap \
  genomeSize=100M minReadLength=100 minOverlapLength=30 -correct \
  -p 4T1_BC06 -d 4T1_BC06 -nanopore-raw \
  workspace/pass/barcode06/fastq_runid_*.fastq

```

This creates a file `4T1_BC06/4T1_BC06.correctedReads.fasta.gz`.

# Chimeric read filtering

The next step I carry out is a basic read-level QC to exclude chimeric reads. Porechop can be used for this, although that removes adapters by default, which is not particularly useful in this case.

I use LAST to search for adapter sequences within the corrected reads, pass it through my conversion script, and extract out duplicated mappings (i.e. where the same read/adapter pair appears more than once in the mapping results):

```
lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | awk -F',' '{print $1,$2}' | sort | \
  uniq -d | awk '{print $1}' | uniq > reads_with_duplicated_adapters.txt

```

I use this file to filter out chimeric reads from the corrected dataset using another fastq filtering script I've created:

```
~/scripts/fastx-fetch.pl -v -i reads_with_duplicated_adapters.txt 4T1_BC06/4T1_BC06.correctedReads.fasta.gz | \
  gzip > 4T1_BC06.correctedReads.uniqueOnly.fasta.gz

```

# Adapter filtering

Along roughly similar lines to the chimeric read filtering, I then look for the strand switch and VNT adapters in the sequences. For a forward-orientated query, I expect the strand switch primer to be in the forward direction, and the VNP primer to be in the reverse direction:

```
lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,+' -e 'ONT_VNP,-' | \
  awk -F',' '{print $1}' | sort | uniq > fwdQry_seqs_BC06.txt
lastal -P 10 barcodes_primerSeqs.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,-' -e 'ONT_VNP,+' | \
  awk -F',' '{print $1}' | sort | uniq > revQry_seqs_BC06.txt

```

And then more filtering to split the reads up into forward and reverse \[genome-relative\] subsets:

```
~/scripts/fastx-fetch.pl -i fwdQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz
~/scripts/fastx-fetch.pl -i revQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz

```

# Mapping

Now that the reads have been oriented, the mapping can be done:

```
~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > fwd_4T1_BC06.CU_vs_mmus.bam
~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > rev_4T1_BC06.CU_vs_mmus.bam

```

# Splitting and recombining

The mapped files are then split based on their mapping direction using `samtools view` flag filtering to exclude or include reverse-mapped reads:

```
samtools view -F 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > fwd_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > rev_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -F 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > rev_rev_4T1_BC06.CU_vs_mmus.bam

```

And finally the files are recombined to identify forward and reverse-encoded transcripts via `samtools merge`. I'm calling the "pos" strand the one that is encoded in the same direction as the genome, and the "neg" strand the reverse-complement direction:

```
samtools merge pos_4T1_BC06.CU_vs_mmus.bam \
  rev_fwd_4T1_BC06.CU_vs_mmus.bam fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools merge neg_4T1_BC06.CU_vs_mmus.bam \
  fwd_fwd_4T1_BC06.CU_vs_mmus.bam rev_rev_4T1_BC06.CU_vs_mmus.bam

```

Alternatively, it's possible to set the BAM flags for read 1 and read 2 for the reverse and forward-encoded adapters respectively, which allows everything to then be combined into a single BAM file (and treated in the same way as a strand-specific Illumina BAM file).

If this has worked properly, then mapping human or mouse to the mitochondrial genome should show most expression appearing on the positive strand, with a small scattering of negative-strand expression, a bit like this:

For a bit of a background on how/why this works, see this answer.

---
Tags: sequence-alignment, nanopore, minion
---


# Question
Title: How to read gene regulatory network edge list files?
The following is an excerpt from an edge list file from the Gene Regulatory Database, YeastNet v3

```
YML100W YMR261C 5.73968189451551
YDR074W YML100W 5.73143623470086
YDR074W YMR261C 5.72385204041689
YML106W YMR271C 5.64698433340441
YGR240C YMR205C 5.64063955646679
YDL131W YDL182W 5.63907349642605
YNL037C YOR136W 5.62782704009976
YHR208W YJR148W 5.61712512462848
YAL062W YOR375C 5.59956647235823
YNL104C YOR108W 5.59046546303091
YGR132C YGR231C 5.58996626839351
YBR140C YOL081W 5.58738025271142
YDR453C YML028W 5.58281825879471
YNL096C YOR096W 5.58060076040965
YGR180C YJL026W 5.57785446902287
YGR056W YLR357W 5.57582611528462
YBR009C YNL030W 5.57261292043205
YMR169C YMR170C 5.55866364259678
YAR071W YBR093C 5.55668852247789
YER081W YIL074C 5.55252000887001

```

What do these rows mean?

# Answer
> 3 votes
You can read about it on the website, the last column is the Log likelihood scores (LLS) between the genes in the first two columns. An LLS of 0 means not better than random link.

On the website they show papers with fuller details about the statistics.

---
Tags: gene, networks, transcriptome-regulation
---


# Question
Title: Merging regions according to their identifier
I have a huge file (20 `GB`) which has a range of genomic locations, and for each location there is an identifier(`4th column`), which is sometimes the same.

```
file1.txt

chr1  10  20  ABC
chr1  13  20  ABC
chr1  14  21  ABC
chr1  22  27  ABC
chr1  29  37  ABC
chr2  15  21  JJJ
chr2  21  31  JJJ
chr2  23  27  JJJ
chr2  35  56  JJJ
chr2  25  26  MMM
chr3  32  42  MMM
chr3  45  76  MMM
chr3  88  101 MMM
chr3  101 105 MMM

```

I am trying to merge the ranges of column `2` and `3`, if they are in the same chromosome and having the same identifier ( having the `1st` and `4th` columns the same).

I have tried to use `bedtools`, as shown in this post, which looks very similar but not the same with mine, as they are looking for different outputs ( range vs merge )

So the answers in that post was with using `groupby` in bedtools2,source code from `github`.

I have applied the same as the following:

```
./groupBy -i ../file1.txt -g 1,4 -c 2,3 -o min,max | awk -v OFS='\t' '{print $1, $3, $4, $2}'

chr1    10      37      ABC
chr2    15      56      JJJ
chr2    25      26      MMM
chr3    32      105     MMM

```

But it only connects the ranges of identifiers.. It groups by, and does not merge.

And, if I use `bedtools v2.26 merge` function :

```
>sort -k1,1 -k2,2n file1.txt > file2.txt
>cat file2.txt
chr1    10      20      ABC
chr1    13      20      ABC
chr1    14      21      ABC
chr1    22      27      ABC
chr1    29      37      ABC
chr2    15      21      JJJ
chr2    21      31      JJJ
chr2    23      27      JJJ
chr2    25      26      MMM
chr2    35      56      JJJ
chr3    32      42      MMM
chr3    45      76      MMM
chr3    88      101     MMM
chr3    101     105     MMM

bedtools merge -i file2.txt
chr1    10      21
chr1    22      27
chr1    29      37
chr2    15      31
chr2    35      56
chr3    32      42
chr3    45      76
chr3    88      105

```

Which, merges without seeing the identifiers as it also gives me this output.

What I want is to `merge` the ranges if they `overlap`, if having the `same 1st and 4th column` as:

```
expected_output.txt

chr1  10  21  ABC
chr1  22  27  ABC
chr1  29  37  ABC
chr2  15  31  JJJ
chr2  35  56  JJJ
chr2  25  26  MMM
chr3  32  42  MMM
chr3  45  76  MMM
chr3  88  105 MMM

```

I only have access to 15GB of RAM, so I can't store the entire file in memory.

# Answer
The following python script will do what you want and should be relatively memory efficient. It processes a single chromosome at a time, so either sort the input or ensure that at least entries in a group are next to each other.

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

# Process a sorted list of tuples (e.g., [(0, 10), (10, 20), ...])
# Overlapping items will be merged and a list of the merged regions is
# returned (e.g., [(0, 20), ...])
def merge(values):
    start, end = values[0]
    o = []
    for v in values[1:]:
        if v[0] <= end:
            if v[1] > end:
                end = v[1]
        else:
            o.append((start, end))
            start, end = v
    o.append((start, end))
    return o

# Process a single chromosome, one group at a time
def processChrom(d, chrom):
    for k, v in d.items():
        v.sort()
        tuples = merge(v)
        for t in tuples:
            print("{}\t{}\t{}\t{}".format(chrom, t[0], t[1], k))

def main():
    currentChrom = None
    d = dict()
    for line in sys.stdin:
        cols = line.split()
        if cols[0] == currentChrom:
            if cols[3] not in d:
                d[cols[3]] = []
            d[cols[3]].append((int(cols[1]), int(cols[2])))
        else:
            if currentChrom:
                processChrom(d, currentChrom)
            d = dict()
            d[cols[3]] = []
            d[cols[3]].append((int(cols[1]), int(cols[2])))
            currentChrom = cols[0]

    if currentChrom:
        processChrom(d, currentChrom)

if __name__ == "__main__":
    main()

```

Usage is `python foo.py < foo.txt` and the output is:

```
chr1    10  21  ABC
chr1    22  27  ABC
chr1    29  37  ABC
chr2    25  26  MMM
chr2    15  31  JJJ
chr2    35  56  JJJ
chr3    32  42  MMM
chr3    45  76  MMM
chr3    88  105 MMM

```

You'll notice that this isn't necessarily sorted, so you may want to pass the output through `bedtools sort` or `sort`.
> 3 votes

# Answer
If I understand correctly, you could create dummy chromosomes made by merging chromosome and identifier, merge with bedtools, split back chromosomes and identifiers. E.g.

```
awk -v OFS='\t' '{print $1"_"$4, $2, $3}' file1.txt \
| mergeBed \
| sed 's/_/\t/' \
| awk -v OFS='\t' '{print $1, $3, $4, $2}'

```

(Assuming "\_" is not part of chromosome or identifier name, adjust accordingly if so).
> 6 votes

# Answer
Just to add this can be done with `bedtools` by processing each identifier separately which may or may not be practical given your large file.

```
for ID in $(perl -ane 'print $F[3] . "\n";' file1.txt | sort | uniq);do
    grep $ID file1.txt | \
    perl -ane 'chomp;print join("\t",@F[0..$#F]) . "\n";' | \
    sort -V | \
    bedtools merge -i stdin -c 4 -o distinct
done | sort -V

```
> 4 votes

# Answer
Here is a fairly simple (and hopefully readable) native Python solution. It assumes that the bed file is sorted prior to parsing:

```
from itertools import groupby

class BedRecord(object):
    def __init__(self, chr, start, end, feature):
        self.chr = chr
        self.start = start
        self.end = end
        self.feature = feature

def parse_records(in_file):
    with open(in_file) as f:
        for line in f:
            yield BedRecord(*line.strip().split())

def merge_groups(group):
    output = [next(group)]
    for record in group:
        if record.start <= output[-1].end:
            output[-1].end = max(output[-1].end, record.end)
        else:
            output.append(record)
    return output

def main(in_file):
    for _, group in groupby(parse_records(in_file), lambda x: (x.chr, x.feature)):
        for r in merge_groups(group):
            print('\t'.join([r.chr, r.start, r.end, r.feature]))

if __name__ == '__main__':
    main('test.bed')

```
> 3 votes

---
Tags: bedtools, merge
---


# Question
Title: LRT or LRT-like test on cyclical (Sleep) data
I have RNA-Seq data from 4 time points (3 hours awake, 9 hours awake, 3 hours asleep, 9 hours asleep). I'm interested in doing something similar to a LRT where genes are found to be significant if they changed at any point in the time course.

However, I'm not sure how to do this considering this is cyclical. Running a LRT test in DESeq2 gives different results depending on which timepoint I choose as the beginning of the time course since I assume it doesn't calculate the last transition back to being awake (Or asleep depending on your reference).
Any ideas on an appropriate way to do this?

```
time = c(3,9,15,21)
ddsLRT <- DESeqDataSetFromTximport(txi.rsem.lrt, sampleTable, ~time)
ddsLRT<- DESeq(ddsLRT, test="LRT", reduced = ~1)

```

# Answer
> 5 votes
I think the problem is that you have `time` as a linear value rather than a factor. While this naively makes sense (after all, you have to go through 3 hours awake before you can get to 9 hours awake), what you end up doing is fitting a single coefficient to each gene of "linearly changes with time". What you want instead is factor, so you can have a system where things like `per` and `clk` can cycle:

```
time = factor(c(3,9,15,21))

```

The rest is then correct. Note that LRTs are not asking, "is there a change at any time point?", since there will likely be DE genes according to this that are not DE at any single time point. Rather an LRT such as this is asking, "What are genes that vary significantly over time?" There's a subtle different between those two questions. The former question is actually answered by a bunch of pairwise comparisons (taking the union of the DE genes).

---
Tags: rna-seq, r, statistics, deseq2
---


# Question
Title: Are there certain alignment methods/tools which perform better with a high density of indels?
I have a set of experiments which should result in WGS reads with a high density of indels.

Question: Are there certain alignment tools/methods which perform respectively "better" with these reads? I'm not very interested in performance, but rather "better" in terms of dealing with the presence of many deletions/insertions.

Options I know of:

* BWA (I think there are special settings whereby I could run bwa for this case, e.g. see this)
* BBMap

# Answer
> 4 votes
Minimap2 attempts to solve a lot of problems that are not properly dealt with by other mapping software. It's written by the same author as BWA, but designed to deal with situations where the reference and query sequences are quite different (e.g. for noisy long reads).

---
Tags: sequence-alignment
---


# Question
Title: How to download gene expression data from NCBI gene database
In the NCBI gene database, I can add the expression tracks (circled in picture blow) through 'Tracks' button, but **How I can download the expression data directly, not just look the picture?**

# Answer
You can download raw data from SRA (NCBI). On this website under RNA-Seq alignments, you'll find the samples. Click on them, for example SAMN05231885, and then click on PRJNA325427 to see the link to the SRA data (click on ). Here you can download the raw fastq files of the experiment.
> 3 votes

---
Tags: rna-seq, database, gene, data-download
---


# Question
Title: How to get results from Homo.sapiens package in bioconductor for a specific reference
I want to use the Homo.sapiens package in Bioconductor to retrieve the chromosome location start and end for each gene symbol in a specific reference (e.g. hg19 or hg38).

Right now, I am using the following code:

```
> select(Homo.sapiens, keys=c('KRAS'), cols <-c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'), keytype='SYMBOL')

  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25358180 25403854
2   KRAS   chr12 25386768 25403863

```

My question is two-fold:

* How can I specify the human genome reference I want to retrieve the data from? (e.g. hg19 or hg38)
* Why do I receive multiple and different mappings through `Homo.sapiens` for some genes? If I check the locations through GeneCards I get:
    + `chr12:25,204,789-25,250,936 (GRCh38/hg38)`
    + `chr12:25,357,723-25,403,870 (GRCh37/hg19)`

# Answer
I couldn't find how is this package build, my guess is that is as in the example of the OrganismDBI package. See that vignette section to build your own data package with the genome you want.

The locations are different because you ask for the transcript start and end, if you add the "TXNAME", you will (probably) see there are different transcripts, so different mappings are reported (untested).

### Update

The example of the OrganismDBI package builds correctly the Homo.sapiens package, which appears to be linked with hg19. I provide here the relevant code to build the package for hg19 and hg38.

```
> source("https://bioconductor.org/biocLite.R")
> biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
> biocLite("TxDb.Hsapiens.UCSC.hg38.knownGene")

> library(OrganismDbi)

# instructions for hg19, for hg38 change accordingly
> gd <- list(join1 = c(GO.db="GOID", org.Hs.eg.db="GO"),
             join2 = c(org.Hs.eg.db = "ENTREZID",
                       TxDb.Hsapiens.UCSC.hg19.knownGene = "GENEID"))
> destination <- tempfile()
> dir.create(destination)
> makeOrganismPackage(pkgname = "Homo.sapiens.hg19", graphData = gd,
                      organism = "Homo sapiens", version = "1.0.0",
                      maintainer = "Maintainer<maintainer@email>",
                      author = "Author Name", destDir = destination,
                      license = "Artistic-2.0")

Creating package in /var/folders/sf/tfrswpc96jl8b4520mn784440000gn/T//RtmpzR4YnY/file111b490af70f/Homo.sapiens.hg19

# installing package
> install.packages(destination, repos = NULL, type="source")

```

Then, you can load and use the packages with:

```
> library(Homo.sapiens.hg19)
> library(Homo.sapiens.hg38)

> select(Homo.sapiens.hg19, keys=c('KRAS'),
     cols = c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'),
     keytype = 'SYMBOL')
'select()' returned 1:many mapping between keys and columns
  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25358180 25403854
2   KRAS   chr12 25386768 25403863
> select(Homo.sapiens.hg38, keys = c('KRAS'), cols =
       c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'), keytype = 'SYMBOL')
'select()' returned 1:many mapping between keys and columns
  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25204789 25250931
2   KRAS   chr12 25209168 25250936
3   KRAS   chr12 25209431 25250803
4   KRAS   chr12 25233819 25250929

```
> 3 votes

---
Tags: gene, reference-genome
---


# Question
Title: How to verify if the stationary states of a biological model are correct?
I'm mathematically modeling yellow fever, and yet, after assembling the model, I need to find the stationary states of the model. I had to do it by hand because wolfram software could not get the necessary values.

My question is, how do I check if my stationary states are in line with my biological model?

I was told that there is a way to verify this by performing a substitution of the s values of the points within my own template code. But I have no idea how to do it.

The code used is taking forever to get back to me as this data behaves in my model

```
function ypsir = ypsirtry(t,y)
betah = .01;
muh = 0.01;
sigmah = 0.05;
alphah = 0.06;
deltah = 0.05;
mui = 0.07;
mum = 0.07;
betai = 0.03;
lambdai = 0.06;
lambdah = 0.5;
k = 500;

ypsir = zeros(5,1);
ypsir(1) = sigmah*y(3) - muh*y(1) - betah*y(1)*y(5) + lambdah*(y(2)+ y(3) +
y(1));
ypsir(2) = -deltah*y(2) - y(2)*(alphah + muh) + betah*y(5)*y(1);
ypsir(3) = deltah*y(2) - sigmah*y(3);
ypsir(4) = -mui*y(4) - mum*y(4) - betai*y(2)*y(4) + lambdai*(y(4) + y(5))*(1
- (y(4) + y(5))/k);
ypsir(5) = betai*y(2)*y(4) - mui*y(5) - mum*y(5);
end

```

Along with this calling program

```
clear;
to = 0;
tf =50;
yo = [99 1 0 1 3];
[t y] = ode45(@(t,y)ypsirtry(t,y),[to tf],yo);
plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),t,y(:,5))
title('First Model of Yellow Fever')
xlabel('time')
ylabel('susceptible, infected, recovered,Mosquito Non-carrier, Carrier
mosquito')

```

My stationary states are:

```
S = (1/betai)*[(muh+alphah)+(sigmah/deltah)*deltah]*(mui+mum)/betai*N*I);

I = sigmah*R/deltah;

R = ((muh/betah + lambdah)*(muh+alphah))/(sigmah+lamdah)-(muh+alphah)*
(deltah/sigmah*deltah)-(betah/muh - lambdah)*(betai*N*I)*
(delta/sigmah*deltah));

N = -k/lambdai((-lamdai+(lambdai/k))*P/-betai*I + lambdai - (2k/lambdai*P)-
mui-mum);

P = betai*N*I/(mui+mum);

```

At where S = susceptible I = infected R = retweeted N = non-carrier mosquito P = carrier mosquito

# Answer
> 1 votes
After much analysis of the model, I came to the conclusion that I will not be able to verify how the equilibrium points behave graphically. Since you will have to create a phase diagram in R ^ 5, there are five equations, and I can not do two-to-two comparison in these cases, because I will always have a function that depends on a variable that will not be placed, but will be in the equation plotted anyway.

But in any case, if it were in three dimensions, for the purpose of information for those who are interested, the verification of the points would be to place the points (which in calculating have both numbers and not functions in function of other functions), within the model created, hence plot the phase diagram and check whether the arrows converge or not to the point. If they are converging, they are correct.

Regarding the model is taking forever, I took the logistics part and ran well!

---
Tags: modelling, matlab
---


# Question
Title: How can I easily get the read size distribution of reads mapping on a certain set of regions?
Suppose I have a BAM file indicating where reads in a library have mapped, and a bed file describing a set of genomic regions.

Is there a way to easily get the size distribution of the reads mapping on this set of regions?

# Answer
> 4 votes
This can be done rather simply by combining this samtools based answer and bioawk:

```
samtools view -L regions.bed aligned_reads.bam \
    | bioawk -c sam '{hist[length($seq)]++} END {for (l in hist) print l, hist[l]}' \
    | sort -n -k1

```

The `samtools` command extracts the alignements in the region given with the `-L` option, in sam format. The `-c sam` option of `bioawk` makes the `$seq` variable (corresponding to the read sequence field) available in the awk command.

The final `sort` ensures the size classes are ordered numerically.

We obtain a two-column `size count` table.

Examples with timings:

```
$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | bioawk -c sam '{histo[length($seq)]++} END {for (l in histo) print l, histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.592s
user    0m0.816s
sys     0m0.036s

```

Since sam format is actually simple to parse, one can try to use `mawk` instead for slightly higher speed (the sequence being in the 10th field):

```
samtools view -L regions.bed aligned_reads.bam \
    | mawk '{hist[length($10)]++} END {for (l in hist) print l"\t"hist[l]}' \
    | sort -n -k1

```

(Contrary to bioawk, the default output field separator is not a tabulation.)

Example:

```
$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | mawk '{histo[length($10)]++} END {for (l in histo) print l"\t"histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.362s
user    0m0.500s
sys     0m0.004s

```

# Answer
> 1 votes
Via BEDOPS `bedmap` and `bam2bed`, along with the `hist` functionality in `bashplotlib`

```
$ bedmap --echo-map-size --multidelim '\n' intervals.bed <(bam2bed < reads.bam) | hist

 77|           o
 73|           o
 69|           o
 65|           o
 61|           o
 57|           o
 53|           o
 49|           o
 45|           o
 41|           o
 37|           o
 33|           o
 29|           o
 25|           o
 21|           o
 17|          oo
 13|         ooo
  9|         ooo
  5|         ooo
  1| o    o oooo
    -----------

------------------------------
|          Summary           |
------------------------------
|     observations: 118      |
|    min value: 16.000000    |
|      mean : 34.415254      |
|    max value: 36.000000    |
------------------------------

```

---
Tags: bed, sam, text-processing
---