# Question
Title: How can I use Arlequin via the command line?
I've got a decent knowledge of programming (incl. bash scripting) but I fail to understand how Arlequin works. Could you please help me with a very simple reproducible example on how to use Arlequin via the command line?

As the .zip file comes with its pre-compiled version, I'll just note that I'm on Mac OS X but of course, a reproducible example with Linux (or at worst with Windows) would also be helpful.

# Answer
> 3 votes
For MAC OSX, the executable found in `arlsumstat_macosx` (`arlsumstatmac_64bit`) seem to work appropriately. The directory for `arlsumstat_macosx` does not contain any examples but you can use example files from `arlecore_macosx/Example files_linux` using this syntax.

The complication is that `arlsumstatmac_64bit` will look in the current directory for the file `arl_run.ars`. There is a file with the same name in `arlecore_macosx` but it is not compatible. So you have to be in the `arlsumstat_macosx`.

In short, `cd` to `arlsumstat_macosx` and do

```
./arlsumstatmac_64bit path/to/arlecore_macosx/Example\ files_linux/Freqncy/cohen.ars path/to/arlecore_macosx/Example\ files_linux/Freqncy/cohen.arp

```

, where the `.arp` file contains the data and the `.ars` file contains the settings (the description of what you want Arlequin to do).

# Answer
> 1 votes
A concise description for using Arlecore on terminal for MAC, LINUX or Windows is provided here: http://cmpg.unibe.ch/software/arlequin35/man/arlecore\_readme.txt

---
Tags: genome, statistics, positive-selection
---


# Question
Title: What is a signal 'chunk' in the context of Nanopore sequencing?
Does one chunk correspond to the raw signal of a single read read? Why does messing with this parameter affect basecalling results?

# Answer
If you have access to the nanopore community website you can read more about this here: https://community.nanoporetech.com/posts/not-chunky-enough (I'm not sure to which extent I can copy paste from there so I won't).

A TLDR version: basecalling doesn't happen simultaneously for an entire read, but it is broken in segments which are processed separately and combined afterward. Those segments are chunks. Basecalling results are affected mostly by low-complexity regions, in which a heavily biased region (in terms of nucleotide content) leads to artifacts in the per-chunk basecalling. Those issues are solved by taking a chunking-size larger than most low-complexity regions.
> 1 votes

---
Tags: nanopore
---


# Question
Title: Short reads vs long reads with regards to break points
My professor in class said the following statements.

1. Short reads can map break points
2. Long reads can disambiguate repeat related issues

I'm very new to genomics, so I don't quite understand the first statement. What is meant by short reads can map break points. What are break points to begin with?

# Answer
> Short reads can map break points

With short reads you do an "indirect observation/inference" of a structural variant: based on your paired-end reads you figure out something is wrong because their distance or orientation is not what you expect. With long reads, you often get a direct observation of a breakpoint, with sufficient long reads on both sides (either split alignment or spanning the SV) allowing almost nucleotide accuracy.

> Long reads can disambiguate repeat related issues

Short reads cannot align unambiguously in many repetitive regions, because the repeat is longer than the read. Long reads, however, span most common repeats entirely, allowing you to unambiguously assign a location for each read, as such resolving the issues.
> 3 votes

# Answer
A "break point" typically denotes one end of a structural variation, such as an inversion or deletion. You can simply think of it as the position where the normal chromosome was broken to allow the variation. One can find these with short (Illumina) reads relatively easily since with paired end sequencing you'll have pairs spanning these breakpoints with odd relative orientations or insert sizes or such. Note that this is not a particular strength of short reads, since long reads are probably even better at this (they are not as good for finding SNPs at the moment though).
> 3 votes

---
Tags: ngs
---


# Question
Title: removing mitochondrial read and unassembled "random" from multiple bam files
This is the one i used for a single bam file to filter its mitochondrial as well as unassembled read

```
samtools view -h SRR2920466.bam | awk '{if($3!= "chrM" && $3 != "chrUn"){print $0}}' | samtools view -Shb - > this.filter.bam

```

Is there any way i can do it on multiple files with parallel as i do to index and sort bam files using samtools .

So i tried this ,which is of course not working as its not able to pipe properly because it shows me error regarding additional fonts

```
ls *.bam | parallel "samtools view -h {} | awk '{if($3 != "chrM" && $3 != "chrUn"){print $0}}'| samtools view -Shb {.}.bam"

```

Any help with my above codes? As i get it there is wrong with my pipe .

After updating the samtools to version 1.9 i ran the small script

i get this error

```
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string

```

I have 30 bam files i tried with 2 to see if its working or not , after the suggestion of @WouterDeCoster i tried the above one after updating samtools . the basic samtool with egrep works but for multiple file i couldn;t do it.

# Answer
Your awk command is inside double quotes (`""`) which means that all of the awk field variables (e.g. `$3`) are being interpreted by the shell. Since `$3` is not set in the shell, your awk is trying to compare an empty string and that's why you get the error:

```
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error

```

Note how the error message shows `if( != chrM` and not `if($3 != chrM`. Since you *need* to use double quotes here because you need both an outer (for `parallel`) and an inner (for awk) quote, and single quotes can't be nested, you will have to escape the awk variables:

```
ls *.bam | parallel "samtools view -h {} |
    awk '{if(\$3 != "chrM" && \$3 != "chrUn"){print \$0}}' | samtools view -Shb {.}.bam"

```
> 1 votes

---
Tags: awk, shell
---


# Question
Title: Converting a VCF into a FASTA given a reference with Python, R
I am interested in converting a VCF file into a FASTA file given a reference sequence with Python or R.

Samtools/BCFtools (Heng Li) provides a Perl script `vcfutils.pl` which does this, the function `vcf2fq` (lines 469-528)

This script has been modified by others to convert InDels as well, e.g. this by David Eccles

```
./vcf2fq.pl -f <input.fasta> <all-site.vcf> > <output.fastq>

```

Given that converting a VCF into a FASTA given some reference is a somewhat standard procedure in bioinformatics, are there standard R or Python implementations used by the community?

# Answer
You can whip up something quite easily in Python using pyfaidx, which allows you to access a FASTA file by genomic coordinates. So you'd just have to pull the coordinates from the VCF to get the sequence, then output in FASTA format with whatever unique identifier for each sequence you want to use.

Edit: The `FastaVariant` class is new in `pyfaidx`, I hadn't seen it before. It makes this even easier, as you won't have to manually mess with the VCF to get the positions. Something like the following should work:

```
consensus = FastaVariant('yourgenome.fasta', 'variants.vcf.gz', het=True, hom=True)

out = open("variants.fasta", "w")

for chrom in consensus.keys:
    for var in consensus[chrom].variants_sites:
        record = consensus[chrom][var-1:var]
        print(record, file=out)

out.close()

```
> 5 votes

# Answer
I tried the `pyfaidx` approach reported by Jared, but I got unhandled exceptions that I could not handle easily within the framework. I had `bcftools` lying around, so I tried `bcftools consensus` and it worked like a charm. However, you do need an indexed VCF.

```
bcftools consensus all-site.vcf.gz < input.fasta > output.fasta

```

NB: `bcftools consensus` has a few options specified with the `--haplotype` argument for choosing which alleles should be incorporated into the FASTA file.
> 3 votes

# Answer
Assuming you have `samtools`-indexed FASTA files for your reference genome, and have BEDOPS and htslib installed, you could use bed2faidxsta.pl to convert BED to FASTA:

```
$ vcf2bed < variants.vcf | bed2faidxsta.pl --fastaDir=/path/to/fasta > variants.fa

```

To index your FASTA files, you could do the following:

```
$ cd /path/to/fasta
$ bgzip -c my_genome.fa > my_genome.fa.bgz
$ samtools faidx my_genome.fa.bgz

```
> 2 votes

---
Tags: r, python, fasta, vcf, samtools
---


title: PLINK - Transposed BED file?

text:
# Question
Title: PLINK - Transposed BED file?
I'm currently working with the PLINK file system to store large amounts of genotype data.

The plain format consists of three files, two files for storing phenotype and marker information as well as one file (.PED) containing the genotype matrix.

There exists a binary version of the .PED file (.BED), which stores the data in a number of consecutive byte blocks, with each block encoding the genotype of all samples for a specific marker.

While there exists a transposed version of the .PED file (.TPED) which is also supported by PLINK, I couldn't find out whether there also exists a transposed version of the binary .BED file.

Does it exist at all? Are there other file formats to compactly store a genotype matrix for fast sample / marker extraction?

# Answer
Yes, this exists, and can be efficiently generated by plink 2.0's `--export ind-major-bed` command. (The third byte is 0 instead of 1, and the specification is otherwise identical to that of regular plink .bed files with samples and variants swapped.)
> 2 votes

---
Tags: gwas, plink, genotyping
---


# Question
Title: Interpretation of hierarchical clustering with bootstrapping
I have data that includes 'cases' and 'controls' and have carried out hierarchical clustering.

I used the `pvclust` package in R to bootstrap the results and significant branches are highlighted with red rectangles (based on au\>0.95):

What is clear is that no clustering occurs that separates 'cases' and 'controls', and this is in fact what we expected and want to show. We want to show that the measured variable does not distinguish between cases and controls.

1. List `itemApart` from saying visually no clear clusters emerge that distinguishes between cases and controls are there any objective measures that can be used to say no significant clustering occurs between two groups?
2. One observation I have here is that the AU values and the BP values are very different, even though both p-values should be interpreted in a similar fashion, am I missing something?
3. **Perhaps `pvclust` (bootstrapping) is not the right option here, is there a better way of showing quantifiably that no significant clustering occurs between two groups? Perhaps some kind of supervised clustering (I am not sure what this even means in this context)?**

# Answer
> 1 votes
Good question. You have performed the first step via unsupervised learning and AU bootstraps give great results, orthodox bootstraps give no significant (or very few) clusters. The contrast between these approaches is unusual because they should be somewhere close to each other. If you can resolve that, your analysis is fine. So at present IMO the unsupervised learning needs further investigation. ... deleted

... **antibodies** ... oh I see. Trees, like you've done, don't work here because of the cross-reactivity creates network phylogenies. This is the antonym of bifurcating trees. If someone has a robust bifurcating phylogeny with antibody data and gets good bootstraps .. utterly amazing. However, if you subject your data to a network analysis you'll find connections all over the place - in practical terms this really messes up bootstraps. Thats the likely reason for you results. There is no point doing a complex bifurcating tree analysis because the antibody cross-reactivity will mess it up. k-means has been implemented as unsupervised learning in these instances and would provide the starting point in using clustering to define a supervised learning problem.

# Answer
> 1 votes
I've got it. The question you are asking is homogeneity between you target group and your general population of patients. What you want to show is "complete" homogeneity between the two groups regardless of any sub-sampling, ie. if there is groups within groups that might break your hypothesis. K-means is perfect here. The issue with K-means is designating the number of groups, but here thats irrelevant because its not the question you are asking. You are asking homogeneity and if your hypothesis is true no matter how many groups you designate for K-means, e.g. 2,3,4,5,6,7,8,9,10 (or maybe a more reduced sampling) homogeneity between target (treatment?) and general population is observed. This would therefore be a definitive analysing ASSUMING that you don't encounter a group where there is exclusively one population and the other is absent. Obviously the number of K-mean groups can't exceed the number of treated/target patients.

It would be a really good application of k-means because it avoids the question, So why was that a priori selected?

Personally I would first present your preferred a priori group size and then say same result for all group sizes. Its simple, recognised and neat and easy for a reviewer to follow. There are ML solutions to K-means, but I think "keep it simple".

---
Tags: r, clustering
---


title: Knowing which positions are variants

text:
# Question
Title: Knowing which positions are variants
Let's say I'm sequencing a donor, and I have 20X coverage. And let's say all the reads covering a particular base has the same allele at the position. Then, is that position considered a "variant" with only 1 allele in that position or is that not even considered a variant?

The reason I'm asking is, I was given some read data, and it's in the form location:variant - location:variant.

But there was often the case where for all the variants at a particular location, there was only one allele. But this data is labelled to only have variants. Hence, how is a particular location deemed a "variant."

For example, lets say I'm told that the following is **"variant reads given in the form read# location:variant - location:variant"**

```
read 1 6523432: A - 6710432 T
read 2 6523432: C - 6710432 G
read 3 6523432: A - 6710432 T
read 4 6523432: C - 6710432 G

```

Then, why are these two locations considered variants? It seems like it's just a heterozygous locations. Why were these locations considered to carry variants? Similarly,

```
read 1 6523432: A - 6710432 T
read 2 6523432: A - 6710432 T
read 3 6523432: A - 6710432 T
read 4 6523432: A - 6710432 T

```

In those reads, it just seems like it's homozygous. Why is it being called a variant? How was it classified as such?

Put another way, if I'm given 100000 reads, and I said I only want to look at **informative reads**, and so I'm given a subset with only 100 reads, and those 100 reads look something like the above, then how were these reads selected?

# Answer
1. Set 1

```
read 1 6523432: A - 6710432 T
read 2 6523432: C - 6710432 G
read 3 6523432: A - 6710432 T
read 4 6523432: C - 6710432 G

```

I understand this as showing the location of 2 variants, supported by 2 reads each. Now, if the reference has an A at position 6523432, your sample could be heterozygous for an A-\>C variant at position 6523432. If the reference has a C at that position, then your sample could be heterozygous for C-\>A. If the reference has something else again, for example G, then your sample could be heterozygous for G-\>A/C. Similarly for location 6710432.
2. Set 2

```
read 1 6523432: A - 6710432 T
read 2 6523432: A - 6710432 T
read 3 6523432: A - 6710432 T
read 4 6523432: A - 6710432 T

```

In this case, if the reference has an A at position 6523432, the reads do not support a variant. However, if the reference has anything else, then the reads would support a homozygous variant e.g. G-\>A. And the same for position 6710432.

Now, in all of these cases, you cannot know whether a variant is real without taking into account the number of reads that support it, the quality of the alignment, and the quality of the reads. That's why we do this using dedicated variant callers, software that is designed to take all these factors into account and decide whether to call a variant at that position.

But in all cases, the variant must be defined with respect to a reference genome. You cannot call a variant just by looking at reads without any knowledge of a reference. A "variant" just means "not what is found in the reference". It is largely a matter of random chance whether something is a "variant" or "normal". If the individual(s) whose genome was taken as a reference has the same sequence, then you call it normal. If not, then you call it a variant. There are many cases where the standard hg19 or hg38 genomes have a "normal" residue which is actually far less common than the known "variant". There's nothing special about the reference genome, after all, it's just the one we chose to call a reference.
> 3 votes

---
Tags: variant-calling, terminology
---


title: Proteins with one SS bond?

text:
# Question
Title: Proteins with one SS bond?
I would like to find proteins with exactly one SS bond. Is there a database where I can search this? I've tried advanced search on https://www.rcsb.org/, but no such option, at least I could not find it.

# Answer
In cases when you cannot find answer through the web interface, you can download all the data and search it locally.

Relying on the PDB annotations of SS bonds, you could search for PDB files that have exactly one SSBOND record:

```
$ cd $PDB_DIR/structures/divided/pdb/
$ zgrep -c ^SSBOND aa/pdb*.ent.gz | grep :1$
aa/pdb2aaq.ent.gz:1
aa/pdb3aad.ent.gz:1
aa/pdb5aap.ent.gz:1
aa/pdb6aa4.ent.gz:1
aa/pdb6aa5.ent.gz:1

```

All PDB files (gzipped) take about 30GB of disk space.

The caveat is that some entries are available only in the mmCIF format. Using mmCIF would be slighty more complex: you would need to count `_struct_conn.conn_type_id` being `disulf`. And you would need to download 10GB more.
> 3 votes

---
Tags: public-databases, protein-structure
---


# Question
Title: Genotyping technologies do not maintain phase?
Why exactly does genotyping not maintain phase? I guess I'm more confused on what's the difference between sequencing and genotyping, and as a result how phase is maintained in one, but not the other. This is an image trying to do phasing using genotyping methods.

I understand how phasing works, I just don't quite get what is the difference between genotyping and sequencing, and why one of them doesn't work.

# Answer
Phasing means that you are able to distinguish the both copies of your target sequence. (Usually you have one from your mother and one from your father).

If you can do this, depends mainly on the sequencing technology you are using. For example in case of sanger sequencing, you cannot do this, because you (pre)amplify all your copies at once and detect the sequence at the same time at the same location.

In case of ngs you have read information at the end. One read represents the sequence of one original dna molecule. If you detect multiple variants in one read, you can be sure that they are from the same chromosome and thus they are in phase. Now you can go on and try to find other reads that have unambiguous overlaps with these variants on your read. If you find some, you can also say that these variants have their origin on the same chromosome.
> 2 votes

# Answer
Many (small) variant calling algorithms use pileups to genotype variants position-by-position. This approach alone is insufficient to phase adjacent variants. However, if the variants of interest are close enough to be covered by a single read, it's certainly possible to phase them even with \<100bp Illumina reads. In fact, I was implementing just such an approach recently when I asked this question about accessing bases by *aligned* position rather than by read position.
> 1 votes

# Answer
Genotyping (there are of course many different methods) is likely to give you information in the format of

Locus 1: 100% A

Locus 2: 50% C, 50% T

Locus 3: 50% C, 50% A

That's not phased.

Short or long read technology will give you data in the form of reads, each of which originates from one molecule, and if a read covers two loci, you'll know how they are phased. You will know if one DNA molecule has ACC, or ACA.

Sanger sequencing will give you mixed peaks at the locus of a mixed SNP, so there will be no phasing data.
> 0 votes

---
Tags: ngs, haplotypes, genotyping, phasing
---


# Question
Title: Why do some assemblers require an odd-length kmer for the construction of de Bruijn graphs?
Why do some assemblers like SOAPdenovo2 or Velvet require an odd-length *k*-mer size for the construction of de Bruijn graph, while some other assemblers like ABySS are fine with even-length *k*-mers?

# Answer
> 31 votes
From the manual of Velvet:

> it must be an odd number, to avoid palindromes. If you put in an even number, Velvet will just decrement it and proceed.

the palindromes in biology are defined as reverse complementary sequences. The problem of palindromes is explained in this review:

> Palindromes induce paths that fold back on themselves. At least one assembler avoids these elegantly; Velvet requires K, the length of a K-mer, to be odd. An odd-size K-mer cannot match its reverse complement.

It is possible to construct graph with palindromes, but then the interpretation will be harder. Allowing only graphs of odd *k*-mers is just an elegant way to avoid writing a code for interpretation of a more complicated graph.

# Answer
> 14 votes
To expand on the answer above, in case it isn't clear, we show:

1. Why palindromic sequences must be even in length
2. Why palindromic sequences induce self-loops in a de Bruijn graph
3. Why self loops in a de Bruijn graph are problematic

---

## 1. Palindromic sequence ⇒ sequence is of even length

> **Idea:** in an odd-length k-mer, its middle nucleotide is 'flipped' in its reverse complement, so the two can never be equal.

Suppose you have a palindromic sequence $X$. Then $X$ is identical to its reverse complement, which we will label $\bar{X}$.

Suppose $X$ is of odd length. Then it is of the form $AbC$, where $len(A) = len(C) = \frac{len(X)-1}{2}$, and $len(b) = 1$.

Then

$X = \bar{X} \implies AbC = \overline{AbC} = \bar{C}\bar{b}\bar{A}$

And hence:

$b = \bar{b}$

(since $len(A) = len(C) = len(\bar{C}))$. But this is a contradiction, since $b$ is a single nucleotide, and cannot be equal to its complement. Therefore k-mers of odd length cannot form palindromes.

Hence the length of a k-mer that forms a palindrome must be even.

---

## 2. Why palindromic k-mers induce self loops

Each node in a traditional de Bruijn graph is a unique string, but in most bio-informatics implementations each pair of reverse-complementary k-1-mers are identified as a single node, e.g for $k=6$:

A palindromic k-mer (of $k \geq 2$) is of the form:

$xAy$

where $len(A) = k-2$, $x=\bar{y}$ and $A=\bar{A}$ (possibly the empty string).

Hence it will contribute two nodes in the de Bruijn graph:

1. its left k-1-mer $xA$
2. its right k-1-mer $Ay$

And an edge going from 1 to 2.

But since this k-mer is palindromic, $xA = \overline{Ay}$ and hence these two nodes are reverse-complementary, and thus the 'same' node, and so this edge is a self loop at this node.

---

## 3. Why are self loops problematic?

Self-loops (if they occur in a node with $in\\\_degree \geq 2$ and $out\\\_degree \geq 1$) increase the number of possible Eulerian paths in a de Bruijn graph (or more specifically, in the connected component containing this node, which represents a *contig*, of which their may be multiple), as you have an additional possible Eulerian path for each time you traverse this node.

This increases ambiguity in reading the graph, as each possible Eulerian path is an extra possible reconstruction of the full sequence.

Consider the example:

There is only one possible Eulerian path:

* $A B C D B E$

However, if we include a self loop at $B$, which is visited twice above, this doubles to two possible Eulerian paths:

* $A B B C D B E$
* $A B C D B B E$

Depending on whether we traverse the self loop during the first time we reach $B$, or the second.

---
Tags: ngs, assembly, k-mer, de-bruijn-graphs
---


# Question
Title: Find similar proteins with similar tertiary structure?
I'm a mathematician studying knotted curves and am trying to find some applications of my work in proteins. Using https://www.rcsb.org/, how do I find proteins that have similar (but perhaps not the same) tertiary structure? Preferably they would have similar primary structure and function... Perhaps what I'm looking for is "similar proteins", if there is such a thing?

Also, is there a way in the database to select more "common proteins" (important or well studied, not just some rare exotic one (if there is such a thing)?

# Answer
As a general question on the board, more detail is often good, because of the answers can be inaccurate due to misunderstanding of the biological context.

There are several ways to do this and at least 3 that I'm aware of.

1. Blast - sequence similarity searching the PDB database
2. Contact map identity, I use this approach, but its complicated - and its precision or utility is context specific, i.e. it could be loads of work and miss the vital detail you need.
3. Direct mapping one protein structure against another, this in particular is true for specific domains, or sites within the protein and the result is angstroms. Its not difficult and I assume this can be scaled to the whole protein.

**Point 1 "Blackboxes"**

By far the easiest is option one and I suspect I-Tasser (here) use this principle before instigating "threading". This package is considered "good" (others appear to consider it the best of its type), I believe it has done well in CASP tests (Critical Assessment of Techniques for Protein Structure Prediction). So the simple answer is dump your sequence in there an watch the output. Pyre is something similar. These are interesting tools but certainly not definitive, really there a starting point.

**Point 1 "Formal Blasting PDB"**

This is the easiest and most mature approach for your question, without knowing the details. I've run a test to ask the following question.

```
Which protein structure available on PBD is most similar to Japanese encephalitis virus (JEV) based on protein homology?

```

**Step A** I've no specific interest in the virus BTW. The most robust Blast engine IMO is NCBI here . I dump a JEV sequence into the box, and click three important options (otherwise it doesn't work).

* Database = protein databank (PDB), i.e. structural data
* Organism "*Japanese encephalitis virus (taxid:11072)*": the taxid will come up automatically (if not just proceed as normal)
* Next to the Organism box I click "*Exclude*", otherwise this whole approach doesn't work
* Obviously I select "*blastp*", otherwise it definitely will not work

The output is/was here (but is now inactive). Summarising what it presented, E is so high its "0" (\>2e-139) protein ID "77.59%" and just states "envelope protein" Chain A. I hit the link and the protein that comes up has a name which ends in "**6A0P\_A**". This is the PDB reference and Percentage ID and PDB reference are the only two bits of info you need. The PDB reference needs the "\_A" removing, thus its 6A0P. The \_A just refers to chain A of a dimer.

**Step B** You then go back to the Web site you mentioned and input the PDB ID you stated into the search bar. This results in this quite interesting page here. The structure is Usutu virus, which surprised me.

**Point 2** . This is much harder and not only involves specific Python code, but a post unsupervised learning analysis and results in an average similarity of structure over the whole protein. There are a number of published algorithms for generating a contact map, or its easy to write your own. Its more advanced, its complicated and unless you've seen loads of them they not immediately understandable. Contact maps angle at a very specialist approach of predicting and understanding protein structure. It is not for the faint hearted and some biological details will be omitted, which may or may not be important for your purpose.

**Point 3** . I don't have research grade expertise here, so I forget the exact button to click, but its quite easy via PyMol and other graphical viewers to perform spatial similarity for specific epitopes and this is something you need to learn within your sphere. Basically you will load two proteins into the image viewer and then they will be a "similarity aligning" function and you zoom in on your output to generate an answer in angstroms. You will need to learn this approach for the later parts of your project.

**Summary**, The orthodox way would be to perform point 1, select several targets and then examine them in detail via point 3. I definitely recommend direct Blast if you go down this route, because it is the most tractable and avoids potential unknown subsampling. I personally would not recommend *Point 2* in a general context because you have to consider an unsupervised learning approach towards processing a heap of contact maps. Ultimately, which ever approach you take you can never avoid Point 3 (which someone else does, in my case).

**Note** This may or may not be obvious to your but just because two sequence have the closest homology doesn't mean they will have the closest structural similarity.
> 1 votes

---
Tags: database, protein-structure
---


# Question
Title: Testing for viral/bacterial infection FASTQ files
What's the best way to test if a pair of Illumina FASTQ files from blood DNA extraction contain any reads are from viral or bacterial infection from human samples? Thanks in advance.

EDIT:

Looking at the available options, I started testing Kraken, the metagenomics classification (Wood and Salzberg Gen Biol 2014) for samples where I have seen an indication of viral reads. I will update with results / further questions.

# Answer
If you're talking RNA, I have a blog post from a while back about Diamond/MEGAN:

http://achri.blogspot.com/2017/04/data-for-free-viral-metagenomics-in.html

In principle this will work for DNA if you have a non-trivial number of pathogen reads (since they tend to be gene dense in their genomes). Otherwise I use Kraken on a database for all the NCBI nt database, but this is slower and doesn't let you sanity check the false positive taxonomy assignments as well as I describe in the post. Of course, ideally, you run multiple methods if you have the resources to do so.
> 2 votes

# Answer
I think you could try pathoscope for this.
> 2 votes

# Answer
To know if there are bacterial genomes in the sequenced sample I have used an approach based on kmers analysis. An useful tool to do that is kmerfinder which uses a kmers database built from bacterial genomes. Quoting from the kmerfinder web site:

> The method is used to find the best match (species identification) to the reads in one or more fastq files or one fasta file in a (kmer) database produced using the KMA program.

I haven't use it for viruses genome analysis, but conceptually it should be work as well. Hope this would be helpful.
> 1 votes

# Answer
You may find this blog post helpful: Identifying viral sequences in TCGA data using Kraken and Centrifuge.

> Our goal was to evaluate the performance of Centrifuge and Kraken for detecting viral transcripts across a cancer cohort. Centrifuge and Kraken produce very similar results and consistently detect HPV strains in cervical cancer.

There is also a review Overview of Virus Metagenomic Classification Methods and Their Biological Applications, but they focus specifically on viruses:

> We evaluated 49 published computational classification workflows for virus metagenomics in a literature review. To this end, we described the methods of existing workflows by breaking them up into five general steps and assessed their ease-of-use and validation experiments. Performance scores of previous benchmarks were summarized and correlations between methods and performance were investigated. We indicate the potential suitability of the different workflows for (1) time-constrained diagnostics, (2) surveillance and outbreak source tracing, (3) detection of remote homologies (discovery), and (4) biodiversity studies.
> 1 votes

# Answer
What type of NGS data is it? If smallRNA or RNAseq data, you can simply assemble the reads into longer contigs and do the blast search using your reference viral or bacterial sequences.
> 0 votes

---
Tags: fastq, ngs, covid-19
---


# Question
Title: Control width of samtools tview "snapshot" when redirecting
I have a large number of loci I would like to examine manually with `samtools tview`. Rather than typing or copy-n-pasting dozens of coordinates, I was hoping to automate the generation of pileup "images" (text files).

With some experimenting just now, I was happy to learn that if I specify the display option `-d H` or `-d T` the output can easily be redirected to a file for later inspection. However, it appears that I cannot control the window captured. Changing the width of the terminal window doesn't change the output, nor does adjusting the `-p` flag (i.e. `-p chr1:501-700`). No matter what, `tview` always reports an 80bp window in text or HTML mode.

Is there any way I can adjust this?

# Answer
I cross-posted this question on the official samtools issue tracker, and finally got a response.

> Looking at the code for the HTML view, I see this and this. So maybe try
>
> `COLUMNS=200 samtools tview -d H -p CHR:X-Y aln.bam ref.fasta`
> 1 votes

# Answer
Partially a shameless plug. See if ASCIIGenome can help you with that.

To review a number of regions that you saved in, say, a bed file, you can either process them in batch or load BAM(s) and bed file and jump from one bed interval to the next using the `next -start` command.
> 2 votes

---
Tags: samtools, read-mapping, visualization
---


# Question
Title: What is the best argument for maximum parsimony method in phylogenetic tree construction?
Beside it (maximum parsimony) being compuationally cheap what other good arguments are there for it? Is there any model behind this principle ? Why would one expect this principle to provide right phylogeny in any situation at all ?

# Answer
> *The principle of parsimony states that the simplest explanation consistent with a data set should be chosen over more complex explanations...* (Stewart 1993)

Sometimes independent events in different lineages cause similarities in molecular sequences (or morphology or behavior or whatever characteristics one happens to be studying). These events are not true evidence of recently shared ancestry, but it is easy to mistake them as such evidence.

The principle of parsimony asserts that we will be better off assuming the simplest case--that is, that identical characters ***are*** evidence of recent shared ancestry, since the simple explanation should be true more frequently than alternative explanations. There are cases for which the maximum parsimony model works very well, and there are cases for which maximum likelihood, distance matrix, Bayesian, or other models perform better.

---

Stewart CB (1993) The powers and pitfalls of parsimony. *Nature* **361**, 603-607, doi:10.1038/361603a0.
> 2 votes

# Answer
I'd say maximum parsimony approaches assume genetic similarity is unlikely to occur by convergence, but likely a result of divergence. If this assumptions holds, finding the tree configuration that minimizes the number of unlikely events makes sense.
> 2 votes

# Answer
Parsimony was a nice idea and developed the algorithms and tree methods towards its successors, maximum likelihood and Bayesian. The fundamental assumption of mutational independence at each alignment position was founded via parsimony. Just the recursion algorithm involved in reading a tree isn't trivial coding. Also Walter Fitch (Fitch Marigold method) was a nice guy.

Weighted parsimony has performed well in simulation studies to prevent 'long branch attraction' artefacts. The key issue that was never resolved was how to objectively weight an MP analysis. It could be resurrected however, because machine learning ML could independently provide the weights in specific contexts. That a lot of "coulds" though.
> 1 votes

---
Tags: phylogenetics, phylogeny
---


# Question
Title: Evidence of emergence of "genuine" novel protein interactions by aa mutation
Anyone knows of an example of a protein that, without coming from a recent duplication event, underwent a mutation(s) that caused it to have a novel interaction with a new ligand, substrate, other protein or molecule? By novel I imply that it doesn't "go back" to an ancestral interaction or doesn’t develops an interaction that it's already performed by a close relative. I assume that these mutations are extremely rare, and that a duplicate is almost always necessary to have them fixed in a lineage.

# Answer
HA Sialic acid receptor binding in influenza virus determining host specificity, viz. Bird flu vs human flu.
> 1 votes

# Answer
That's easy, red blood cells and the single mutation causing the switch in charge resulting in sickle cell anaemia. There are five different convergent mutations although the concept of convergence here has been challenged. The mutation results in the physical collapse of the protein structure and massively reduces its binding with oxygen. The receptor interactions with the malarial parasite render infection of the red blood cell impossible. Even in the heterozygote state, oxygen binding with the heme group is not massively disrupted, but the ability of the malaria parasite in particular plasmodium falciparum, to invade the red blood cell is again greatly impaired. It is an incredible example of human evolution and demonstrates the selective pressure of holoedemica malaria. There is no gene duplication involved here and again due stress there are potentially 5 convergent events resulting in the same phenotype and the scale of protein collapse from circular red blood cell to sickle cell shape is the result of what would otherwise appear to be a small change in amino acid charge for a single amino acid substitution.
> 0 votes

---
Tags: sequence-annotation, phylogenetics, protein-structure, variation
---


# Question
Title: Authoritative source on human cytogenetic regions?
I am looking for a database that would keep track of human cytogenetic regions and genomic coordinates per genome assembly. I had expected the Genome Reference Consortium to have it, or Ensembl, or NCBI, but I could not find it.

Is there a place with a translation table between human cytogenetic regions and genomic coordinates?

## EDIT

After benn's answer I found that you can also get this information from Ensembl's REST API: https://rest.ensembl.org/info/assembly/homo\_sapiens?content-type=application/json&bands=1.

I have transformed the data provided by Ensembl and generated a CSV file with the cytogenetic regions:

https://github.com/ramiromagno/gwasrapidd/blob/master/data-raw/cytogenetic\_bands.csv

# Answer
> 5 votes
Do you mean cytoband coordinates? For the human genome you can find them for example at UCSC.

```
hgdownload.cse.ucsc.edu/goldenPath/hg38/database/cytoBand.txt.gz

```

---
Tags: assembly, human-genome
---


title: Repeatmasker output table

text:
# Question
Title: Repeatmasker output table
I'm working with an unassembled genome of cephalopod and I would to find transposable elements. My idea is to use repeatscout to generate a library and repeatmasker to find the elements. After creating a library with repeatscout, I run Repeatmasker but it generates a table like this:

```
 file name: OB_100DEC.fa
sequences:           100
total length:  209466439 bp  (184235452 bp excl N/X-runs)
GC level:         35.24 %
bases masked:   73100024 bp ( 34.90 %)
==================================================
               number of      length   percentage
               elements*    occupied  of sequence
--------------------------------------------------
SINEs:                0            0 bp    0.00 %
      ALUs            0            0 bp    0.00 %
      MIRs            0            0 bp    0.00 %

LINEs:                0            0 bp    0.00 %
      LINE1           0            0 bp    0.00 %
      LINE2           0            0 bp    0.00 %
      L3/CR1          0            0 bp    0.00 %

LTR elements:         0            0 bp    0.00 %
      ERVL            0            0 bp    0.00 %
      ERVL-MaLRs      0            0 bp    0.00 %
      ERV_classI      0            0 bp    0.00 %
      ERV_classII     0            0 bp    0.00 %

DNA elements:         0            0 bp    0.00 %
     hAT-Charlie      0            0 bp    0.00 %
     TcMar-Tigger     0            0 bp    0.00 %

Unclassified:    560941     66374119 bp   31.69 %

Total interspersed repeats: 66374119 bp   31.69 %

Small RNA:            0            0 bp    0.00 %

Satellites:           0            0 bp    0.00 %
Simple repeats:  290711     15358198 bp    7.33 %
Low complexity:       0            0 bp    0.00 %
==================================================

* most repeats fragmented by insertions or deletions
  have been counted as one element

The query species was assumed to be homo
RepeatMasker Combined Database: Dfam_Consensus-20181026

run with rmblastn version 2.6.0+
The query was compared to unclassified sequences in ".../OB_100DEC_repeats_filtered1.fa"

```

This was my command to run repeatmasker:

```
./RepeatMasker -s -lib /home/RepeatScout-1.0.5/OB_100DEC_repeats_filtered1.fa /home/Workdirectory/OB_100DEC.fa

```

The file `OB_100DEC_repeats_filtered1.fa` is the output of repeatscout and looks like this:

```
>R=3 (RR=4.  TRF=0.000 NSEG=0.000)
TAAGGCGGCGAGCTGGCAGAATCGTTAGCACGCCGGGCGAAATGCTTAGCGGTATTTCGTCTGTCTTTACGTTCTGAGTT
CAAATTCCGCCGAGGTCGACTTTGCCTTTCATCCTTTCGGGGTCGATAAAATAAGTACCAGTTGAGCACTGGGGTCGATG
TAATCGACTTACCCCCTCCCCCAAAATTTCTGGCCTTGTGCCTATATTAGAAACGATTATT
>R=4 (RR=5.  TRF=0.122 NSEG=0.226)
ACACACACACACACACACACACACACATATATATATATATACATATATACGACGGGCTTCTTTCAGTTTCCGTCTACCAA
ATCCACTCACAAGGCTTTGGTCGGCCCGAGGCTATAGTAGAAGACACTTGCCCAAGGTGCCACGCAGTGGGACTGAACCC
GGAACCATGTGGTTGGTAAGCAAGCTACTTACCACACAGCCACTCCTGCGCCTATATATAT
>R=6 (RR=7.  TRF=0.134 NSEG=0.247)
TTGTTTCAGTCATTTGACTGCGGCCATGCTGGAGCACCGCCTTTAGTCGAGCAAATCGACCCCAGGACTTATTCTTTGTA
AGCCTAGTACTTATTCTATCGGTCTCTTTTGCCGAACCGCTAAGTTACGGGGACGTAAACACACCAGCATCGGTTGTCAA
GCGATGTTGGGGGGACAAACACAGACACACAAACACACACACACACATACATATATATATATATATATATA
>R=8 (RR=9.  TRF=0.391 NSEG=0.000)`

```

and it is composed by all the TE (not classified) find in all the scaffolds of the genome.

Can anyone explain why there are no TE? The fasta file that I used is only a part of the whole genome, it is about 200 mb and it is composed by the 100 greatest scaffolds of my genome. Maybe I made a mistake with the use of libraries. Anyway I installed RepeatMasker with rmblastn version 2.6.0+ and RepeatMasker Combined Database: Dfam\_Consensus-20181026. thank you

# Answer
As far as I can tell, you use the `-lib` option for species where there is no existing RepeatMaster database. Therefore, all it does is attempt to locate versions of the repetitive elements from the `-lib` file in the target genome file. However, since your lib file doesn't have any information on what type of element each of your hits is, it looks like RepeatMasker is treating all of them as "Unclassified".

Note the following from the `repeatmasker.help` file you should have downloaded along with RepeatMasker:

> The recommended format for IDs in a custom library is:
>
> \>repeatname#class/subclass
>
> or simply
>
> \>repeatname#class
>
> In this format, the data will be processed (overlapping repeats are merged etc), alternative output (.ace or .gff) can be created and an overview .tbl file will be created. Classes that will be displayed in the .tbl file are 'SINE', 'LINE', 'LTR', 'DNA', 'Satellite', anything with 'RNA' in it, 'Simple\_repeat', and 'Other' or 'Unknown' (the latter defaults when class is missing).

This makes sense. RepeatMasker cannot know if an arbitrary sequence you give it is a SINE or a LINE or anything else. It will therefore rely on the categories present in the fasta headers. Since your file has no categories, it identifies everything in your lib file as "Unknown", which is presumably what the "Unclassified" section in your output refers to:

```
Unclassified:    560941     66374119 bp   31.69 %

Total interspersed repeats: 66374119 bp   31.69 %

```
> 1 votes

---
Tags: transposable-elements, repeatmasker
---


# Question
Title: NCBI retrieve genome assembly with splicing or "modifications"
I am trying to fetch genome assemblies that include splicing information.
E.g. Here, is a genome assembly of a yeast, where genes show variants like this: `join(5928..7626,7685..7785,7857..7982)`.

I am retrieving the information from NCBI with the following commands:
`efetch -db nuccore -id CP003820.1 -format tabular`. I tried different formats and it changes but I don't get the `join(...)`

1) Am I correct and the `join(...)` indicates "splicing-like" modifications in yeast?
2) How can I get the `join(...)` (any other format is fine too) with efetch?

# Answer
> 3 votes
The page you link to is a GenBank Flat file, the default format used by GenBank. So, to get that format, which includes the `join` lines, use:

```
./efetch -db nuccore -id CP003820.1 -format gb

```

The `join()` lines give the different ways in which a sequence can be built from this gene. For example, one mRNA sequence produced by the locus you mention (CP003820) is:

```
  mRNA            complement(join(100..1497,1556..1767,1823..2274,
                 2322..2782,2846..3126,3206..3890,3958..5263,5322..5422,
                 5494..>5645))

```

The corresponding coding sequence is for the mRNA above is:

```
 CDS             complement(join(168..1497,1556..1767,1823..2274,
                 2322..2782,2846..3126,3206..3890,3958..5263,5322..5422,
                 5494..5645))

```

These aren't necessarily splicing "modifications", and aren't "variants", they are simply the way that the exons are joined to make an mRNA.

Note that since CP003820.1 is "Cryptococcus neoformans var. grubii H99 chromosome 1, complete sequence", there are multiple genes described in there and each will have its own (or multiple, if there is alternative splicing) `join` statements. Specifically, I see 1014 genes described in the file:

```
$ grep -Pc '^\s+gene\s+' CP003820.gb
1014

```

950 mRNAs:

```
$ grep -c 'mRNA.*join(' CP003820.gb
950

```

And 955 proteins:

```
$ grep -Pc 'protein_id=' CP003820.gb
955

```

---
Tags: phylogenetics, data-retrieval
---


title: Extracting the number SNP in each range

text:
# Question
Title: Extracting the number SNP in each range
I have called copy number and I have 2 files (I have shared the link of my files ); One contains some ranges

```
> head(cndata[,c(2,3,4)])
     Chromosome   Start      End
4470          1   51479   817980
4471          1  818499  1136753
4472          1 1138735  2558308
4473          1 2558740  5724264
4474          1 5724940  5749083
4475          1 5749226 12529544
>

```

https://www.dropbox.com/s/wei3i787c9nrz8q/WTSI-OESO\_121\_1pre.copynumber.caveman.txt.csv?dl=0

I have another file like below

```
              Chromosome         Position
rs62635286       1                 51479
rs75454623       1                 114930
rs806731         1                 30923

```

https://www.dropbox.com/s/tiqjhkrpfr9d9gu/WTSI-OESO\_121\_1pre.copynumber.txt?dl=0

In first column I have some SNPs, how I can count the number of SNP in each range? For example based on the second table how many SNP are in range of `51479` to `817980`?

# Answer
> 3 votes
There is like a thousand different ways how to achieve this. You could use a specialised software for this (like bedtools) or calculate it simply in R.

R solution: You can make a function that calculates the number of SNPs in a range and than you can apply it on all the ranges in the table with genomic ranges.

```
snp_table <- read.table('WTSI-OESO_121_1pre.copynumber.txt')
genomic_ranges <- read.table('WTSI-OESO_121_1pre.copynumber.caveman.txt.csv', sep = ',', col.names = c('Chromosome', 'Start', 'End'), stringsAsFactors = F)

how_many_in_range <- function(coords){
    coords = as.vector(coords)
    sum(snp_table$Chromosome == coords[1] & snp_table$Position > as.numeric(coords[2]) & snp_table$Position < as.numeric(coords[3]))
}

genomic_ranges$number_of_snps <- apply(genomic_ranges, 1, how_many_in_range)

```

Now the table `genomic_ranges` has a variable that sais how many SNPs are in the range.

---
Tags: snp, bedtools, wgs, copy-number, text-mining
---


# Question
Title: No refseq transcripts in 2nd half of Y chromosome for HG38
Obviously there is a gap in my understanding of how sex chromosomes are annotated. I've been working on some CNV calls and noticed that most of the 2nd half of ChrY is missing annotations and no reads are being mapped there. Yq12 seems to be all NNNNNNN in the hg38 reference.

Anyone know why?

# Answer
Nevermind, Answered my own question. It's a repetitive region that isn't very well understood

> The human Y chromosome, because it is enriched in repetitive DNA, has been very intractable to genetic and molecular analyses. There is no previous evidence for developmental stage- and testis-specific transcription from the male-specific region of the Y (MSY). Here, we present evidence for the first time for a developmental stage- and testis-specific transcription from MSY distal heterochromatic block. We isolated two novel RNAs, which localize to Yq12 in multiple copies, show testis-specific expression, and lack active X-homologs. Experimental evidence shows that one of the above Yq12 noncoding RNAs (ncRNAs) trans-splices with CDC2L2 mRNA from chromosome 1p36.3 locus to generate a testis-specific chimeric β sv13 isoform. This 67-nt 5′UTR provided by the Yq12 transcript contains within it a Y box protein-binding CCAAT motif, indicating translational regulation of the β sv13 isoform in testis. This is also the first report of trans-splicing between a Y chromosomal and an autosomal transcript.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1832090/
> 2 votes

---
Tags: sequence-annotation, read-mapping
---


# Question
Title: Library for computing BWT-based alignments
I am writing a software tool to which I would like to add the ability to compute alignments using the efficient Burrows-Wheeler Transform (BWT) approach made popular by tools such as BWA and Bowtie. As far as I can tell, though, both of these tools and their derivatives are invoked strictly via a command-line interface. Are there any libraries that implement BWT-based read alignment with a C/C++ API?

Python bindings would also be great, but probably too much to expect.

# Answer
> 11 votes
First, let us remark that there exist several hundred read mappers, most of which have been even published (see, e.g., pages 25-29 of this thesis). Developing a new mapper probably makes sense only as a programming exercise. Whereas developing a quick proof-of-concept read mapper is usually easy, turning it into a real competitor of existing and well-tuned mappers can last for years.

It is not clear from the provided description how long is your reference, how many alignments you need to compute, etc. In certain situations, it may be useful to write a wrapper over existing mappers (e.g., using subprocess.Popen for running the mapper and PySam for parsing the output); while in some other situations, standard dynamic programming may be sufficient (e.g., using SSW with its Python binding).

Let assume that you want to develop a toy read mapper. Most of read mappers are based on a so called seed-and-extend paradigm. Simply speaking, first you detect candidates for alignments, usually as exact matches between a read and the reference (using either a hash table or some full-text index – e.g., BWT-index). Then you would need to compute alignments for these candidates, typically using some algorithm based on dynamic programming, and report the best obtained alignments (e.g., the ones with the highest alignment score).

There exist two big, powerful and well debugged libraries implementing BWT-indexes which can be easily used for building a read mapper:

# Answer
> 9 votes
BWA-MEM can be used as a library. File bwa/example.c shows the basic functionality for single-end mapping. It should give identical mapping to the bwa-mem command line. Header bwa/bwamem.h contains basic documentation. Paired-end mapping is doable, but not well exposed. Several teams, including GATK and MS genomics, have been using bwa-mem as a library.

# Answer
> 4 votes
SeqAn supports BWT tables for use with their parametrisable alignment algorithms.

To use it, follow the general outline for building a SeqAn short-read aligner, and use the `FMIndex` specialisation instead of — as in the example — the `IndexQGram`.

# Answer
> 1 votes
If you're looking for a python API, nanopore have created one (although I am not sure how "maintained" it is) https://github.com/nanoporetech/bwapy

---
Tags: sequence-alignment, bwa, read-mapping, api
---


title: cooler out of memory

text:
# Question
Title: cooler out of memory
I ran into the following error when using `cooler` to load a contact pairs file:

```
ubuntu@ip-172-31-8-42:/data/hg-tmp$ FILE=GSE109344_monkey_fibro_allValidPairs.txt.gz; sudo docker run -v $(pwd):/data $TAG cooler cload pairs -c1 2 -p1 3 -c2 5 -p2 6 /data/chromInfo.rheMac2.
simple.txt:${BINSIZE} /data/$FILE /data/${FILE}.cool;
Traceback (most recent call last):
  File "/usr/local/bin/cooler", line 11, in <module>
    sys.exit(cli())
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 764, in __call__
    return self.main(*args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 717, in main
    rv = self.invoke(ctx)
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 1137, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 1137, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 956, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/usr/local/lib/python2.7/site-packages/click/core.py", line 555, in invoke
    return callback(*args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/cooler/cli/cload.py", line 477, in pairs
    map(pipeline, reader),
  File "cytoolz/functoolz.pyx", line 493, in cytoolz.functoolz.Compose.__call__
  File "/usr/local/lib/python2.7/site-packages/cooler/create/_ingest.py", line 425, in _aggregate_records
    .aggregate(agg)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/generic.py", line 1315, in aggregate
    return super(DataFrameGroupBy, self).aggregate(arg, *args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/generic.py", line 186, in aggregate
    result, how = self._aggregate(arg, _level=_level, *args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/base.py", line 498, in _aggregate
    result = _agg(arg, _agg_1dim)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/base.py", line 449, in _agg
    result[fname] = func(fname, agg_how)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/base.py", line 432, in _agg_1dim
    return colg.aggregate(how, _level=(_level or 0) + 1)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/generic.py", line 760, in aggregate
    return getattr(self, func_or_funcs)(*args, **kwargs)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/groupby.py", line 1227, in size
    result = self.grouper.size()
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/ops.py", line 225, in size
    ids, _, ngroup = self.group_info
  File "pandas/_libs/properties.pyx", line 36, in pandas._libs.properties.CachedProperty.__get__
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/ops.py", line 252, in group_info
    comp_ids, obs_group_ids = self._get_compressed_labels()
  File "/usr/local/lib/python2.7/site-packages/pandas/core/groupby/ops.py", line 272, in _get_compressed_labels
    return compress_group_index(group_index, sort=self.sort)
  File "/usr/local/lib/python2.7/site-packages/pandas/core/sorting.py", line 379, in compress_group_index
    comp_ids, obs_group_ids = table.get_labels_groupby(group_index)
  File "pandas/_libs/hashtable_class_helper.pxi", line 1255, in pandas._libs.hashtable.Int64HashTable.get_labels_groupby
  File "pandas/_libs/hashtable_class_helper.pxi", line 199, in pandas._libs.hashtable.Int64Vector.resize
MemoryError: cannot allocate memory for array

```

The machine has 8Gb of RAM. Any ideas?

# Answer
> 5 votes
It seems to be some kind of a memory leak in pandas's `read_csv`, maybe due to it doing gzip-decompression.

Doing the decompression in a separate process and piping the result into `cooler cload pairs` should do the trick.

```
gzip -dc GSE109344_monkey_fibro_allValidPairs.txt.gz | cooler cload pairs -c1 2 -p1 3 -c2 5 -p2 6 /data/chromInfo.rheMac2.
simple.txt:${BINSIZE} - /data/${FILE}.cool

```

---
Tags: python, errors
---


# Question
Title: State of the art mutation simulation software
There are many features affecting mutation probabilities, e.g. CpG mutations are 10-fold more likely than other types of mutations.

Is there a model (preferably with software) which can take two aligned genomic regions, **estimate** parameters of the neutral mutation procces and then **simulate** mutations for another region? Obviously, there are models like HKY and GTR, but they are somewhat simplistic and do not take the context into account.

One model I found is this, but it seems that there is no software implementing the model. Is there something more recent?

I am also looking for software which can account for **indels**.

# Answer
> 2 votes
**Question** I suspect this question is about the phylogenetics of methylation and the approach the investigator is proposing would be the last approach to use.

**Summary** The approaches to assess the phylogenetics of methylation in order of preference are:

* dN/dS between CpG sites and non-CpG sites,
* An explicit molecular clock between CpG sites and non-CpG sites
* A likelihood ration test (LRT) based on a null distribution of mutations generated by a Monte Carlo algorithm

**Background/Rationale** The approaches you are using have been around for a long while and generate the random distribution of mutations for a given phylogeny. Its phylogenetics in reverse and the technique is known as Monte Carlo, you start with a randomised probability and send it through a parameterised model to predict the amino acid/nucleotide. Thus it is ML(max likelihood), Bayesian, "phylogenetic-HMM" in reverse. It is used within likelihood ratio test calculations. The model is determined by a standard ML, Bayesian phylogeny algorithm, i.e. its circular mathematics because there is no independent means of calculating the mutation behaviour, so the precise context you are using this for needs to be carefully considered.

**Packages** There are lots of packages that do this type of Monte Carlo simulation, SOWHat being a good one (LRT) and you generate a large number of replicate data sets (100 or 1000). The one of the base algorithms is "seq-gen", although PAML may have implemented this.

**Considerations** To use this style of approach you need to carefully consider your question. The "mutation simulations" when analysed via a phylogenetic program will produce the same parameters that you initially set and trees of basically the same length. If you are using this to generate a null distribution for a phylogenetics test these approaches are useful, you then compare the observed likelihood against the null distribution. If you are using it to work out whether the mutation rate at CpG sites is higher than at other sites, its one approach amongst a number of alternatives.

**Drawbacks** Calculating a de novo null distribution of the mutation rate is computationally very expensive and therefore tends to the last calculation performed. It will address a single hypothesis, usually regarding topology.

**Crux of the problem** To really do what (I think) you want to do you need an independent measure of mutation and mutational behaviour and its not trivial to achieve. You would have to consider it like a machine learning calculation, with a formal training set, I don't know whether such an approach has been implemented.

**Strengths/Summary** In summary, Monte Carlo phylogenetics simulation has a limited application because of circularity of the calculation, BUT, BUT, BUT what I have omitted is when it is appropriate it is a very powerful test indeed.

---

**Points of interest** I do like the population dynamic simulations within Beast, this certainly "a priori" and its potential IMO needs further exploration. However, I don't think you looking at how molecular epidemiology might impact the mutational behaviour methylation. I have not looked at Indelible, but looks interesting.

I was going to talk about "processor war simulations" but its probably off-topic.

---
Tags: phylogenetics, phylogeny, mutations, nucleotide-models
---


# Question
Title: How do you download stuff from NCBI fast?
I have 2 .csv files that contains a list of accession codes. For example for this experiment a .csv file will contain:

```
SAMN10870244
SAMN10870243
SAMN10870242
SAMN10870241
SAMN10870240
SAMN10870239
SAMN10870238
SAMN10870237
SAMN10870236
SAMN10870235
SAMN10870234
SAMN10870233

```

Then I run the following

```
for i in $(cat *.csv); do echo ${i}; /usr/bin/sratoolkit/bin/fastq-dump --gzip ${i} & done

```
1. Some of my downloads takes a lot of time days. Is there a faster way to get stuff from NCBI?
2. Are there mirrors of NCBI to accelerate the process?
3. Why is it so slow?

# Answer
Answers to this question can be found on biostars. To give you a short summary:

1. Faster way might be to use the `parallel-fastq-dump`, as suggested in this answer. I never tested that tool though, my own experience is that `prefetch` is more stable than `fastq-dump` command.
2. In that same post it seems that the Japanese mirror seems faster, don't know if that is still valid.
3. The reason this is never optimized might be that downloading from `SRA` is not a task you'll have to do repeatedly, as described here in this post.
> 6 votes

# Answer
You can try `wget` to download SRA files from NCBI server. This seems to be 2x faster in terms of CPU time than `prefetch` and has the advantage of not clogging your `/home/ncbi/public/sra/` directory in case you have space limitation. The tutorial for `wget` downloading method is available at https://github.com/johanzi/ncbi\_tutorial.

Here my timing of the `prefetch` method:

```
time prefetch SRR3300113
real    3m5.609s
user    0m6.332s
sys     0m2.076s

```

Versus the `wget` method:

```
time wget ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR330/SRR3300113/SRR3300113.sra
real    3m4.696s
user    0m0.456s
sys     0m3.332s

```

That makes a user+sys time of 8.402s for `prefetch` and 3.788s for `wget`. Not so much difference for real time though (about a second).
> 1 votes

---
Tags: phylogenetics
---


# Question
Title: UCSC packages/tools ucsc-netchainsubset and ucsc-fatotwobit not working on GNU/Linux
I have installed some ucsc tools on my Mac computer and it worked fine:

### Installation

```
conda create --name ucsc-netchainsubset377 --channel bioconda --yes ucsc-netchainsubset=377

```

### Conda environment activation

```
source activate ucsc-netchainsubset377

```

### Checking the tool's help

```
    netchainsubset
netChainSubset - Create chain file with subset of chains that appear in the net
usage:
   netChainSubset in.net in.chain out.chain
options:
   -gapOut=gap.tab - Output gap sizes to file
   -type=XXX - Restrict output to particular type in net file
   -splitOnInsert - Split chain when get an insertion of another chain
   -wholeChains - Write entire chain references by net, don't split
    when a high-level net is encoundered.  This is useful when nets
    have been filtered.
   -skipMissing - skip chains that are not found instead of generating
    an error.  Useful if chains have been filtered.

```

The same commands on GNU/Linux produced the following error:

```
netchainsubset
netchainsubset: command not found

```

even though the tool seems to be installed successfully:

```
Downloading and Extracting Packages
ucsc-netchainsubset- | 310 KB    | ############################################################################ | 100%
mysql-connector-c-6. | 4.4 MB    | ############################################################################ | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

```

The same happened with `ucsc-fatotwobit` and `ucsc-liftup` and I would not be surprised if this behaviour affected the whole suite of ucsc programmes installed through conda.

Anyone can suggest a solution?

### Update 5 minutes after posting

Thanks to terdon, I printed the PATH when the conda environment was active and listed what was in the folder where conda installed the UCSC tool.

```
ls /mnt/home2/miska/cr517/anaconda3/envs/ucsc-netchainsubset377/bin
c_rehash         libpng-config      mysql_config    openssl  pngfix
libpng16-config  my_print_defaults  netChainSubset  perror   png-fix-itxt

```

Then I typed `netChainSubset` and got the help page, suggesting the tool is installed and working. Interestingly, on Mac, the tool is also called `netChainSubset` but it seems it can be called with `netchainsubset`or `NETCHAINSUBset` for that matter. The case doesn't seem to matter on mac. To make my code more portable, I am going to follow the case in the `bin` of where the tool is installed.

Cheers terdon.

Cheers.

# Answer
## Short answer: respect the case of the tool name

## Longer story below

Thanks to terdon, I printed the PATH when the conda environment was active and listed what was in the folder where conda installed the UCSC tool.

ls /mnt/home2/miska/cr517/anaconda3/envs/ucsc-netchainsubset377/bin c\_rehash libpng-config mysql\_config openssl pngfix libpng16-config my\_print\_defaults netChainSubset perror png-fix-itxt

Then I typed netChainSubset and got the help page, suggesting the tool is installed and working. Interestingly, on Mac, the tool is also called netChainSubset but it seems it can be called with netchainsubsetor NETCHAINSUBset for that matter. The case doesn't seem to matter on mac. To make my code more portable, I am going to follow the case in the bin of where the tool is installed.
> 1 votes

---
Tags: software-installation, ucsc, bioconda, conda, anaconda
---


title: How do I find split reads?

text:
# Question
Title: How do I find split reads?
How can I detect a split read in a BAM file? Is there any sign in the CIGAR string that describes split read?

# Answer
The answer to this question depends on the aligner used to generate SAM/BAM file. For BWA method that is part of the pipeline mentioned in comments, split reads will be flagged as supplementary alignment (2048). You can modify this behavior by using option -M as described in the BWA-MEM manual page.

So, in case of split alignment, you will end up with multiple aligned segments that have identical query name and they are flagged as 2048. These segments may only have part of the original read sequence, but the complete read sequence is present in at least one of the segments.

Also, there is no unique sign in CIGAR of an individual segment to identify split reads. Each aligned segment of the original read may have long deletion in its CIGAR or aligned pairs. However, without looking at all segments with the read ID it is not clear whether the segment has a novel sequence in it that doesn't align with reference or that part matches different chromosome/coordinates.
> 3 votes

---
Tags: sequence-alignment, sam, structural-variation, cigar
---


title: Repeatscout installation problem

text:
# Question
Title: Repeatscout installation problem
I've installed Repeatscout from https://bix.ucsd.edu/repeatscout/. I've also installed TRF and nseg required by the program. I followed the steps:

```
build_lmer_table -l 16 -sequence ~/file_fasta/od_1000.fa -freq od_1000.freq

```

```
 RepeatScout -sequence od_1000.fa -output od_1000_repeats.fa -freq od_1000.freq -l 16

```

and the I put this command:

```
cat od_1000_repeats.fa| filter-stage-1.prl >od_1000_repeats_filt1.fa

```

but this message appears:

```
bash: filter-stage-1.prl: command not found...

```

Why? It is all correctly installed, I can't understand the problem

# Answer
If you're in the RepeatScout directory you need to execute the Perl script with a `./` prefix, like so.

```
cat od_1000_repeats.fa | ./filter-stage-1.prl > od_1000_repeats_filt1.fa

```

If you're not in the RepeatScout directory, you need to provide the full relative or absolute path when executing the script. Alternatively, you can put the RepeatScout directory in your `$PATH` environmental variable.

These are all standard UNIX/Linux concepts that many bioinformaticians take for granted, for better and for worse. If you're lost in all the terminology, a crash course in UNIX basics is probably where you would want to start.
> 2 votes

---
Tags: software-installation
---


# Question
Title: I get the following error when using the GSEA software 'Enrichment Map Visualization' tool using the Cytoscape 3.7.1 application
I am trying to visualize a GSEA analysis I conducted using the Broad Institute's GSEA Enrichment Map Visualization, however I get the following error and am not sure what the issue is: any help would be greatly appreciated.

```
INFO - successfully connected to cytoscape rest
INFO - successfully connected to cytoscape rest
INFO - Found enrichment map command
INFO - Get URL:http://localhost:1234/v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt
DEBUG - CookieSpec selected: default
DEBUG - Auth cache not set in the context
DEBUG - Connection request: [route: {}->http://localhost:1234][total kept alive: 0; route allocated: 0 of 2; total allocated: 0 of 20]
DEBUG - Connection leased: [id: 4][route: {}->http://localhost:1234][total kept alive: 0; route allocated: 1 of 2; total allocated: 1 of 20]
DEBUG - Opening connection {}->http://localhost:1234
DEBUG - Connecting to localhost/127.0.0.1:1234
DEBUG - Connection established 127.0.0.1:51968<->127.0.0.1:1234
DEBUG - Executing request GET /v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt HTTP/1.1
DEBUG - Target auth state: UNCHALLENGED
DEBUG - Proxy auth state: UNCHALLENGED
DEBUG - http-outgoing-4 >> GET /v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt HTTP/1.1
DEBUG - http-outgoing-4 >> Host: localhost:1234
DEBUG - http-outgoing-4 >> Connection: Keep-Alive
DEBUG - http-outgoing-4 >> User-Agent: Apache-HttpClient/4.5.2 (Java/1.8.0_201)
DEBUG - http-outgoing-4 >> Accept-Encoding: gzip,deflate
DEBUG - http-outgoing-4 >> "GET /v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt HTTP/1.1[\r][\n]"
DEBUG - http-outgoing-4 >> "Host: localhost:1234[\r][\n]"
DEBUG - http-outgoing-4 >> "Connection: Keep-Alive[\r][\n]"
DEBUG - http-outgoing-4 >> "User-Agent: Apache-HttpClient/4.5.2 (Java/1.8.0_201)[\r][\n]"
DEBUG - http-outgoing-4 >> "Accept-Encoding: gzip,deflate[\r][\n]"
DEBUG - http-outgoing-4 >> "[\r][\n]"
DEBUG - http-outgoing-4 << "HTTP/1.1 500 Internal Server Error[\r][\n]"
DEBUG - http-outgoing-4 << "Date: Sun, 17 Mar 2019 06:55:04 GMT[\r][\n]"
DEBUG - http-outgoing-4 << "Content-Type: text/plain[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Origin: *[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Headers: origin, content-type, accept, authorization[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Credentials: true[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Methods: GET, POST, PUT, DELETE, OPTIONS, HEAD[\r][\n]"
DEBUG - http-outgoing-4 << "Content-Length: 20[\r][\n]"
DEBUG - http-outgoing-4 << "Server: Jetty(9.4.10.v20180503)[\r][\n]"
DEBUG - http-outgoing-4 << "[\r][\n]"
DEBUG - http-outgoing-4 << "Failed: empty String"
DEBUG - http-outgoing-4 << HTTP/1.1 500 Internal Server Error
DEBUG - http-outgoing-4 << Date: Sun, 17 Mar 2019 06:55:04 GMT
DEBUG - http-outgoing-4 << Content-Type: text/plain
DEBUG - http-outgoing-4 << Access-Control-Allow-Origin: *
DEBUG - http-outgoing-4 << Access-Control-Allow-Headers: origin, content-type, accept, authorization
DEBUG - http-outgoing-4 << Access-Control-Allow-Credentials: true
DEBUG - http-outgoing-4 << Access-Control-Allow-Methods: GET, POST, PUT, DELETE, OPTIONS, HEAD
DEBUG - http-outgoing-4 << Content-Length: 20
DEBUG - http-outgoing-4 << Server: Jetty(9.4.10.v20180503)
DEBUG - Connection can be kept alive indefinitely
INFO - status:Internal Server Error
INFO - Unable to create Enrichment Map:Internal Server ErrorResponseEntityProxy{[Content-Type: text/plain,Content-Length: 20,Chunked: false]}
DEBUG - Connection manager is shutting down
DEBUG - http-outgoing-4: Close connection
DEBUG - Connection manager shut down

```

# Answer
This is a tough one, because the only relevant lines from your GSEA debug log were these two:

```
DEBUG - http-outgoing-4 << "Failed: empty String"
DEBUG - http-outgoing-4 << HTTP/1.1 500 Internal Server Error

```

My recommendations to you at the time would've been

1. try closing all the programs, restarting them, possibly rebooting and restarting them, and even analyzing a completely different set of data to see if you get the same results, **—or—**
2. update to the latest Java 8, reinstall Cytoscape 3.7.whatever and the EnrichmentMap plugin from scratch, and re-download GSEA, **—or—**
3. lastly, try reinstalling everything on a completely different machine/OS and see if you get different results (*e.g.,* try on a Mac or Linux if you had problems with Windows originally, or vice-versa)

I know that's lame, but sometimes that's how you solve weird computer errors, when the error message the application itself gives you is non-specific and unhelpful.

If you had seen the error message

```
None of the gene sets have passed the filter. Try relaxing the gene set filter parameters.

```

in the GSEA application logs, then you would've had the problem described in EnrichmentMap's GitHub issue #358, and the solution would be to use less stringent p- and FDR Q-value cutoffs, as prescribed in the "Tips on Parameter Choice" section of the EnrichmentMap documentation.

*My* problem was that GSEA's default cutoffs (p-value \< 0.005 and FDR Q-value \< 0.1) were entirely too stringent for my results set, and using the "Very permissive" settings of p-value \< 0.05 and FDR \< 0.25 from the documentation solved it for me.

Hopefully this benefits someone who ends up here after googling for "Unable to create Enrichment Map:Internal Server ErrorResponseEntityProxy," which is just a stupendously unhelpful error message when "None of the gene sets have passed the filter" is *actually* what went wrong. Maybe go thumbs-up this comment in the EnrichmentMap issue tracker, if you think the maintainers need to work on that.
> 1 votes

---
Tags: errors, gsea
---


# Question
Title: PAML branch-site model: Identical lnL value between null and alternative models?
I have been running PAML, via codeml, to try to identify positive selection in certain branches of my tree for one gene (but will eventually be trying to automate the analysis over multiple genes). This was done by setting `model = 2` and `NSsites = 2` in the codeml.ctl file. I attempted to run both the null and alternative models by setting `fix_omega = 1` in the null.ctl and `fix_omega = 0` in the alternative.ctl.

However, after running both models, they have returned exactly identical log-likelihood (lnL) values. I wasn't expecting to find anything significant, but the fact that they were identical struck me as odd. Out of curiosity, I altered the branch of interest in the newick tree file and reattempted the analysis, but again, the lnL values remain the same between the null and alternative models (although slightly different to the value returned with the different branch selected).

I have since repeated this many times, choosing a different branch(es) each time to see if I can get a difference between the null and alternative, but no luck. Am I missing something obvious? Should I really be getting the same exact result between null and alternative models, regardless of where in the tree I place `#1`? To reiterate, I wasn't expecting to find significant results, but was at least expecting (non-significantly) different numbers for lnL...

# Answer
> 1 votes
The issue with PAML is that it doesn't throw many exceptions (programming jargon for trapping and reporting errors) if something is wrong, generally the "exceptions" are identified in the failure to report an expected value in the output.

**Possible explanation** fix\_omega is not turning the omega parameter on or off it is simply defining what the value of omega will be. If the default value, which is set at .4, but can be changed, is very close to the estimated value the likelihoods will be similar - or even identical.

**Goto debugging** As computer/data scientists, the very first thing we should do is just check there is run time variation between the options you are assessing via,

```
time codeml ./file

```

The thing to remember is fix\_omega=0 is where Omega is being estimated from the data, whilst =1 means it will be supplied with a value. There should be longer run time on the "0" option.

**Nitty gritty debugging**

1 Please check you have

```
seqtype = 3

```

If not, then that is your answer, omega was never being set.

2 For debugging only make sure you have,

```
model = 1

```

This is global dN/dS for all the data not for a specific clade within the tree, the reason for this is given in the final step below

3, Then try a different NSsites model,

```
NSsites = 8

```

This should produce real differences. The NSsites options are,

> * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs; * 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
> * 10:beta&gamma+1; 11:beta&normal\>1; 12:0&2normal\>1; * 13:3normal\>0

w here is for omega

4 The proof of your system is that when you

```
fix_omega = 1
omega = .4 * try lots of different values, e.g. 0.3, 0.4, 0,7, 1, 1.5, 2

```

This should produce likelihood variation

---

Finally

```
model = 2

```

Then go to your tree file and make sure that the nomenclature used follows the manual, I think from memory the # is deployed. The output file should change to specify that you have two values for dN/dS from your treefile, thus the output from model=2 should vary from the output of model=1 if not the delineation has been ignored. If you do see two rates, great and now you are into parameter exploration. Parameter exploration is basically the options 1 to 4 described above, plus the stuff you mentioned.

---

Oh you mentioned using Newick format. I haven't used PAML for a while (dNdS stuff is ancient) but it it only used phylip format when I used it. I suspect if it got the wrong tree format PAML will throw an exception, but it's worth checking.

---
Tags: positive-selection
---


title: What does PCA mean on GWAS

text:
# Question
Title: What does PCA mean on GWAS
I understand what GWAS is and I'm able to perform certain tests with the p-values, etc. But what I am having a hard time wrapping my head around is what PCA on GWAS means.

So let's say I have 100,000 individuals and genotype data for 10 million SNPs. And I'm analyzing this for 10 different diseases. I also have p-value information for each of the 10 diseases for the 10 million SNPs.

At my work, they were throwing around the following phrases. I'm not quite sure what exactly it means.

1. "Perform PCA on GWAS." Does this mean, I'm performing PCA on the genotype data or p-value data? I've done PCA on classic machine learning problems like "iris." But I'm not able to quite wrap my head around what PCA information tells me in terms of this genetic data.
2. "PCA on individual diseases vs. PCA on entire GWAS" What exactly is the difference. I realize this is a broken statement, but they said we have two sets of PCs. One on the individual diseases which are unique, and then PCs on the entire GWAS for all 10 diseases.

Sorry if my questions is not properly asked, I realize it's probably heavily abusing notation.

**Added clarifications:**

1. So, my matrix has 100k rows and 10 million columns. Then, PC1, for example, is a vector with 100k values, each value corresponding to each individual. What exactly does PC1 value that's associated with the first individual tell me with regards to the PC1 value that's associated with 2nd individual, and so on. My current understanding if you think of a number line that represents PC1, then each of the PC values (corresponding to each individual) represents where on the number line that individual falls.
2. How exactly does this eliminate/account for variation in population stratification?

# Answer
From my memory of what a statistician told me, a PCA aims to determine independent linear combinations of variables (i.e. genotypes) that account for the most variation in the dataset. With 10 million SNPs, the vectors describing the genotypes of all individuals reside somewhere inside a 10-million-dimension hypersphere, and the first principle component determines the best vantage point / angle from outside that sphere that will create the greatest spread of data when data is projected onto a perpendicular / orthogonal axis. The second principle component is the perpendicular / orthogonal axis that will create the second-greatest spread of data.

*Aside: There's some wonderful linear algebra (which I've forgotten) that makes it very quick (computationally) to do this angle calculation, exploiting the nature of linear transformations and orthogonality.*

The point of the PCA for GWAS is to identify (and preferably exclude) differences that might be due to the structure of populations, rather than the disease of interest. It's fairly common for researchers to do a PCA on the entire set of genotype data, then use the first 5, 10, or 20 principal components as covariates in the association model.

There's an underlying assumption that the majority of genetic variation is due to population structure - if this is not the case, then the disease would have to be extremely omnigenetic, with almost the entirety of the genome contributing to the disease. The approach also won't deal with *all* the population structure (because there may be non-linear contributions of genotypes to population structure), but it will get rid of the largest problems.

Additional to this, I'd recommend the following things:

# Structure Visualisation

Carry out a STRUCTURE analysis or PCA on the combined cases and controls, and display the results of the first two principle components, labelling, colouring, or otherwise indicating both groups. As an additional check, it may be useful to do addititional PCAs for cases and controls separately.

If they appear different, it suggests that the cases and controls were drawn from separate population clusters, which should be commented on in the discussion of the paper. There's no need to panic and find more data, just make the observation.

# Subsampling

Carry out association tests on an equal-sized subsample of both datasets. This will mean that area-under-the-curve statistics can be properly interpreted. If you've got unbalanced datasets, use a little bit less than the smaller of the two groups (e.g. for 5k diseased and 95k healthy, use a random subsample of 4.5k from both groups).

# Multiple testing / Internal validation

Run the association tests multiple times (with different subsamples), recording the individuals that were subsampled for future reference / reproducibility.

If an association is only present in one of the many different subsamples, it should be ignored as a spurious result. It means that there is some association that is not generalisable to all (or most) of the cases and controls. These variants might be interesting in themselves on an individual level, but not in a general population.

I explain some of this in more detail in the following preprint (which I'm not planning to submit to any journal, due to it being a bit stale and needing a lot of work to fix up):

Bootstrap Distillation: Non-parametric Internal Validation of GWAS Results by Subgroup Resampling
> 5 votes

# Answer
PCA = principle component analysis and a multivariate statistic, today it is trendily retermed "unsupervised learning" and here is likely being deployed for individuals within your data set. It works by identifying the maximum variance within multidimensional space, shearing it and describing this as the first principle component. The second principle component is sections the variance at 90o to the first in all axes (orthogonal is the posh word), and so on ... The maximum information is within the first two components and they are sometimes plotted on x and y axes to look for potential clusters within your data. In this case the first few components are all that is needed because you not identifying the genes of interest just checking the data looks sound.

**Point 1** PCA on all GWAS is attempting to assess that the genetic variation between diseased and healthy people is evening spread across the genome. For example if your diseased population came from e.g. the indigenous people of Western Samoa (purely as an example), which has a very defined genetic heritage (at least for HLA) then most of the SNPs will be unrelated to the disease your are attempting to identify. This will mean you have to redefine the healthy population, the null distribution, to screen out the artefacts. In this case a better representation would be to select healthy individuals from Western Samoa as your comparison. So within the first two components you want to see lots of mixing between your patients vs healthy people. If the patients are clumping (clustering) thats bad here and you have to find out why.

PCA has weaknesses BTW for complicated reasons which are boring, but it is often not assessed directly but via a tSNE analysis and then plotted on two axes.

**Point 2** If PCA from point 1 looks cool, then you might split the data set into diseased and another one of healthy and perform a PCA on each but now look at the SNPs ... not the people. I hope that makes sense. If you compare the PCA plot of SNPs between these two groups, the differences MIGHT and I stress might highlight the SNPs of interest, i.e. the SNPs causing your disease. However and its a big however the number of SNPs you are dealing with is huge and if you are manually going through all that data to assess anomalous clusters, or outliers to make a list of candidate SNPs for further analysis (@gringer describes some of the further options above) - I wish you good health and long life, and you will need it because you will be very old at the end of that optical analysis.

What you really do is an additional analysis (I would deploy tSNE) to assess the differences between the two PCA plots and end up with a third plot indicating the differences between all the other PCA plots. In this context it is complicated to deploy and would be nested within a tSNE analyse of all components which is what I alluded to in point 1.

**Point 2 - issues** Generally if someone thought I was going to eye-ball 50000 SNPs between two graphs, or throw them back a reanalysis which conceptually made their brain hurt ... their gonna have conceptualise real hard.

**Summary** Having said all that PCA is merely an exploratory tool and is not the defining significance test required, to find the "one SNP" which will explain your disease.

---

The imbalance of the data set makes advanced SNP identification statistically difficult, but the 30:70 split is okay. However a subsample of the healthy population, which is what Point 1 is all about, would likely correct for this imbalance.

---

On the question of mixing Point 1, what you don't want to see is your diseased patients forming a "blob", or even a circle, or any sort of defined cluster against the spread of healthy patients on a scatterplot of PC1 Vs PC2. Perhaps I should have said if your patients don't cluster on the output to Point 1 that's good. There are statistics that identify clusters on a PC1 PC2 scatter plot, and you want it to fail that test.
> 2 votes

---
Tags: phylogenetics, phylogeny, snp, gwas, pca
---


# Question
Title: Generate table for total number of SV events per sample
I have called and annotated structural variants and I have a table by name of samples and structural events like below

```
Sample  Type
t_005   DEL
t_005   INS
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DUP
t_005   DUP
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DUP
t_005   BND
t_005   DEL
t_005   DEL
t_005   BND
t_005   DEL
t_005   DEL
t_005   DUP
t_005   DEL
t_005   DEL
t_005   DEL
t_005   DEL

```

How can I count the number of events to get the below table?

```
> head(SV_Counts[,1:4])
               t_005      o_005       t_013         o_013
    DEL          27          40          25          30
    DUP          11          10          13          22
    BND          23          20          10          12
    TRA           3           2           4           4
    >

```

I tried this from google but I failed to adapt my data to that

```
    > melt(table(a))
   Sample Type value
1   o_036  BND   280
2   o_005  BND     8
3   o_013  BND   160
4   o_021  BND   144
5   o_036  BND     0
6   t_005  BND    15
7   t_013  BND   108
8   o_036  DEL    63
9   o_005  DEL    15
10  o_013  DEL   121
11  o_021  DEL   243
12  o_036  DEL     0
13  t_005  DEL    40
14  t_013  DEL    76
15  o_036  DUP    34
16  o_005  DUP     4
17  o_013  DUP    57
18  o_021  DUP    88
19  o_036  DUP     0
20  t_005  DUP    14
21  t_013  DUP    48
22  o_036  INS    58
23  o_005  INS     0
24  o_013  INS     4
25  o_021  INS    10
26  o_036  INS     1
27  t_005  INS     1
28  t_013  INS     2

```

But hereafter I don't know what to do

# Answer
So you have data in "long format" and want to convert it to "wide format". You can use `spread()` from `tidyr` for that:

```
require(tidyr)
# Assuming your input is called `foo`
spread(foo, Sample, value)

```
> 2 votes

---
Tags: r, structural-variation, wgs
---


# Question
Title: Delly finds identical deletions with different consensus sequences
I am using delly to call structural variants on a human WGS sample. I noticed that it finds certain variants that are essentially duplicates. For example:

```
chr4    177594033       DEL00005339     G       <DEL>   .       PASS    PRECISE;SVTYPE=DEL;SVMETHOD=EMBL.DELLYv0.8.1;CHR2=chr4;END=177594095;PE=0;MAPQ=0;CT=3to5;CIPOS=-3,3;CIEND=-3,3;SRMAPQ=60;INSLEN=0;HOMLEN=3;SR=5;SRQ=0.993243;CONSENSUS=GACAAAAACATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAAC;CE=1.75466     GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-51.3797,0,-42.791:10000:PASS:1115:1821:1820:1:0:0:23:27
chr4    177594033       DEL00005340     G       <DEL>   .       PASS    PRECISE;SVTYPE=DEL;SVMETHOD=EMBL.DELLYv0.8.1;CHR2=chr4;END=177594095;PE=0;MAPQ=0;CT=3to5;CIPOS=-3,3;CIEND=-3,3;SRMAPQ=60;INSLEN=0;HOMLEN=3;SR=12;SRQ=0.994186;CONSENSUS=ATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAACAATGTAATAAGTATGTTTGAATAGTATTATGAT;CE=1.72637    GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-51.3797,0,-42.791:10000:PASS:1115:1821:1820:1:0:0:23:27

```

The two variants above have the same start position and the same end position (`END=177594095`). However, they differ in the `SR`, `SRQ` and `CONSENSUS` fields. These are defined in the VCF header as:

```
##INFO=<ID=SR,Number=1,Type=Integer,Description="Split-read support">
##INFO=<ID=SRQ,Number=1,Type=Float,Description="Split-read consensus alignment quality">
##INFO=<ID=CONSENSUS,Number=1,Type=String,Description="Split-read consensus sequence">

```

Of course, the consensus sequences are quite similar, but they aren't identical (I have simply added some spaces to the beginning of the 2nd sequence here to align them):

```
GACAAAAACATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAAC
         ATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAACAATGTAATAAGTATGTTTGAATAGTATTATGAT

```

So what am I looking at here, exactly? If this is a deletion on chromosome 4 from position 177594033 to position 177594095, then how can the length of the consensus sequence be different? This is the "split read consensus", so presumably the sequence that appears to have been deleted based on how the reads covering that region are split. This means that these reads align to the genome with a gap and this gap indicates the deletion (image source):

Therefore, if the start and end position are the same, the consensus should also be identical. And, in any case, if the start and end positions are the same, then this should be the same structural variant, so why is it being reported twice? How can I interpret this result?

# Answer
It turns out that this was just a bug which has been fixed by this commit. The current `master` version of delly no longer exhibits this behavior.

If you are facing this issue, either build the current `master` version or wait a bit for the next release which should also include this fix.
> 1 votes

---
Tags: structural-variation
---


# Question
Title: Frequency of specific viral sequence in .BAM or .fastq
I was wondering if anyone had any experience 'counting' the frequency of a particular viral sequence in an individual's fastq sequence. So basically, counting the number of occurrences of a particular sequence given someone's fastq file.

The thing is, we must speak with the sequencing company regarding what format we'd like the files and must choose between `.bam` and `.fastq` \- both would be extra $. Usually I'd choose `.bam` out of convenience but I was wondering, given that it has already been aligned the human genome I might not be able to proceed with this 'counting' process. So, my question is: if I were to use `samtools` \- or any other program - to convert the `.bam` files to `.fastq`, is there any possibility that I could lose such information - e.g. viral seqs - that would have otherwise been available should I have just gotten the `.fastq` files directly?

# Answer
Most read aligners will report unaligned reads as well, which presumably will include your viral sequences. I would ask them to formally confirm that the BAM files will contain unaligned reads before choosing that option.
> 6 votes

# Answer
If you are interested in a particular viral sequence, it seems wasteful to request bams that are aligned to the human genome, but you could request bams aligned to the human genome with unaligned reads retained. As an alternative, you can generate unaligned bams, where all of the same reads are present that would be present in a fastq but in a bam format. You can revert this to a fastq or align the fastq directly to a viral reference yourself. Aligning to small references is fast.
> 3 votes

---
Tags: sequence-alignment, samtools, fastq, covid-19, sars-cov-2
---