title: Output from vcftools missingness

text:
# Question
Title: Output from vcftools missingness
I'm new to data filtering on vcf data and vcftools.
I performed variant calling on my dataset, CHR22, homo sapiens. I'd like to remove sites that are missing in more than 5% of individuals.

```
vcftools --missing-site --vcf updated_ids68.vcf

```

#

This gives me a file with out.lmiss. That is missingness per locus.
I had 68 individuals in my input vcf file.

```
CHR     POS     N_DATA  N_GENOTYPE_FILTERED     N_MISS  F_MISS
22      16848278        106     0       30      0.283019
22      16848492        68      0       68      1
22      16849180        69      0       67      0.971014
22      16849229        68      0       68      1
22      16849376        133     0       3       0.0225564
22      16849476        132     0       6       0.0454545
22      16851734        126     0       36      0.285714
22      16852588        123     0       13      0.105691

```

#

I'm unable to understand the output here. What does third column - *N\_DATA* tell?
133, 126 doesn't tell that this is number of individuals if my input had 68 individuals only.

VCFtools (0.1.15)
fileformat=VCFv4.1

# Answer
> 4 votes
Humans are diploid, so you can expect to see up to 2\*N (2\*68=136) alleles, so N\_DATA is the number of observances of that allele.

---
Tags: variant-calling, variation, filtering, vcftools
---


title: Pathway plot in GeneTrail data base

text:
# Question
Title: Pathway plot in GeneTrail data base
So I'm doing a pathway analysis using GeneTrail database. One thing I am not able to get is this similarity measure being calculated and plotted:

I want to know how this is being made, what is the output being plotted because each gene is may be involved in multiple pathways.

I want to be able to make it in R again. Any suggestion or help would be appreciated. How to plot the enrichment result in this form ?

# Answer
> 3 votes
This is a a symmetric matrix with the pathways in both axes. So for each pair a similarity index is calculated (usually it takes into account the number of genes in each pathway and the total number of genes in both pathways). There are several similarity indexes available. A typical one is Jaccard or Sørensen-Dice index, usually 0 means no similarity and 1 equal pathways.

You can do this kind of analysis in R using BioCor (Disclosure: I am the author of this package) which uses the Dice similarity. You will need the name of the genes of each pathway, after calculating the similarity of pathways using `mpathSim`, you can then plot it with a heatmap (using the `heatmap` function).

---
Tags: r, rna-seq, pathway, figure-explanation
---


# Question
Title: Do any computational phylogenetic methods enable the specification of ancestral states?
Various phylogenetic algorithms estimate ancestral states of a phylogenetic dataset. Is there a way in either maximum parsimony, distance-based methods, or Bayesian inference to indicate what the ancestral states of characters were?

# Answer
> 2 votes
For maximum parsimony, Fitch's method seems simple and intuitive, and can be done in the R phangorn package: Fitch, Walter M. (1971). "Toward Defining the Course of Evolution: Minimum Change for a Specific Tree Topology". Systematic Zoology. 20 (4): 406. doi:10.2307/2412116. ISSN 0039-7989. JSTOR 2412116

There are various methodological criticisms of parsimony that you can find in the literature, but its appeal is you can easily understand why the algorithm came to the conclusions it did.

RAxML can do maximum likelihood reconstructions based on a tree and tip values, and I think Mr. Bayes can do a Bayesian version.

# Answer
> 1 votes
Yes this approach is done in parsimony, i.e. you can indicate the outgroup and the ancestral state is calculated from there. I assume this is also performed using 'Beast' because it is a non-reversible matrix. Regular ML and Bayes approaches use reversible mutation matrices therefore the tree is calculated first, you can then specify an outgroup to calculate the ancestral state. However, in a reversible matrix the ancestral state is not used in the calcuation. I hope that makes sense.

---
Tags: software-recommendation, phylogenetics, algorithms, phylogeny
---


# Question
Title: Plotting coverage of annotation over collection of region
I'm trying to plot "meta" coverage of annotation: i.e. features (eg. gene class) over certain regions. It is similar to read coverage plots over gene body, except my input is two bed files (both in BED6 format) - (A) one containing the regions for which the meta plot is needed, and another being the list of genes or any feature whose coverage over the regions in (A) is needed.

Is there any package or tool which can create such plots (my domain is limited to python but I can try to work with R) ?

Something akin to this but for whole gene body\[1\]:

\[\]- Kelley, D, Rinn, J (2012). Transposable elements reveal a stem cell-specific class of long noncoding RNAs. Genome Biol., 13, 11:R107.

Edit: The BED file containing list of regions looks like this:

```
chr3    39218734    40053659    region1 0.92426187419769    +
chr4    140163762   140453127   region2 0.896103896103896   -
chr7    40549151    41205036    region3 0.986072423398329   +
chr8    81291743    81963246    region4 0.94184168012924    -
chr9    12284032    12539789    region5 0.95539033457249    -

```

And the bed file containing features to be plotted looks like this:

```
chr3    39218100    40053200    LINE    1   +
chr4    140163962   140453027   LINE 1  -
chr7    40549002    41204999    SINE 1  +
chr8    81291143    81963846    LTR 1   -
chr9    12284332    12539720    LTR 1   -

```

# Answer
> 3 votes
If one assumes that repeat elements of a given type (e.g., LINEs) don't overlap each other, then the following will work:

1. Split your BED file by repeat element, such that you have a LINE.bed, SINE.bed, etc.
2. Convert those to bedGraph (e.g., `awk 'BEGIN{OFS="\t"}{print $1,$2,$3,"1.0"}' LINE.bed > LINE.bedGraph`).
3. Use UCSC tools to convert those bedGraph files to BigWig.
4. Install deepTools and run `computeMatrix reference-point -b 2500 -a 20 -S LINE.bigWig SINE.bigWig LTR.bigWig -R Regions_Of_Interest.BED -o foo.mat.gz`
5. Make a profile plot with plotProfile (`plotProfile -m foo.mat.gz -o foo.png --perGroup`)

If you do have overlapping regions you'll need to first make a disjoint set of intervals (you can use bedops for this), find the coverage of them (e.g., `bedtools coverage`) and then continue on with that.

# Answer
> 0 votes
There is a reasonably nice looking tutorial on metaplots in R here: https://rpubs.com/achitsaz/94710

---
Tags: python, software-recommendation, visualization, coverage
---


# Question
Title: Which software is used to run Molecular Dynamics simulation of a small RNA hairpin?
I am currently working in a project that needs MD simulation of a small RNA hairpin. I am currently trying it on GROMACS but it seems to be protein based, hence am facing several issues. Is there any software that is designed for the same?

# Answer
> 2 votes
Google found this review that might be of interest: https://pubs.acs.org/doi/10.1021/acs.chemrev.7b00427

Here is a link to a free full-text version: https://www.ncbi.nlm.nih.gov/pmc/articles/pmid/29297679/

The authors "here provide a comprehensive overview of the fast-developing field of MD simulations of RNA molecules. We begin with an in-depth, evaluatory coverage of the most fundamental methodological challenges that set the basis for the future development of the field, in particular, the current developments and inherent physical limitations of the atomistic force fields and the recent advances in a broad spectrum of enhanced sampling methods. We also survey the closely related field of coarse-grained modeling of RNA systems. After dealing with the methodological aspects, we provide an exhaustive overview of the available RNA simulation literature, ranging from studies of the smallest RNA oligonucleotides to investigations of the entire ribosome. Our review encompasses tetranucleotides, tetraloops, a number of small RNA motifs, A-helix RNA, kissing-loop complexes, the TAR RNA element, the decoding center and other important regions of the ribosome, as well as assorted others systems. Extended sections are devoted to RNA–ion interactions, ribozymes, riboswitches, and protein/RNA complexes. Our overview is written for as broad of an audience as possible, aiming to provide a much-needed interdisciplinary bridge between computation and experiment, together with a perspective on the future of the field."

Hope this helps.

# Answer
> 2 votes
AMBER has force fields specifically for viewing RNA free energy and structural dynamics. Specifically, AMBERff99 is useful for RNA analysis. It is fairly easy to use and a comprehensive software for your analysis. It has plenty of tutorials from basic to advanced usage, as well as many publications. A specific publication related to your question can be found at

doi: 10.1016/j.bpj.2009.10.040

But in general, the J Phys Chem B is relevant to MD studies.

---
Tags: rna
---


# Question
Title: How I know which gene is a good predictor in this neural network or not?
I have 25 highly differentially expressed genes among and patients to chemotherapy.

I have made a neural network of these genes. The accuracy of the model is `0.73` but I don't know from 25 genes how I could select the best predictors or remove the weaker ones.

```
> # Convert column from /to 1/0
> = as.numeric(data$)-1
> data1 = cbind(,df)

> head(data1[,1:4])
   ALB CALML5 CALML6
A2       0 4.64   5.77   3.16
A3       0 7.12   9.13   5.13
A4       0 7.84   7.09   4.85
A6       0 7.90   9.64   6.52
A7       0 5.84   6.26   0.00
A8       0 4.83   5.77   2.15
>
> library(caTools)
> set.seed(101)
> # Create Split (any column is fine)
> split = sample.split(data1$, SplitRatio = 0.70)
> # Split based off of split Boolean Vector
> train = subset(data1, split == TRUE)
> test = subset(data1, split == FALSE)
> feats <- names(data[,-1])
> # Concatenate strings
> f <- paste(feats,collapse=' + ')
> f <- paste('TG ~',f)
> # Convert to formula
> f <- as.formula(f)
> f
~ ALB + CALML5 + CALML6 + CDK6 + CHGA + CREB3L3 + CSF3 +
    CXCL5 + CXCL6 + CXCL8 + DKK1 + EGLN2 + FGF20 + FGF4 + IL11 +
    IL1B + IL6 + ISG15 + KLK5 + KRT17 + MMP3 + MMP7 + PKP1 +
    PTGS2 + S100A7A
> #install.packages('neuralnet')
> library(neuralnet)
Warning message:
package ‘neuralnet’ was built under R version 3.5.2
> nn <- neuralnet(f,train,hidden=c(25,2,1),linear.output=FALSE)
> View(test)
> # Compute Predictions off Test Set
> predicted.nn.values <- compute(nn,test[2:26])
> # Check out net.result
> print(head(predicted.nn.values$net.result))
             [,1]
A4  0.26890819169
B5  0.94836576493
B7  0.19753554824
B8  0.99790838807
B9  0.99764782073
B12 0.05230351043
> predicted.nn.values$net.result <- sapply(predicted.nn.values$net.result,round,digits=0)
> table(test$Private,predicted.nn.values$net.result)

    0 1
  0 4 3
  1 1 7

> accuracy<-(table(test$Private,predicted.nn.values$net.result)[1,1]+table(test$Private,predicted.nn.values$net.result)[2,2])/(sum(table(test$Private,predicted.nn.values$net.result)))

> accuracy
[1] 0.7333333333
>

```

# Answer
> 1 votes
Yes, there is a formal solution. You want a heat map of the layers weighting and this is the information sought. It is more difficult to interpret than other supervised learning algorithms, but that is the nature of neural networks.

I know Sci-test Learn (Python), but I don't know R for neural networks, so I can't advise you on the code, only that it will be a few lines of code and a graphical (relative) output. I imagine Tensorflow will have a really brilliant heat map layer weighting.

---
Tags: r, rna-seq, transcriptome, clustering, machine-learning
---


# Question
Title: Restricting match output for multiple sequence alignment using mafft?
So I aligned roughly 5k sequences and I got my output using mafft. However, I want to restrict the output to only present conserved segments of a specific length (let's say 25bp). Does mafft have parameters to set this that I seem to be missing?

# Answer
> 2 votes
## Q - Is there some parameter for this?

No, there is no parameter for this. Getting the conserved regions is something that you will have to perform post-alignment.

## How to get the conserved regions?

You will need to define what exactly you mean by conservation and calculate a conservation score for each position of the alignment.

For instance you can look at each column of the alignment and calculate the percent of queries that had the same base/AA at that position. If using AAs you can define conservation more loosely depending on what properties are important to you, or using a scoring matrix.

Once you have a co-occurrence matrix you can define the regions of interest based on a conservation threshold that you define.

## Is there software to find conserved columns?

There is a program MstatX which outputs conservation scores of individual columns.

## How MstatX works

The example of MstatX provides this example alignment. Each column is progressively less conserved than the last. The last column is conserved, but has some gaps -- this lowers the conservation score.

```
   col          1
 seq1  1234567890
 seq2  DDDDDDIPDL
 seq3  DDDDDDIPVL
 seq4  DDDDDDIPYL
 seq5  DDDDDDIPA-
 seq6  DDDDDDLWT-
 seq7  DDEDEELWK-
 seq8  DDEDEELWP-
 seq9  DDEDEFVSR-
 seq10 DEEFFFVSH-

```

The conservation scores for each column were:

```
1       2.71272e-08
2       0.243433
3       0.312789
4       0.243433
5       0.483364
6       0.480391
7       0.498757
8       0.498757
9       0.973514
10      0.267319

```

## How to get conserved regions

The python program below will output all windows of conserved regions when using a fasta alignment and the default output from `MstatX`.

* Positional args
    + **1.** `MstatX_alignment_example.mali`, is the fasta alignment file.
    + **2.** `output.txt` is the default output of `mstatx` when called with `mstatx -i MstatX_alignment_example.mali`.
    + **3.** 0.4 is the minimum conservation score
    + **4.** 3 is the minimum number of bases that must have that conservation score.

Calling this program with the following command on the above example:

`python print_conserved_regions.py MstatX_alignment_example.mali output.txt 0.4 3` prints:

```
{1,4}

```

Meaning that there was only one window at least 3 bases long in which all positions had a conservation score below 0.4. This window was positions 1 through 4 in the alignment.

```
#call this file cons_windows.py
import sys
from Bio import AlignIO

align_filename = sys.argv[1]
cons_filename = sys.argv[2]
cutoff = float(sys.argv[3])
minchar = int(sys.argv[4])

def parse_conservation(filename):
    cons = {}
    with open(filename, "r") as f:
        for line in f:
            if line.strip():
                temp = line.split()
                cons[int(temp[0])] = float(temp[1])
    return cons

align = AlignIO.read(align_filename, "fasta")
cons = parse_conservation(sys.argv[2])

i = 1
while i < max(cons):
    if cons[i] <= cutoff:
        j = i
        done = False
        while j < max(cons) and not done:
            if (cons[j] > cutoff) and (j-i >= minchar):
                print("{{{},{}}}".format(i,j-1))
                done = True
                i = j
            j+=1
    i+=1

```

---
Tags: sequence-alignment
---


# Question
Title: Any fast options to query large VCF bed intervals?
I'm doing some analysis and I need to subset a large VCF file (~8GB gziped) given a bed interval and identify within a subset of rsid.

Unfortunately, both my normal choices to do this analysis (`snpSift` and `bedtools`) are taking way to long or failing due to memory issues in my local computer and a remote server.

Do you guys know any other options or suggestions to speed up this process?

Follow the commands I use:

```
bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

```

or

```
gzcat <myvcf>.vcf.gz | \
    java -Xmx10g -jar snpSift.jar intervals <mybed>.bed | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

```

The `bedtools` command usually fail due to unknown reason and `SnpSift` runs over 6 hours even given 10GB of ram. My local machine have 8GB of RAM, but the server have 32GB.

# Answer
> 4 votes
this:

```
bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

```

can be replaced with one GATK SelectVariants https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org\_broadinstitute\_gatk\_tools\_walkers\_variantutils\_SelectVariants.php

```
java -jar gatk.jar -T SelectVariants -R ref.fa \
  -L myinterval.bed --keepIDs myrsid.txt -V myvcf.vcf.gz -o out.vcf

```

and furthermore, you can also work in parallel by splitting your bed file (e.g. by chromosomes), and merge the results.

# Answer
> 4 votes
Generically, with BEDOPS:

```
$ vcf2bed < <(gunzip -c snps.vcf) | bedops -e 1 - myRegions.bed > answer.bed

```

Or:

```
$ vcf2bed < <(gunzip -c snps.vcf) | bedmap --echo --echo-map-id --delim '\t' myRegions.bed - > answer.bed

```

Etc.

BEDOPS supports input streams; use them where you can.

That said, if you're going to query the VCF data frequently, convert once and use it in downstream operations.

```
$ vcf2bed < <(gunzip -c snps.vcf) > snps.bed

```

Then:

```
$ bedops -e 1 snps.bed myRegions.bed > answer.bed

```

Etc.

If you have a computational cluster, you can further parallelize work by splitting SNPs on chromosomes with `bedextract`:

```
$ for chr in `bedextract --list-chr snps.bed`; do echo ${chr}; bedextract ${chr} snps.bed > snps.${chr}.bed; done

```

Then submit a batch job for each chromosome:

```
$ bedmap --chrom ${chr} --echo --echo-map-id --delim '\t' myRegions.bed snps.${chr}.bed > answer.${chr}.bed

```

Using `bedextract` and adding `--chrom ${chr}` to your `bedmap` call will yield large speed improvements over other toolkits, as these options focus work on the chromosome of interest and nothing else.

# Answer
> 3 votes
I was tinkering with the command and was able to complete the execution inverting the order of rsid and bed intervals filtering. The command is as follows:

```
gzcat <myvcf>.vcf.gz | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]" | \
    java -Xmx10g -jar snpSift.jar intervals <mybed>.bed

```

Maybe my rsid subset is smaller than my BED intervals and thus the process is quicker and more efficient.

Anyway thanks everyone

# Answer
> 1 votes
Here is another way to do it in 2019:

```
$ bcftools view -T regions.bed -i 'ID=@<myrsid>.txt' input.vcf > output.vcf

```

Things can be sped up if we compress and index the `input.vcf`. Then `bcftools` will have near-immediate random access to the position.

```
$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz
$ bcftools view -T regions.bed -i 'ID=@<myrsid>.txt' input.vcf.gz > output.vcf

```

Internally `bcftools` works with `bcf`. Whenever the input and output is `vcf` it converts `vcf`\>`bcf`\>`vcf`. When dealing with large files, this produces a large overhead. Converting `vcf` to `bcf` first can reduce runtime for further downstream analyses dramatically. In my experience 2-3x faster.

```
$ bcftools view -Ob input.vcf > input.bcf
$ bcftools index input.bcf
$ bcftools view -Ob -T regions.bed -i 'ID=@<myrsid>.txt' input.bcf > output.bcf

```

---
Tags: vcf, bed, filtering, interval
---


# Question
Title: Seeking explanation of the hg38 files downloaded from bowtie 2 website
I downloaded the H. sapiens, NCBI GRCh38 files from Bowtie's website. After unzipping, there are 6 files, 4 that end in set 1.ebwt, set2.ebwt, set3.ebwt, and set4.ebwt and two that end in set.rev.1.ebwt and set.rev.2.ebwt.

I am unable to find any explanation for what these files are either online.

which one do I use for bowtie2 mapping?

Thanks

# Answer
> 3 votes
Bowtie is providing you with the index files for hg38. This are the result of the bowtie2-build indexer. However you have downloaded the indexes for Bowtie1 not bowtie2

> Bowtie 2’s .bt2 index format is different from Bowtie 1’s .ebwt format, and they are not compatible with each other.

Then according to the manual you pass the basename of the index for the reference genome. So this should be the name prior to the extension, .1.bt2 for Bowtie2.

```
bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r> | --interleaved <i>} -S [<sam>]
Main arguments
-x <bt2-idx>
The basename of the index for the reference genome. The basename is the name of any of the index files up to but not including the final .1.bt2 / .rev.1.bt2 / etc. bowtie2 looks for the specified index first in the current directory, then in the directory specified in the BOWTIE2_INDEXES environment variable.

```

---
Tags: sequence-alignment, genome, human, bowtie2
---


# Question
Title: Interpreting bayesian phylogenetic model using Tracer?
I am doing phylodynamic analysis of virus sequences using bayesian method using BEAUti and BEAST. As a part of my analysis I am using Tracer to check how my model is working.

I can interpret some of the results like from the Tracer like :

**clock.rate - rate of nucleotide substitution**

**treemodel.rootheight- time when divergence occured**

**ESS- measures sampling quality (ideally has to be greater that 200)**

But there are other many terms on this analysis which I have no idea how to explain.

Please can anybody help me to find out what other aspects of the evolution can be explained using this table.

posterior -33531.076
prior -481.982
likelihood -33049.094
treeModel.rootHeight 32.271
tmrca(Brazil) 3.853
tmrca(Canada) 16.644
tmrca(China ) 11.752
tmrca(Columbia) 1
tmrca(Thailand) 1.093
tmrca(USA) 32.271
tmrca(allseq) 32.271
constant.popSize 43.595
CP1.kappa 18.716
CP2.kappa 17.16
CP3.kappa 22.796
frequencies1 0.224
frequencies3 0.221
frequencies4 0.26 5001
CP1.mu 0.76 CP2.mu 1.221
CP3.mu 1.019
clock.rate 1.354E-3
meanRate 1.354E-3
CP1.treeLikelihood -9462.258
CP2.treeLikelihood -12520.305
branchRates 0E0
coalescent -467.376

# Answer
> 1 votes
When I have a moment I will fill in the blanks

You must have an RNA virus (nothing else will mutate that fast). The basal group and root is the USA. It estimates the time to that basal event is 32.3 (likely years).

"frequencies" frequency of ACGT where the missing one will be 0.25.

Its a simple model of constant population size.

ESS is the most important parameter to ensure convergence of the Bayesian simulation .

Likelihood is the second most important parameter. I'll edit again shortly.

---

I'll fill in one more blank and summarise ...

**Kappa** ... This will be the transition transversion ratio in HKY. I suspect it also will be transition transversion ratio in Kimuras 3-parameter model, this is because ... the transversions are further split into two mutation groups, ie. purine to pyrimidine has 2 further combinations and you end up with 2 ratios. that makes sense, transition/transversion (kappa 1), transition/transversion 1 (kappa 2), transition/transversion2 (kappa 3) ....

* CP1.kappa 18.716
* CP2.kappa 17.16
* CP3.kappa 22.796

CP2 and CP3 are either side of CP1 .. which is what you would expect on my interpretation.

These are very high values by the way, essentially there are very few transversions in this data, so ... I dunno. You need to check the data, if you see lots (and lots) of A-G and T-C mutations, then thats what I believe it is saying.

Briefly posterior and prior are the fundamental method of a Bayesian calculations, the output is probability expressed as a log, so the more negative the higher the probability (the prior looks like a random tree, or a reasonable tree with really rubbish parameters, likely the former).

---

**Summary** Essentially I am going to call it a day here. The real issue is that the output summarise the entire modern theory in phylogenetics, so its a very long response required to address all of the parameters and I honestly don't think it helps your analysis. What I am clear about is:

*** You've omitted a gamma distribution (I think Tracer would summarise this) ***

So my advice is:

* ESS MUST be homogeneous between runs - big trouble if not
* gamma distribution MUST be run

Beast can be tricky getting gamma to converge via ESS, just a word of warning.

I can't see the full robustness, because its absent... However, from the output, I think its pretty clear, that you (well me) can more or less see the model you ran.

---
Tags: sequence-alignment
---


# Question
Title: Running Snakemake in one single conda env
I have been experimenting a lot lately with Snakemake, I love it. Recently I also switched to using conda (--use-conda) in the way that is advertised. However, I have some issues with it, mostly related to the way we work. I work in a lab that processes data from our sequencer but also external datasets, moreover, multiple people use my pipeline.

So, with respect to conda this means: currently I have 5 envs, averaging 1.5 GB per env, this takes about 20 min to set up and it's quite some data to pull down. Especially since it is pulled down every time someone analyses data (sometime multiple times a day). This feels like quite a waste of resources. Also, I fix some package versions in the yaml files, but certainly not all (this may lead to issues when repos get reshuffled I read), so I also can't guarantee the exact environment long term (say over years, dependencies may change/update).

So I'm considering switching to a central conda env which I version and of which I export lists of all packages installed for reproducibility reasons. This will then contain all tools.

My question is, is there a nice way to do this? I had issues with putting a "source activate" line in every shell command (and how to do this for Python code under "run" anyway) before, see this: https://stackoverflow.com/questions/54024862/snakemake-conda-env-parameter-is-not-taken-from-config-yaml-file.

So, is there a way, also when using a qsub cluster, to always use a specific conda env, for everything you do?

# Answer
Using only the `--use-conda` directive, snakemake copies environments somewhere in the local `.snakemake` directory. This is repeated for each directory you run it in and replicates all environments effectively for each dataset.

You can however tell snakemake to use a shared environment location using `--conda-prefix` as described here:

> Specify a directory in which the ‘conda’ and ‘conda-archive’ directories are created. These are used to store conda environments and their archives, respectively. If not supplied, the value is set to the ‘.snakemake’ directory relative to the invocation directory. If supplied, the –use-conda flag must also be set. The value may be given as a relative path, which will be extrapolated to the invocation directory, or as an absolute path.

This way you could only download/copy the environment once and share it, e.g. this way:

```
snakemake --use-conda --conda-prefix /path/to/shared/environment

```
> 5 votes

# Answer
A good practice is to provide `yml` file along with your `Snakefile`. This file should contain all the enviroment definitions that are needed to run your pipeline. The minimum content of this file would be to install `snakemake` itself. Any other program used in your pipeline should be included as well. For reproducibility version numbers must be used.

So, if someone would like to use your pipeline, (s)he first creates a new enviroment based on your definition:

```
$ conda env create -f my_pipeline.yml

```

Afterwards activate this environment and run the pipeline:

```
$ conda activate my_pipeline
$ snakemake

```

The `conda` directive in your `Snakefile` is only neccessary if you need to use program which dependencies cannot be resolved within the main environment. If you start snakemake with the `--use-conda` parameter, snakemake will create a new environment based on the provided definition file and activate it for the rule.

> Especially since it is pulled down every time someone analyses data (sometime multiple times a day).

The environment creation is only done once. If you run the pipeline a second time `snakemake` can activate the environment diretcly without downloading and installing something again.

The tutorial I wrote some time ago on biostars about `snakemake` and `conda` might be useful for you.
> 5 votes

# Answer
Inside your Snakefile you can provide the same `conda:` directive to each rule. If you install your pipeline in a central location and have everyone use the same conda when running it (hint: use a wrapper script that sets the path to conda) then everyone will be using the same environment. This has the added benefit that you can make this readable to everyone, but disallow them to make any modifications (so users can't mess things up).

The tack we've taken is to have a basic conda env with snakemake and little else and then separate envs for various pipelines (and versions thereof if software versions in your YAML files change) or parts thereof, in cases where you have rules using tools with conflicting dependencies. This ends up being fairly convenient for more complicated pipelines, though if you can get away with just a single environment then finswimmer's answer is simpler.

Note that all of this should work with any scheduler, since snakemake will insert the environment activation stuff in the script it submits to your scheduler.

As an aside, this is exactly how snakePipes works.
> 1 votes

---
Tags: python, snakemake, bioconda
---


title: Building RNAz on MacOS

text:
# Question
Title: Building RNAz on MacOS
Having trouble building RNAz 2.1 on MacOS. Following instructions at https://www.tbi.univie.ac.at/software/RNAz/

```
tar -xzf RNAz-2.1-tar.gz
cd RNAz-2.1
./configure
make

```

`make` fails with

```
/usr/bin/clang++ -o RNAz RNAz.o rnaz_utils.o svm_helper.o zscore.o cmdline.o strand.o svm.o -L../librna -lRNA -lm
Undefined symbols for architecture x86_64:
  "_HairpinE", referenced from:
      _alifold in libRNA.a(alifold.o)
      _fill_arrays in libRNA.a(fold.o)
      _backtrack in libRNA.a(fold.o)
      _stack_energy in libRNA.a(fold.o)
  "_LoopEnergy", referenced from:
      _alifold in libRNA.a(alifold.o)
      _fill_arrays in libRNA.a(fold.o)
      _backtrack in libRNA.a(fold.o)
      _stack_energy in libRNA.a(fold.o)
ld: symbol(s) not found for architecture x86_64

```

Thinking that compiling for 32-bit might work I tried setting

```
CXXFLAGS = -m32
LDFLAGS = -m32

```

without success (same result). Same result with `gmake`.

**Configure** Configure checks for `gcc`

```
...
checking for gcc... gcc
checking whether the C compiler works... yes
...

```

**Makefile**

Makefile specifies the following:

```
CC = gcc
CCDEPMODE = depmode=gcc3
CFLAGS = -g -O2
CPP = gcc -E
CPPFLAGS =
CXX = g++
CXXDEPMODE = depmode=gcc3
CXXFLAGS = -g -O2

```

**Compiler**

Apple LLVM version 10.0.0 (clang-1000.11.45.5)

Target: x86\_64-apple-darwin18.2.0

Thread model: posix

InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin

Has anyone built this successfully on MacOS? Or perhaps suggestions as to how I might build this successfully? I'm at a loss.

# Answer
> 3 votes
Compile with gcc or remove `inline` from `librna/fold.c:776` and `librna/fold.c:846`.

# Answer
> 2 votes
The other option would be to use RNAz via a docker container instead on installing it directly. The docker container can be downloaded from biocontainers. If you are new to docker here is a tutorial for mac.

# Answer
> 1 votes
The first thing to try is to perform the install as sudo.

```
sudo ./configure # definitely
sudo make ... sudo make install # possibly (?)

```

The clang++ compiler may be using a path that is not setup in your profile and sudo will overcome this. If sudo works you know what the problem is and its a case of changing permissions and configurations accordingly. There is a good chance this will resolve the matter.

There are other possibilities, it could be the compiler. It looks like the compiler is fine however looking at the config file and the check. The final issue is it could be the program ;-)

> Tried this. Same result. Will continue scratching head

Try an installation on Linux

Otherwise its usual stuff, namely upgrading Xtools

---
Tags: rna
---


# Question
Title: How to find/build the evolutionary history of a protein from its sequence?
I'd like to build the evolutionary history of a protein, given its sequence. Namely, given a FASTA entry how can I build an evolutionary tree? Here is the 5wxy protein as an example:

```
>5WXY:A|PDBID|CHAIN|SEQUENCE
MGHHHHHHMMKTKLPILGVLGGMGPVVTAEFLKSIYEYNPFIDKEQESPNVIVFSFPSAPDRTGSIDSGKEREFIDFIQV
NLEHLNKLADCIVIGSCTAHYALPQIPENLKDKLISLIKIADQELQEYNEPTLLLASTGTYQKKLFQEGCTTADLIISLS
ESDQKLIHEMIYKVLKRGHDPLSILRDIEALLEKYNTRSYISGSTEFHLLTKSLKLKGIDSIKAIDPLSTIAQNFSQLII
KQAQVDLVTDCHQPSNPKSP

```

Do I need to do some queries to Entrez in order to build the evolution tree?

I saw this paper, but I am not sure how I can use it with a single mammal's protein primary sequence.

---

**Question** given a protein name, is it possible to look for its evolutionary history or build evolution tree using other databases beside PDB on top of bioinformatics algorithms.

# Answer
> 3 votes
The general procedure is:

1. Find the sequence for the same (or as similar as you can find (see mentions of blast in the comments above)) gene in other species. You can use uniprot or any other large sequence database for this.
2. You now have a large set of sequences related to yours, which you'll need to compare. Search the methods section of your favorite paper to see how they reconstructed likely evolutionary history from this, there are a number of packages out there.

# Answer
> 2 votes
There you go, its Microcystis aeruginosa aspartate/glutamate racemase. Easy does it. If you want to know how you can do this in 15 seconds let me know.

---

Its one of the new features of NCBI's Blast

1. Go to Blast here, https://blast.ncbi.nlm.nih.gov/Blast.cgi
2. Enter your sequence into the box (it doesn't accept PDB codes alone)
3. Enter the protein database - when I first did this calculation I used SwissProt, thinking there would be alot of sequences - I then used "nr"
4. Under the algorithm parameters enter "50" (default is too many)
5. Hit "Blast"
6. Once the search is complete at the top of the page are the hyperlinks: "Other reports: Search Summary \[Taxonomy reports\] \[Distance tree of results\]"
7. Click on "Distance tree of results"
8. The following page will load automatically, automatically aligning your sequences and producing, in this case a parsimony based tree, but there is also the option of a nj tree (recommended) .... here: https://www.ncbi.nlm.nih.gov/blast/treeview/treeView.cgi?request=page&blastRID=5A1ZBRU1014&queryID=lcl|Query\_210080&entrezLim=&ex=&exl=&exh=&ns=50&screenWidth=1280&screenHeight=800
9. Click "Tool", "Download", "PDF" ...

---
Tags: phylogenetics, sequence-alignment, proteins, phylogeny, pdb
---


# Question
Title: Gene Ranking - signal to noise ratio used in GSEA-P algorithm?
I'm looking at Broad Institute's orignal GSEA-P algorithm R script which I downloaded here: http://software.broadinstitute.org/gsea/downloads.jsp.

I'm trying to adapt their GSEA.1.0.R script to process datasets that have 1 gene expression profile as for class 1, and N gene expression profiles for class 2. The bottleneck seems to be the GSEA.MetaboliteRanking() method. I looked and it appears they are using "signal to noise ratio" (SNR) metric to rank genes in patient profiles.

**I was wondering how the signal to noise is calculated from gene expression profiles?**

Note: I am using/modifying their script because it's the only implementation of the GSEA algorithm that I could find where they expose their code.

i.e., what is the mathematical reasoning behind the following lines of code? Note, there are little to no comments in the code, but the following variables, from what I can tell, are:

```
A: Matrix of gene expression values (rows are genes, columns are samples)
class.labels: Phenotype of class distinction of interest. A vector of
              binary labels having first the 1s and then the 0s

P = class.labels1 * subset.mask
n1 = sum(P[,1])
M1 = A %*% P
M1 = M1/n1
A2 = A*A
S1 = A2 %*% P
S1 = S1/n1 - M1*M1
S1 = sqrt(abs((n1/(n1-1)) * S1))

P = class.labels2 * subset.mask
n2 = sum(P[,1])
M2 = A %*% P
M2 = M2/n2
A2 = A*A
S2 = A2 %*% P
S2 = S2/n2 - M2*M2
S2 = sqrt(abs((n2/(n2-1)) * S2))

```

# Answer
> 1 votes
The noise is comming from the subset.mask, which is created above in a loop with the number of permutations.

```
for (r in 1:nperm) { #L90
....
        subset.mask[, r] <- as.numeric(c(subset.class1, subset.class2)) # L107

```

So by multiplying the random selection of subsets by the expression we get the "noise". Later on line 251 we get the ratio of noise/signal by dividing the signal from the contrast and the signal by the noise.

Those transformation on the middle seem to come from normalizing the signal in order to be able to compare the resampling, but I don't fully understand why they do this. Also note that the R scripts are outdated and no longer maintained, so they can contain errors.

---
Tags: r, gene-expression, gsea
---


title: Python string editing in Snakemake

text:
# Question
Title: Python string editing in Snakemake
I'm running a custom perl script using Snakemake. I use this rule:

```
rule complexity_20mer_counter:
    input:
       os.path.join(fastq_trimmed_dir, '{sample}_{read}' + fastq_extension.split('.')[0] + '_val_' + '{read}'.strip('R') + '.fq.gz')
    output:
        os.path.join(custom_qc_dir, '{sample}_{read}_20mer_counter.txt')
    shell:
        '''
        #!/bin/bash
        perl scripts/20mer_counter.pl {input} > {output}
        '''

```

It fails (I'll get to why later) with this error:

```
Building DAG of jobs...
MissingInputException in line 1 of /rst1/2017-0205_illuminaseq/scratch/swo-406/2019-01-24/rules/20mer_counter.smk:
Missing input files for rule complexity_20mer_counter:
fastq_trimmed/xs0067_P2018SEQE27R01_S1_R1_001_val_R1.fq.gz

```

... this is because I use the output of trim\_galore, this rule looks like this:

```
rule trim_galore:
    input:
       [os.path.join(fastq_raw_dir, '{sample}_'+ read + fastq_extension)
       for read in reads]
    output:
        timmed_fastq_files = [os.path.join(fastq_trimmed_dir, '{sample}_' + read + fastq_extension.split('.')[0] + '_val_' + read.strip('R') + '.fq.gz')
        for read in reads],
        trimming_reports = [os.path.join(trim_galore_qc_dir, '{sample}_'+ read + fastq_extension + '_trimming_report.txt')
        for read in reads]
    conda:
        "../envs/fastqc.yaml"
    shell:
        '''
        #!/bin/bash
        trim_galore --fastqc --gzip -o {fastq_trimmed_dir} --paired --retain_unpaired {input}
        # Move all qc reports to the fastq_trimmed_dir
        mv {fastq_trimmed_dir}/*.zip \
           {fastq_trimmed_dir}/*.html \
           {fastq_trimmed_dir}/*.txt \
           {trim_galore_qc_dir}
        '''

```

As you can see, trim\_galore not only uses my {read} parameter in the output filename but again adds it at the end, only now without the "R".

I tried solving it in my first rule by applying a python .strip('R') on my {read} wildcard, but this does not seem to work. How can I correctly get '1' and '2' from 'R1' and 'R2' (which are my {read} wildcards)?

(Extra bonus points for a nice output.trim\_galore.something as input for my complexity\_20mer\_counter rule, which does not work because the output is a list (it needs 2 fastq files at the same time), but my complexity\_20mer\_counter rule *can* be run in parallel again.)

# Answer
It seems that it is not possible to apply `strip()` directly on the wildcard placeholders. Luckily it is possible to executed functions in the `input` directive.

So try changing your input of the `complexity_20mer_counter` rule to this:

```
lambda wc: os.path.join(fastq_trimmed_dir, wc.sample + '_' + wc.read + fastq_extension.split('.')[0] + '_val_' + wc.read.strip('R') + '.fq.gz')

```

The wildcard object is passed to our anonymous function. So have access to the wildcards value via `wc.sample` and `wc.read`, which are now handled as normal strings and we can apply `strip` on them.
> 2 votes

---
Tags: python, snakemake
---


# Question
Title: How can I read FCS files using open source libraries?
FCS is a patented data format used for storing flow cytometry data. The most recent version is FCS3.1. There is some documentation on the format, but there is no information on how to read these files. There are some R packages and a MATLAB code to read an FCS file, but I am looking for standard libraries developed either by the FCS consortium or any other group. I also wish to know if FCS is a subset of an existing standard data format that can be read by a standard library using any programming language.

Finally, I would want to convert these files to an easily readable format like HDF5.

# Answer
> 5 votes
A few years ago, I wrote a python script to convert FCS files into tab-separated format. It was far from handling all the possibilities that the format description offers, but at least it worked for some of the files produced on one of our machine: http://www.igh.cnrs.fr/equip/Seitz/en\_equipe-programmes.html

The format documentation I found enabled decoding (see section 3 of the pdf you mention), but it requires reading data in binary mode.

The general idea of this format (and, I guess, many other binary formats), is that there is a header zone at the beginning of the file with a defined number of fields encoding numbers indicating how the rest of the file is structured. So a first phase is to parse this header, following the description given in the documentation of the format. The information extracted from the header tells where to find the data and how it is encoded, still according to rules described in the format documentation.

In case this may be useful, and for the record, here is the code from the above-mentioned script (after stripping comments, some of which are merely copied from the format documentation, and adding a few ones):

```
#!/usr/bin/env python
"""This script tries to read FCS flow cytometry data.
Format parsing inspired by information found here:
http://isac-net.org/Resources-for-Cytometrists/Data-Standards/Data-File-Standards/Flow-Cytometry-Data-File-Format-Standards.aspx
"""

import re
# To decode binary-encoded data
import struct
import sys

class Parameter(object):
    """This object represents one of the parameter types that are present in a DATA segment of a FCS file."""
    __slots__ = ("p_name", "p_bits", "p_range", "p_ampl", "parser")
    def __init__(self, p_name, p_bits, p_range, p_ampl):
        self.p_name = p_name
        self.p_bits = p_bits
        self.p_range = p_range
        self.p_ampl = p_ampl
        # Function for parsing a value of the parameter in the data segment
        self.parser = None

##############################################
# Here starts the parsing of the header part #
# which tells where the other parts are.     #
##############################################

f = open(sys.argv[1], "rb")
# The format name is encoded in 6 letters
# An ASCII letter is coded with one octet
file_format = "".join([f.read(1) for __ in range(6)])
sys.stdout.write("Format: %s\n" % file_format)
# The format descriptions reserves 4 octets that we skip
skip = f.read(4)
# 8 octet chunks encode the start and end positions
# of different parts of the data
text_start = int(f.read(8).strip(" "))
text_end = int(f.read(8).strip(" "))
data_start = int(f.read(8).strip(" "))
data_end = int(f.read(8).strip(" "))
analysis_start = int(f.read(8).strip(" "))
analysis_end = int(f.read(8).strip(" "))

if (analysis_start and analysis_end):
    sys.stderr.write("Cannot deal with ANALYSIS segment of an FCS file.\n")

####################################################
# Here starts the parsing of the "TEXT" portion    #
# which describes how the data proper is organized #
####################################################
f.seek(text_start)
# The first character in the primary TEXT segment is the ASCII delimiter character.
sep = f.read(1)
if sep not in ["_", "@"]:
    alt_sep = "_@_"
elif sep not in ["_", "|"]:
    alt_sep = "_|_"
else:
    assert sep not in ["+", "|"]
    alt_sep = "+|+"
text_segment = f.read(text_end - text_start)

fields = text_segment.split(sep)

info = {}

i = 0
while i < len(fields) - 1:
    key = fields[i]
    i += 1
    val = fields[i]
    i += 1
    # Keywords are case insensitive, they may be written in a file in lower case, upper case, or a
    # mixture of the two. However, an FCS file reader must ignore keyword case. A keyword value may
    # be in lower case, upper case or a mixture of the two. Keyword values are case sensitive.
    info[key.upper()] = val
print "%s events were detected." % info["$TOT"]
print "Each event is characterized by %s parameters" % info["$PAR"]

if info["$NEXTDATA"] != "0":
    sys.stderr.write("Some other data exist in the file but hasn't been parsed.\n")

# L - List mode. For each event, the value of each parameter is stored in the order in which the
# parameters are described. The number of bits reserved for parameter 1 is described using the
# $P1B keyword. There can be only one set of list mode data per data set. The $DATATYPE
# keyword describes the data format. This is the most versatile mode for the storage of flow
# cytometry data because mode C and mode U data can be created from mode L data.
assert info["$MODE"] == "L"

parameters = []

# indices of the parameters
p_indices = range(1, int(info["$PAR"]) + 1)
for i in p_indices:
    p_name = info["$P%dN" % i]
    p_bits =  info["$P%dB" % i]
    p_range = info["$P%dR" % i]
    p_ampl =  info["$P%dE" % i]
    parameters.append(Parameter(p_name, p_bits, p_range, p_ampl))

sys.stdout.write("The parameters are:\n%s\n" % "\t".join([par.p_name for par in parameters]))

# How are 32 bit words organized
if info["$BYTEORD"] == "4,3,2,1":
    endianness = ">"
else:
    endianness = "<"
    assert info["$BYTEORD"] == "1,2,3,4"

    # I stripped a long comment which is just a copy of the documentation
# Type of data:
if info["$DATATYPE"] == "I":
    for par in parameters:
        nb_bits = int(par.p_bits)
        assert nb_bits % 8 == 0
        nb_bytes = nb_bits / 8
        # Determine format string for unpacking (see https://docs.python.org/2/library/struct.html)
        if nb_bytes == 1:
            c_type = "B" # unsigned char
        elif nb_bytes == 2:
            c_type = "H" # unsigned short
        elif nb_bytes == 4:
            c_type = "L" # unsigned long
        elif nb_bytes == 8:
            c_type = "Q" # unsigned long long
        else:
            raise ValueError, "Number of bytes (%d) not valid for an integer (see https://docs.python.org/2/library/struct.html#byte-order-size-and-alignment)." % nb_bytes
        fmt = "%s%s" % (endianness, c_type)
        p_range = int(par.p_range)
        def parser(data):
            value = struct.unpack(fmt, data.read(nb_bytes))[0]
            try:
                assert value < p_range
            except AssertionError:
                print "Value %s higher than %d" % (str(value), p_range)
            return value
        par.parser = parser
    pass
else:
    raise NotImplementedError, "Only the parsing of integer value has been implemented so far."

out_file = open(sys.argv[2], "w")
out_file.write("#amplification_types\t" + "\t".join([par.p_ampl for par in parameters]) + "\n")
out_file.write("parameters\t" + "\t".join([par.p_name for par in parameters]) + "\n")
i = 1
##############################################
# Here starts the parsing of the data proper #
##############################################
f.seek(data_start)
while f.tell() < data_end:
    values = []
    for par in parameters:
        values.append(par.parser(f))
    out_file.write("%d\t" % i + "\t".join(map(str, values)) + "\n")
    i += 1
out_file.close()
f.close()

```

# Answer
> 4 votes
I maintain several Python libraries for interacting with flow cytometry data (FCS files). The most recent of which, FlowKit, provides functionality to export the event data as a NumPy array, Pandas DataFrame, or CSV text. Using the `h5py` library, it's trivial to export as an HDF5 file

# Answer
> 3 votes
R/Bioconductor has a number of different flow cytometry processing packages. One place to start for looking at cytometry data from a high level would be openCyto (or its vignette), which is a large set of tools for basic extraction and analysis of FCS files.

I have looked in the past at the FCS files as an R structure using flowCore. Loading a single FCS file is fairly straightforward and follows a familiar R pattern:

```
file.name <- "/dir/file.fcs"
x <- read.FCS(file.name, transformation=FALSE)
summary(x)

```

Asking about converting FCS files into an "easily readable format like HDF5" doesn't seem like the right question. HDF5 is a container format and shares a lot of similarity to file systems. I've found it best to keep FCS files as they are, as it is a compact, standardised binary format.

---
Tags: file-formats
---


# Question
Title: Is there a way to import tables in PubMed Central to MySQL or other Relational databases?
I need to use these PMC tables from the medical journals for the purpose of integrating the data with a decision-making software. Though, I can see that one can access full-text for an article from PMC here (which includes the article table), but is there a way to extract the table from this XML document and export it to a local MySQL database or any other relational database?

# Answer
> 1 votes
If you can use R, consider using web scraping.

See here : http://bradleyboehmke.github.io/2015/12/scraping-html-tables.html

---
Tags: public-databases, text-processing
---


# Question
Title: How to score how densely bases are distributed along a given oligo?
Suppose I have a short chain string of oligos, and I want to rate them from 0 to 1 based on how clustered the "A"s are in respect to the chain. For example:

```
"AAAAACCCCC" #should score as 1 for being as clustered as it possibly could
"ACACACACAC" #should be scored as 0 for being as thinly distributed as it could

```

Current idea is to strip both ends of A's, then remove first and last characters, then

```
gaps = ''.join(str(x) for x in deck).split("A")
listofCs = [x for x in gaps if x]

```

to join it all together as a list of how big the gaps are (i.e. the "C"s), and simply count the average size

```
average = sum(map(len,listofCs))/len(listofCs)

```

But I feel like there should be a much more elegant solution to this. It feels messy. How can I improve the scoring?

# Answer
You could count the number of transitions between an A and a 'not A', and divide this number by the total number of As. You might also need to scale by the sequence length, if you have different length sequences that you want to compare.

Edit: also, there might be a correction to make if one or both ends of the sequence are an A
> 2 votes

# Answer
To define the base density score, or clustering score, of a sequence $S$ we can think of each position in the sequence as being part of a Markov process. For the sake of symmetry around a single base we consider odd-order Markov models from 3-onward where we ask what is the probability that we see the character at $S\_i$ in the $j$ bases from $\[S\_{i-j}, S\_{i+j}\]$. For each column we have `length(S) -1`, $|S|-1$, possible probabilities to calculate per window.

We can't look left or right of the edges of the sequences so the probabilities are skewed, but we can get a rough idea of what is going on.

* For each base position we calculate a clustering score of $ \frac{1}{w}\Sigma\_{j=1}^{j=w} -log\_2 \frac{m(S\_{i-j}, S\_{i+j})-1 + \frac{1}{|S|}}{|S\_{\\{i-j, i+j\\}}| -1 }$ where:
    + $w$ is the window size we can look at. The max is $|S|-1$ and the min is $1$. This just means that we can calculate a score in the context of the whole sequence or we maybe only care within some smaller window, like 5 bp flanking either side of the base in question.
    + The function $m(S\_{i-j}, S\_{i+j})$ is just the number of counts of matching bases in the window $j$ on either side of the base $S\_i$. For example in the sequence `gatAgga`, if $S\_i$ is the uppercase `A` and $j$ is 3, the function $m(S\_{i-3}, S\_{i+3})$ returns 3, since there is the `A` base at $S\_i$ and the two additional `a` bases at $i-2$ and at $i+3$. We have to subtract one from this number because we don't want to factor in the matching base at $S\_i$.
    + We then add $\frac{1}{|S|}$ to that number. This is a pseudocount that we need to add to prevent from taking the log of zero. This number is pretty arbitrary but I chose $\frac{1}{|S|}$ so that non-clusting regions' scores aren't too high.
    + We then divide that numerator by the length of $S$ from the little window that we just looked at minus one, $|S\_{\\{i-j, i+j\\}}| -1 $. For example in the sequence `gatAgga` from above we only want to consider the bases `gat gga`. This is 6 bases.
    + We then take the negative log to convert the proportional scale to a log scale.
    + The summation just says that we do this for every symmetrical window size around the base in question then sum them up. For example we would calculate the ratio for `tAg`, `atAgg`, and `gatAgga`.
    + The $\frac{1}{w}$ just gives the average log clustering score of all the windows that we looked at. One nice thing about this is that the impact of each increasing window size has less effect on the ratio, so there is a sort of convergence that happens.

Implemented in python:

```
from math import log as mlog
s = "GGATAGGACATTTTTTTTTTGAGCTTTTATTTTCCGA"
assert len(s) > 0
pseudocount = 1/len(s)
#maxdepth = 5
maxdepth = len(s)
score = []
for i in range(len(s)): # iterate positions
    temp = []
    for j in range(1, maxdepth): # iterate depth
        rmin = max(i-j, 0)
        rmax = i+j+1
        num_matches_min_one = sum([1 for k in s[rmin:rmax] if k == s[i]])-1
        len_s_min_one = len(s[rmin:rmax]) -1
        sim = abs(-1 * mlog((num_matches_min_one/len_s_min_one) + pseudocount, 2))
        #print("i={}, j={}, rmin={}, rmax={}, sim={}".format(i, j, rmin, rmax, sim))
        temp.append(sim)
    print("{0:d}\t{1}\t".format(i, s[i]), end='\t')
    score.append(sum(temp)/len(temp))
    print("{0:.4f}".format(score[-1]))

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
x = [i+1 for i in range(len(s))]
plt.bar(x, score)
plt.xticks(x, [k for k in s])
ax.set_yscale("log", nonposy='clip')
plt.xlabel('sequence')
plt.ylabel('-log2(cluster-ness)')
plt.title("cluster-ness score")
plt.savefig('output.png', dpi=150)

```

This calculates the matrix, prints out a tab-delimited file of the position, base, and log clusting score, then saves a plot of the results.

For the test sequence `GGATAGGACATTTTTTTTTTGAGCTTTTATTTTCCGA` and using the whole sequence as context, I got the following plot:

As you can see, the log clustering score is low around the things that are cluster-y and high around things that aren't cluster-y.

**Note:** OP requested that the score be scaled from 0 to 1. I instead used the log scale to make visualization easier.

**Warning 1:** If you plan on doing this for large sequences, it is really important to limit the `maxdepth`, or $w$, parameter. The time complexity for this algorithm is $O(|S|\*w)$, so the number of calculations can blow up to $O(|S|^2)$.

**Warning 2:** It is also important to limit the $w$ parameter because the edge effects propagate through the entire sequence when $w=|S|$. Look at these examples built off of OP's initial suggestions of `AAAAACCCCC` and `ACACACACAC`. The clustering score never stabilizes when $w=|S|$, however it does stabilize away from the edges when $w=5$.

For 50 characters of `ACACAC...`:

For 50 characters of `[25xA][25xC]`:
> 0 votes

# Answer
The Markov chain approach is very interesting. You could deploy a maximum likelihood approach where instead of scoring the genetic distance within a homologue, the algorithm looks at its physically linked neighbour. The genetic distance of the transition point i.e. AAAA then CCCC is deducted from the final genetic distance so that AAAAACCCCCC would be zero. It would be fairly easy to script and use PAML or equivalent where the default phylogeny would be a star tree.

The challenge is to establish a Deep Learning approach. Its not impossible, you just assign pre-determined classes of oligos which fall into the definition you want and use these as your training set.
> 0 votes

---
Tags: ngs, algorithms
---


title: Perfect Phylogeny vs Maximum parsimony

text:
# Question
Title: Perfect Phylogeny vs Maximum parsimony
I am searching various sources about phylogenetics. I saw some materials about perfect phylogeny and also phylogenies acquired from maximum parsimony constraint. They seem very similar to me. Are they the same?

# Answer
> 3 votes
No. A perfect phylogeny is one such that all characters evolve on the tree with no homoplasy, i.e. for a binary character, changes occur once from 0 -\> 1 but never from 1 -\> 0. Maximum parsimony inference *may* produce a perfect phylogeny, but typically for real datasets some degree of homoplasy is required to explain character patterns.

# Answer
> 2 votes
Long discussion that took many years to resolve.

MP doesn't account for back mutation, which is a HUGE problem for nucleotide data because in theory 1:4 mutations is a back mutation. Maximum likelihood has resumed is crown here, however Beast deployment of a Bayesian calculation is also very widely, particularly for molecular dating. ML uses a reversible matrix A \<-\> T, whilst Beast uses a directional matrix A-\>T and will score T-\>A separately.

Happy to discuss the approaches and modern interpretations at length.

The model which accounts for back-mutation is the Jukes-Cantor correction (JC correction), this is a basic model present in every phylogenetic algorithm except p-distances). Essentially JC extrapolates back-calculation for the true number of mutations against the observed number of mutations against divergence time. However, when the P-distance (uncorrected observed distance) exceeds 0.75 for nucleotides the model is not viable. Basically, nucleotide divergence \>0.75 is saturated and the phylogenetic information is essentially random.

---
Tags: phylogenetics, ngs, phylogeny
---


# Question
Title: Hardware Requirements (specs) for Bioinformatics-dedicated desktop
I know this is a somewhat general or vague question, but I’m interested in your opinions. I must build a desktop for general bioinformatics activities with human genomes. I will work with Python and R libraries, I will also need to install denovo sequence assembly software, alignment tools, variant viewers/browsers, sql, among other software. I don’t have much experience customizing hardware that will be dedicated to bioinformatics and genome data analyses. I’m pretty sure I work with linux and/or iOS. What is the recommended RAM, processor speed, HD capacity, graphics cards, etc. I appreciate your opinions on the matter. Thank you.

# Answer
> 4 votes
The RAM requirements for some bioinformatics analyses like assembly can be quite high (in the hundreds of Gigabytes). My recommendation is to get a fast laptop. Something with an i7 quad-core processor, 16 GB of RAM, and 1TB of storage should do. Do your analyses on this laptop. If you run out of RAM or storage, then use a cloud service to get access to bigger machines or buy an external hard drive. Once you have a feeling for what your actual needs are, start looking at a dedicated rack server solution in your local institution's server farm or again, in the cloud.

You don't mention any analyses that clearly benefit from a GPU, but as you get more familiar with your analyses, perhaps that will become useful as well.

# Answer
> 2 votes
If you deal with genomics NGS data then heavy RAM will be essential. I don't think you will get this from a reasonable priced laptop.

If you are using a Deep Learning based algorhithms or Bayesian or maximum likelihood heavy duty processors are needed.

I would go for an allround machine. Two hard drives:

* 0.5 to 1 Tb SSD (I use 1TB),
* 2 to 4 Tb Stata data storage,
* 4-core (light) or 8-core, i7 (10-core is common now),
* 16GB-32GB RAM.

**Important**: For your purpose if you are doing de novo NGS assemblies use:

* ALOT of RAM.
* Sufficient harddrive - make sure this is backed up

If you want evidence why RAM is important simply look at the threads here and on Stackexchange which read ... de novo NGS assembly PC has stopped working?(!).

The following threads give an insight into the sort of problems RAM creates

The following looking like a particularly bad case:

---
Tags: sequence-alignment, ngs, sequence-analysis
---


title: Error in installing Vep

text:
# Question
Title: Error in installing Vep
I am trying to installing Vep on Linux but an error is permanent

```
[fi1d18@cyan01 ensembl-vep-release-95]$ perl INSTALL.pl

Hello! This installer is configured to install v95 of the Ensembl API for use by the VEP.
It will not affect any existing installations of the Ensembl API that you may have.

It will also download and install cache files from Ensembl's FTP server.

Checking for installed versions of the Ensembl API...done

Setting up directories
Destination directory ./Bio already exists.
Do you want to overwrite it (if updating VEP this is probably OK) (y/n)? y
 - fetching BioPerl
 - unpacking ./Bio/tmp/release-1-6-924.zip
 - moving files
Attempting to install Bio::DB::HTS and htslib.

>>> If this fails, try re-running with --NO_HTSLIB

 - checking out HTSLib
fatal: destination path 'htslib' already exists and is not an empty directory.
 - building HTSLIB in ./htslib
In /home/fi1d18/ensembl-vep/ensembl-vep-release-95/htslib
gcc -g -Wall -O2 -fPIC -Wno-unused -Wno-unused-result -I.  -c -o cram/cram_io.o cram/cram_io.c
cram/cram_io.c:61:18: error: lzma.h: No such file or directory
cram/cram_io.c: In function 'lzma_mem_deflate':
cram/cram_io.c:679: warning: implicit declaration of function 'lzma_stream_buffer_bound'
cram/cram_io.c:685: error: 'LZMA_OK' undeclared (first use in this function)
cram/cram_io.c:685: error: (Each undeclared identifier is reported only once
cram/cram_io.c:685: error: for each function it appears in.)
cram/cram_io.c:685: warning: implicit declaration of function 'lzma_easy_buffer_encode'
cram/cram_io.c:685: error: 'LZMA_CHECK_CRC32' undeclared (first use in this function)
cram/cram_io.c: In function 'lzma_mem_inflate':
cram/cram_io.c:695: error: 'lzma_stream' undeclared (first use in this function)
cram/cram_io.c:695: error: expected ';' before 'strm'
cram/cram_io.c:701: error: 'LZMA_OK' undeclared (first use in this function)
cram/cram_io.c:701: warning: implicit declaration of function 'lzma_stream_decoder'
cram/cram_io.c:701: error: 'strm' undeclared (first use in this function)
cram/cram_io.c:701: warning: implicit declaration of function 'lzma_easy_decoder_memusage'
cram/cram_io.c:719: warning: implicit declaration of function 'lzma_code'
cram/cram_io.c:719: error: 'LZMA_RUN' undeclared (first use in this function)
cram/cram_io.c:720: error: 'LZMA_STREAM_END' undeclared (first use in this function)
cram/cram_io.c:732: error: 'LZMA_FINISH' undeclared (first use in this function)
cram/cram_io.c:743: warning: implicit declaration of function 'lzma_end'
At top level:
cc1: warning: unrecognized command line option "-Wno-unused-result"
make: *** [cram/cram_io.o] Error 1
Compile didn't complete. No libhts.a library file found at INSTALL.pl line 820.
[fi1d18@cyan01 ensembl-vep-release-95]$

```

And

```
[fi1d18@cyan01 ensembl-vep]$ perl INSTALL.pl --NO_HTSLIB

Hello! This installer is configured to install v95 of the Ensembl API for use by the VEP.
It will not affect any existing installations of the Ensembl API that you may have.

It will also download and install cache files from Ensembl's FTP server.

Checking for installed versions of the Ensembl API...done

Setting up directories
Destination directory ./Bio already exists.
Do you want to overwrite it (if updating VEP this is probably OK) (y/n)? y
 - fetching BioPerl
 - unpacking ./Bio/tmp/release-1-6-924.zip
 - moving files

Downloading required Ensembl API files
 - fetching ensembl
 - unpacking ./Bio/tmp/ensembl.zip
 - moving files
 - getting version information
 - fetching ensembl-variation
 - unpacking ./Bio/tmp/ensembl-variation.zip
 - moving files
 - getting version information
 - fetching ensembl-funcgen
 - unpacking ./Bio/tmp/ensembl-funcgen.zip
 - moving files
 - getting version information
 - fetching ensembl-io
 - unpacking ./Bio/tmp/ensembl-io.zip
 - moving files
 - getting version information

Testing VEP installation
ERROR: Testing VEP script failed with the following error
Can't locate Time/HiRes.pm in @INC (@INC contains: /home/fi1d18/ensembl-vep/modules /home/fi1d18/ensembl-vep ./Bio /usr/local/lib64/perl5 /usr/local/share/perl5 /usr/lib64/perl5/vendor_perl /usr/share/perl5/vendor_perl /usr/lib64/perl5 /usr/share/perl5 .) at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/Utils/SqlHelper.pm line 93.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/Utils/SqlHelper.pm line 93.
Compilation failed in require at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/DBSQL/DBConnection.pm line 81.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/DBSQL/DBConnection.pm line 81.
Compilation failed in require at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/DBSQL/DBAdaptor.pm line 69.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/DBSQL/DBAdaptor.pm line 69.
Compilation failed in require at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/Registry.pm line 137.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/Registry.pm line 137.
Compilation failed in require at /home/fi1d18/ensembl-vep/modules/Bio/EnsEMBL/VEP/BaseVEP.pm line 58.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/modules/Bio/EnsEMBL/VEP/BaseVEP.pm line 58.
Compilation failed in require at (eval 7) line 3.
        ...propagated at /usr/share/perl5/base.pm line 94.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/modules/Bio/EnsEMBL/VEP/BaseRunner.pm line 56.
Compilation failed in require at (eval 6) line 3.
        ...propagated at /usr/share/perl5/base.pm line 94.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/modules/Bio/EnsEMBL/VEP/Runner.pm line 71.
Compilation failed in require at ./vep line 20.
BEGIN failed--compilation aborted at ./vep line 20.

[fi1d18@cyan01 ensembl-vep]$

```

# Answer
What you can try is stated out in the message:

> If this fails, try re-running with --NO\_HTSLIB
> 1 votes

---
Tags: ensembl, linux, perl
---


title: wget for links inside html pages

text:
# Question
Title: wget for links inside html pages
I am trying to download a file from the following repository: https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR7276474

As you can see, there are several layers to the webpage. For example, clicking on the download tab doesn't change the URL and, the link is not a 'download link' per ce - where simply clicking the link automates a download. I have tried some of the answers on the forum, where they have advised using quotes and the operation: `wget "url/?target=link"`. This however does not work in the following instance.

Basically, I want to download the file labeled 'P1TLH.bam', in the download tab of the link provided:

# Answer
> 1 votes
If you want to download a file, you need to use the link *to that file*. Your original attempt, `wget https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR7276474` wouldn't work since that's a link to the trace page of the relevant run. If you want to download something else, just right click on the link (the one in the screenshot in your question), copy the URL and use that:

```
wget https://sra-download.ncbi.nlm.nih.gov/traces/sra64/SRZ/007276/SRR7276474/P1TLH.bam

```

# Answer
> -1 votes
NCBI has its own digger for this purpose. No, my bad the traditional Perl approach of LWP and HTTP::Request::Common to interrogate NCBI, isn't going to work no how much data munging or spider crawling is done, because its just an NGS data dump. ... well its diggable, but it would be easier to dump the data and interrogate it locally.

---
Tags: phylogenetics
---


# Question
Title: Parallel version of cn.mops' bamDataRanges
I am trying to run cn.MOPS to detect CNVs in multiple bam files. Following the example code, I have run:

```
source("https://bioconductor.org/biocLite.R")
biocLite("cn.mops")
options(width=75)
set.seed(0)
library(cn.mops)
library(Biobase)
library(GenomicRanges)
library(GenomeInfoDb)
cn.mopsVersion <- packageDescription("cn.mops")$Version
library(cn.mops)
BAMFiles <- c('bam1.bam', 'bam2.bam', ..., 'bamN.bam')
bamDataRanges <- getReadCountsFromBAM(BAMFiles)

PLEASE BE PATIENT... this might take a while. Consider using the parallel version of this function

```

Since this is indeed taking a while, I'd like to follow the tool's suggestion and use "the parallel version of this function". But how do I do that? I don't see any mention of "parallel" in the manual. Is there some standard Bioconductor/R trick to make a function parallel? Am I missing something in the manual?

How can I run `bamDataRanges` parallelized?

# Answer
This is an option for the function documented in its help page. Just run `?bamDataRanges` and you will see:

```
Usage:

     getReadCountsFromBAM(BAMFiles, sampleNames, refSeqNames, WL = 25000,
       parallel = 0, ...)

[...]

parallel: The number of parallel processes to be used for this
          function. Default=0.

```

So all you need to run with 10 cores is:

```
bamDataRanges <- getReadCountsFromBAM(BAMFiles, parallel=10)

```
> 0 votes

---
Tags: cnv
---


# Question
Title: General framework for fusing biological networks
I know that there are different types of biological networks. Protein-protein interaction, gene interaction, metabolic networks, gene expression networks, etc. Each of these networks can be represented as a graph.

1) Is there a principled way to integrate these graphs into a single large graph? For example, probably most of these networks have nodes that are genes, or gene related things, and this could allow us to match nodes in one network to nodes in another network.

2) Is there a principled way to ask joint questions about these networks?

# Answer
I believe you're referring to something related to multi-layer networks. Here you can find an interesting reading about this kind of networks:

https://academic.oup.com/comnet/article/2/3/203/2841130

Also you can try with this one:

https://www.sciencedirect.com/science/article/pii/S0370157314002105

As a matter of fact, I'm not an expert of multi-layer networks but, as far as I can tell, the main problem (as you mentioned) is to find *correlations*, or *connections*, between the different layers in order to be able to draw some meaningful conclusion.

I am not sure you can call this a multi-layer network but this can be an example of what you're looking for: suppose you have a protein-protein interaction network. You can enrich your network using biological information. You can associate biological functions, or cellular locations or whatever you can find in GeneOntology (just to make an example of a biological database), as a further layer of information regarding your protein-protein interaction network.

This will help you in elucidating the role a node (or more nodes) is playing with respect to its neighbours: let's say five nodes form a clique and you also know they share the same cellular location...then you're able to say something about these nodes (maybe they form a complex?) that has a dual nature: biological and *graph-based*.

Hope it helps!
> 2 votes

# Answer
You are talking network statistics and this is a complex area of stats.

For generalised simplicity, Google fusion tables have a simple to use but impressively interactive network tool. This is a great starting point for exploration. Have a go and see what you think.

If you are looking at a haplotype map, there is a formal solution to this (parsimony network), which is standalone program written some time ago but still useful today. In my area split-decomposition is an important important technique.

The easiest question is the number of connecting nodes and Google fusion tables is as good as any routine interactive networking tool. As I say beyond this network stats is a specialised field.
> 1 votes

---
Tags: machine-learning
---


# Question
Title: Percentage distribution of cells in all clusters based on their treatment condition?
I have 2151 cells, I clustered them by Seurat to 5 clusters. With the code below, I am able to have the number of cells per cluster and per condition:

```
number_perCluster<- table(object@meta.data$conditions,
                          object@meta.data$clusterID)

```

Could I please get a hint on how I could proceed to produce a figure where I say: X% of cells treated with condition "Y" are located in cluster "C" and so forth?

Thank you in advance.

# Answer
Here is a solution using `dplyr` and `ggplot2`:

```
library(Seurat)
library(dplyr)
library(ggplot2)

meta.data <- pbmc_small[[]]

# create random classifications for the sake of this example
meta.data$condition <- sample(c('A', 'B', 'C'), nrow(meta.data), replace = TRUE)

counts <- group_by(meta.data, condition, res.1) %>% summarise(count = n())

ggplot(counts, aes(res.1, count, fill = condition)) +
  geom_bar(stat = 'identity')

```
> 4 votes

---
Tags: r, scrnaseq, seurat, ggplot2
---


# Question
Title: How can I interpret gene expression data from Bioconductor packages?
I am currently looking at microarray data from a bioconductor microarray dataset. Specifically, I have data (a snippet) which looks like the following:

```
   GSM794885   GSM794886   GSM794887   GSM794888   GSM794889
   -1.67170290 -0.07803822  0.99990560 -1.94692850 -1.78671360

```

I am wondering how to interpret the values in this microarray dataset. Is it generally the case that larger values mean a more significant effect?

# Answer
Processed microarray data is commonly represented as log\_2(expression). This data transformation is used because the data more closely fit a normal distribution in log space. With such a transformation, negative values represent very small amounts, rather than the absence of something.

Larger values mean that more stuff is being produced (presumably transcripts), but its significance can only be assessed when combined with other data (e.g. through differential expression analysis).
> 2 votes

---
Tags: rna-seq, gene, rna, gene-expression
---


title: High Phred Quality score VCF file

text:
# Question
Title: High Phred Quality score VCF file
I am working on ION Torrent S5 XL and I am facing to an issue with VCF output generated by the variant caller of thermofisher. It stating on VCF 4.1 that Phred Quality Score is on a range of 0 to 255. In my VCF I have values over than 255. Do you know if there is a way to transform these values?

# Answer
> 2 votes
See this post on the GATK forums.

It is normal to see values above 255. The QUAL score is the posterior probability that all samples in your callset are homozygous reference. I quote Sheila's response from the thread posted above:

> The QUAL score is the Phred-scaled posterior of AC = 0. We use the AC priors and the PLs to get the likelihood of the data given each AC, then use those to get the posterior probability for each AC. From there, the calculation is 1 - Pr{AC \> 0}.

Please see the GATK post in Methods and Algorithms for the nitty gritty details. Beware, Bayesian statistics abound.

# Answer
> 0 votes
GATK's `HaplotypeCaller` and `freebayes` are two more variant caller that output QUAL values greater than 255. The reason for that is, that each read have its part on calculation that value. This results in the higher the read depth at this position the higher the QUAL value.

To get something useful out of this value, the usual thing that is done, is to calculate the quality per depth (QD) by just divide the QUAL value by the read depth.

Note whenever you try a different variant caller, you have to find out whether the QUAL value is read depth dependend or not.

---
Tags: vcf
---


# Question
Title: How to correctly call a VCF file using damaged DNA? (IonTorrent & FFPE)
EDIT: I am updating this question to make it more specific to my issue.

### For context - original question prior to edit:

How do I obtain a deamination metric when doing the variant calling using the IonTorrent variant caller, and secondly, how do I correct my called variants for deamination to ensure that these don't provide false positives for the downstream analysis?

## Background

I'm calling variants using the IonTorrent TorrentSuite on DNA which has been sequenced from formalin-fixed paraffin-embedded (FFPE) tissue. This has a major issue in that without addition of uracil-N-glycosylase, some of the Ts in the original DNA are deaminated to Cs, which upon sequencing and calling variants can show up as mutations, either as T\>C transitions or G\>A (from the opposite strand, due to PCR in the library prep). I do not have any idea how long these samples were stored without UNG before sequencing.

TVC (Torrent Variant Caller) gives a deamination metric (essentially, sum of T\>C and G\>A variants over all variants called), and for our samples, the highest value seen is ~0.92. Naively postprocessing the variants show that for these samples, C\>T/T\>C transitions overwhelm the remaining variants among my samples.

## Application

My use case is this: these are medical samples, which have been inspected by a pathologist (hence the FFPE treatment), and I want to determine which variants are *predictive* of outcome, hence I have two potentially contradictory goals: reduce false positives and capture the rarer variants which may hold predictive power.

## Question

* Given the IonTorrent variant calling pipeline (sequencing \> BAM file \> TVC \> VCF file with deamination statistic):
    + **(Q1)** Is there a way of correcting the output VCF to remove transition and deamination errors?
    + **(Q2)** Alternatively, is there a set of filters to use in bcftools to reduce the effect of the errors on how variants are called?

# Answer
Sites that were deaminated are often called abasic sites.

## Q1: How to obtain DNA damage metrics?

You might want to check out mapDamage \- this software is used to perform these sorts of analyses and produce plots of types of DNA damage along sequencing reads.

## Q2: How to correct for abasic sites and other DNA damage errors?

You can use your choice of kmer-based or alignment based error correction software. The errors will be random and if the sequencing depth is sufficient the erroneous kmers at the ends of reads will be corrected. However, depending on what you are trying to do downstream with your reads you may be better served by not correcting the reads, then dealing with transitions/damage later.
> 2 votes

# Answer
Having inspected one of my samples (with estimated deamination metric, the proportion of called SNPs consisting of C\>T or G\>A transitions, of 0.837), this looks like a straightforward filtering procedure. Looking at the minor allele fraction (AF) and the quality score (-10 log10(p-value)) per type of SNP: it can be seen that the C\>T and G\>A transitions likely caused by deamination (and PCRing up the deaminated DNA to get the libraries) all cluster in the bottom right of the plot - low AF and low QUAL, so setting thresholds to the left (higher QUAL) and above (higher AF) the groups will remove these from the called variants.
> 1 votes

---
Tags: variant-calling
---


# Question
Title: Calculate the percentage of each unique phylogenetic tree in a BEAST output
I have a nexus formatted BEAST output containing 20,000 phylogenetic trees of seven taxa. Is there any way to get the percentage of each unique phylogenetic tree contained in this output?

I already made an unsuccessful attempt with .

# Answer
I finally managed to do it in . Here is my code:

```
install.packages('devtools')
library(devtools)
install_github('santiagosnchez/rBt')
library(rBt)

beast_output <- read.annot.beast('beast_output.trees')
beast_output_rooted <- root.multiPhylo(beast_output, c('taxon_A', 'taxon_B'))
unique_topologies <- unique.multiPhylo(beast_output_rooted)

count <- function(item, list) {
  total = 0
  for (i in 1:length(list)) {
    if (all.equal.phylo(item, list[[i]], use.edge.length = FALSE)) {
      total = total + 1
    }
  }
  return(total)
}

result <- data.frame(unique_topology = rep(0, length(unique_topologies)),
                     count = rep(0, length(unique_topologies)))
for (i in 1:length(unique_topologies)) {
  result[i, ] <- c(i, count(unique_topologies[[i]], beast_output_rooted))
}

result$percentage <- ((result$count/length(beast_output_rooted))*100)

```
> 4 votes

# Answer
You require a consensus phylogeny.

This is available in RAxML, PAUP, the old Phylip package and possibly the stats companion to Beast, Tracer. I'd need to dig up the code for RAxML and PAUP. The easiest is Consense here because you just download it and input the sequence. NOTE their is possibly a format difference (RAxML and Consense use Phylip format), I suspect you've got nexus, or Newick format.

Beast is likely to have its own consensus formating OR the Beast companion "Tracer" may have this (I can't remember). Tracer is here and would be cool if it did consensus phylogeny because there will be no conversion problems (it doesn't use Phylip format). The alternative would be to use PAUP (no conversion problems), which is currently in a state of flux from various stages of freeware and updates. PAUP is contained within the Genious package, which is available as a 1 month free-ware.

Have a look and let me know: because if you want to use RAxML or Consense you'll likely need to remind me of the Beast output, if it is what I think it is a format conversion app is needed. The calculation is easy, its easy to script a conversion program, but I think tree formating will be an issue.

You could use Ape, within R, but as you have gathered there are alot of other solutions. There are more recent, impressive, solutions as stand-alones which I've not used yet.
> 2 votes

---
Tags: r, sequence-alignment, phylogenetics, phylogeny, beast
---


# Question
Title: Python - Finding a motif - input: a txt file with 10 sequences and 10 motifs
When I run my BruteForce function with only one input it works and the result is correct.

```
def BruteForce(s, t):
    occurrences = []
    for i in range(len(s)-len(t)+1): # loop over alignment
        match = True
        for j in range(len(t)): # loop over characters
            if s[i+j] != t[j]:  # compare characters
                match = False   # mismatch
                break
        if match:   # allchars matched
            occurrences.append(i)

    print(occurrences)

t = 'CAACTTCCA'

s = 'GACAACTTCCAACTTCCAACTTCCCGTCCCAACTTCACAACTTCGGCCCAACTTCCATGCAACTTCACCATCAACTTCGCTCGAAGCTGCCTTCCACTCCAACTTCACAACTTCCTCAACTTCCTCACCAACTTCAGCAACTTCTCTAGGGCCAACTTCCAACTTCTCAACTTCTCAACTTCCAACTTCCGACAACTTCTCCTGGCAACTTCCAACTTCCAACTTCAATACAACTTCGCAGACAACTTCCGCAACTTCGAACAACTTCCAACTTCCCCAACTTCCAACTTCCAACTTCGCCAACTTCCAACTTCCAACTTCCCAACTTCAGATAGCAACTTCGATCTTACACAACTTCACGCAACTTCTCCAACTTCCAACTTCTGTGCAACTTCTCTGAACAACTTCCTCAACTTCCAACTTCGCAACTTCCCCAACTTCCTCAACTTCATGCAACTTCGAGGCAACTTCCCAACTTCGCAACTTCCTATTCCCAACTTCTGTGGCAACTTCTCAACTTCTGGACAACTTCTATGCCCAACTTCACAACTTCCCCAACTTCTTTACAACTTCGACAACTTCATCAACTTCTAGTCAACTTCTGGTCCAACTTCCAACTTCCCCAACTTCCAAAGTGCCGCAACTTCGTAACAACTTCACGCGCTCAACTTCAACCAACTTCTTTTCCCGCAACTTCGCAACTTCACAACTTCTAATCAACTTCCAACTTCGGATCAACTTCCAACTTCGCCAACTTCCAACTTCCAACTTCTCCAGGGACAACTTCAAGTACAACTTCCAACTTCGCAACTTCACAACTTCCCAACTTCGCAACTTCTACACGCAACTTCCAACTTCTGGTCCCAACTTCATCAACTTCAGTCAACTTC'

BruteForce(s, t)

[2, 9, 48, 152, 175, 205, 212, 261, 277, 284, 300, 307, 370, 410, 607, 623, 717, 735, 751, 758, 792, 844]

```

But when I use it in the following loop it returns nothing! I have used the print() to see if it reads the strings and motifs from the input\_2.txt file and it does but it returns an empty list of occurrences \[\].

```
with open('input_2.txt', 'r') as input_2:
    nt = input_2.readline() #nt is the no of strings+motifs in 1st line
    nt = int(nt) #in this case it is 10

    for i in range(nt):
        s = input_2.readline()
        print(s)
        t = input_2.readline()
        print(t)
        BruteForce(s, t)

```

Could you please let me know what changes are needed here? Many thanks, Beeta

input\_2.txt:

```
10

GACAACTTCCAACTTCCAACTTCCCGTCCCAACTTCACAACTTCGGCCCAACTTCCATGCAACTTCACCATCAACTTCGCTCGAAGCTGCCTTCCACTCCAACTTCACAACTTCCTCAACTTCCTCACCAACTTCAGCAACTTCTCTAGGGCCAACTTCCAACTTCTCAACTTCTCAACTTCCAACTTCCGACAACTTCTCCTGGCAACTTCCAACTTCCAACTTCAATACAACTTCGCAGACAACTTCCGCAACTTCGAACAACTTCCAACTTCCCCAACTTCCAACTTCCAACTTCGCCAACTTCCAACTTCCAACTTCCCAACTTCAGATAGCAACTTCGATCTTACACAACTTCACGCAACTTCTCCAACTTCCAACTTCTGTGCAACTTCTCTGAACAACTTCCTCAACTTCCAACTTCGCAACTTCCCCAACTTCCTCAACTTCATGCAACTTCGAGGCAACTTCCCAACTTCGCAACTTCCTATTCCCAACTTCTGTGGCAACTTCTCAACTTCTGGACAACTTCTATGCCCAACTTCACAACTTCCCCAACTTCTTTACAACTTCGACAACTTCATCAACTTCTAGTCAACTTCTGGTCCAACTTCCAACTTCCCCAACTTCCAAAGTGCCGCAACTTCGTAACAACTTCACGCGCTCAACTTCAACCAACTTCTTTTCCCGCAACTTCGCAACTTCACAACTTCTAATCAACTTCCAACTTCGGATCAACTTCCAACTTCGCCAACTTCCAACTTCCAACTTCTCCAGGGACAACTTCAAGTACAACTTCCAACTTCGCAACTTCACAACTTCCCAACTTCGCAACTTCTACACGCAACTTCCAACTTCTGGTCCCAACTTCATCAACTTCAGTCAACTTC

CAACTTCCA

ATTGTGTATCGCTATGTATCGTGTATCGGATTTTGTATCGTGTATCGGTGTATCGTGTATCGTGTATCGTATCACTGTATCGTTGTATCGTAGCGTTGTATCGTAATGTATCGCTCTGTATCGTGTATCGGGTTTGTATCGATGTGTATCGCTGTATCGGTGTATCGGTGTATCGCTTGTATCGACAGCGCTTGTATCGTGTATCGACCTGTATCGGTGTATCGTGTATCGAATGTATCGTTGTATCGAATTGTGTATCGTGTATCGTGTATCGTGTATCGTATGTATCGACTGTATCGCTGTATCGTAGCCTGTATCGTGTATCGGTGTATCGGTGTATCGGCGTGTATCGAATGTATCGTTGTATCGCTGTATCGTGTGCTGTGTATCGGTGTATCGACCCGTTGTATCGTGTATCGTGGGTAAGTGTATCGTTGTATCGTAGTGTATCGTGTATCGTGTATCGTCATGGTATTGTATCGTTATAGCTGTATCGCTCGGCTGTATCGTATGTATCGTGTATCGTGTATCGCATGTATCGGTGTATCGCTGTATCGACCATGTATCGTGTATCGGTGTATCGATGTATCGCTTCCATTAGAAATGTATCGTATGTATCGTGTATCGCTGTATCGTTTGTATCGCATGTATCGATTGTGTATCGTGTATCGTTGTATCGTGTATCGTGTATCGTTGTATCGTAATAACGATGTATCGAATGATGTATCGTTGTGATGTATCGTGTATCGACTGTATCGATGTATCGGGTTGTATCGTGTATCGTGTATCGTAGATGTATCGAGAGCCATTGTATCGCCTTGTATCGCGTGTATCGGTGTGTATCGTTTGTATCGTGTTTGAATGTATCGCTGTATCG

TGTATCGTG

TTTCAGCGTGTTTGGCTTTCAGCACTACTTTCAGCGGTTTCAGCTTTTTCAGCTTTTCAGCTATTTTCAGCTTTCAGCCATCTTTCAGCCCAGGTTTCAGCTGTTTCAGCTTTCAGCGAGTTTCAGCTTTCAGCGCGGGGATTTCAGCGTATTTCAGCGACATTTCAGCAGTTTCAGCTGGAGTGAAGCCGTTTCAGCGCTTTTCAGCACCATTTCAGCTAGTAGGTTTCAGCTTTTCAGCCTTTTCAGCTTTCAGCTATTTCAGCTCTAGAAGTCGTTTCAGCTTTCAGCTTCCTTTTCAGCTTTCAGCTTTCAGCCTTTCAGCTTTCAGCAAACTTTCAGCGATTTCAGCTTTCAGCTTTCAGCTTTCAGCTTTCAGCATTTCAGCTTTCAGCCGTTTTCAGCGACTCTTTCAGCGTTTCAGCTTTCAGCTGATCGTTTTCAGCTTTTCAGCCGGTTTCAGCGAAAGTTGGTCTTTTCAGCCAAATTTTCAGCTTTCAGCTTTCAGCGTTTCAGCTTTTCAGCCATTTCAGCCTATTTTCAGCAATCTTTCAGCAATTTTCAGCGCAATTTCAGCCAAAATTTCAGCATTTCAGCTTTCAGCGTAGTTTCAGCTTTCAGCTGAGCCTTTCAGCTTTTCAGCTCTCTTTTCAGCTTTCAGCGATTTCAGCCTTTTCAGCGCTTTCAGCACCACGCCTTTCAGCTTTCAGCTTTTCAGCCTTTCAGCCTATTTTCAGCGGTTTCAGCGTTTCAGCTTTCAGCACTTTCAGCTTTCAGCTTTTCAGCGGAAGGTTTTCAGCA

TTTCAGCTT

AGCATCCGACTGCATCCGTGCATCCGACGGCATCCGGCATCCGCAGCATCCGCACAGCCGCATCCGAGCCCGCATCCGAAGCATCCGGCATCCGGGCATCCGGCCGCATCCGGCATCCGGGCATCCGTGCATCCGGCATCCGTGTGCATCCGGCATCCGGGGTGCTTGCATCCGCGCATCCGTGCATCCGCGCATCCGGCATCCGCAGTCCGCATCCGCGCATCCGGCATCCGAGCATCCGTATTCGCATCCGTCGAGCATCCGTACGTCGCATCCGTGCATCCGGGAGAGGCATCCGGGCATCCGCGCATCCGGTGGACTATAACGCGCATCCGGCGAGCATCCGCCGCGGCATCCGTGCATCCGGTACGGCATCCGGGCTGAGGCATCCGCTCCCGGGGCATCCGCGGCGCATCCGCACGCGCATCCGGCATCCGAGCATCCGTTGCATCCGGAGCATCCGTAGCATCCGATAGCATCCGCAAGCTTGCATCCGGCATCCGCTACAGCATCCGGTGGCATCCGAAGCATCCGCGCATCCGCGGATTGCATCCGCACACGGCATCCGTAGCATCCGTGCATCCGGAGTGCATCCGCCGGGCGCATCCGCTGCATCCGCAGCGCATCCGAGCATCCGATCTTGCATCCGGCATCCGGGAATGCATCCGGGACTGCATCCGGTCTTAAGGGTGCATCCGAATGCATCCGCTGTAGCATCCGGCATCCGAGTAAGCATCCGAGTTCTGCATCCGGAAGCAGCATCCGGCATCCGGACACCAGCATCCGCGCATCCGGGGCGAGCATCCGAGCATCCGGCATCCGGGGCAAGTGGCATCCGGCATCCGTCGTGCATCCGGGCATCCGAA

GCATCCGGC

GACGAGCTGACGAGCGTGACGAGCAGGCGACGAGCGACGAGCATGCGGTGACGAGCGGACGAGCGACGACGAGCGCGGACGAGCGACGAGCGGACGAGCTCCGACGAGCGACGAGCACGACGAGCGCCCGACGAGCGACGAGCATGATGACGAGCGACGAGCTTAGACGAGCATGCTGACGAGCTGACGAGCTGACGAGCTCGTACATCGACGAGCGGACGAGCAGCCCGATAAGCCTTCGACGAGCTGACGAGCCACGACGAGCGACGAGCGACGAGCGGACGAGCGGGACGAGCGTCGGACGAGCGGACGAGCAGGACGAGCACCTCAATCGACGAGCTGACGAGCGGACGAGCCGACGAGCTGAACTGAAGGACGAGCCTACGTTCTAACGTGCCGTCACTGACGAGCGACGAGCGGACGAGCGATGGACGAGCTGACGAGCGCAGGACGAGCGAGAAGGCTGACGAGCAGACGAGCTGACGAGCTCGACGAGCAGACGAGCGGACGAGCTTTGATAAGACGAGCACGTCGACGAGCCTCGTGACGAGCTGACGAGCGACGAGCTAGACGAGCGACGAGCTGACGAGCGACGAGCACTAAACGCGACGAGCTCGACGACGAGCATGGACGAGCGACGAGCTAGGACGAGCGGCGACGAGCAGACGAGCGACGAGCTGACGAGCTTATAGACGAGCCTGACGACGAGCCAAGACGAGCGACGAGCGGACGAGCGACGAGCACGACGAGCCGACGAGCCGCGGACGAGCCGTAGACGAGCCAATCATTAAGACGAGCCGACGAGCCATTTGGGGACGAGCCTCCTCGACGAGCCATAAGACGAGCTGACGAGCCATGACGAGCATGCCCGACGAGC

GACGAGCGA

GAGGACCCCGGACCCCTGGACCCCAGGACCCCGACGGTGGGGCGGACCCCGAGAAGGACCCCTGGACCCCCGCAGGACCCCTTTATCGGACCCCGGACCCCGGACCCCCGGACCCCGGCTGGACCCCGTCGTAAGGACCCCGAATCGGACCCCTAGGACCCCAGGGACCCCTCCCCGGTTGGACCCCCTGGACCCCCTTGAGAGGGACCCCGGACCCCACGCCGTGCTTAAGGACCCCACTATGGACCCCATACGGACCCCGGACCCCGATCAGAGCGACCAGGACCCCCGGCCTGGACCCCTCGGACCCCAGAAGGACCCCTCGGACCCCTGGACCCCGGACCCCACAGGGACCCCGGACCCCTAGCCGCGGTGGACCCCCGGACCCCCGCAGATGGGGACCCCTCATGCGGACCCCCTGACGGACCCCTGGACCCCCGGGACCCCCAGGGACCCCATTTCGAGGACCCCATGGGGACCCCAAGCTGGACCCCGGACCCCTGGGACCCCCGGACCCCGGACCCCGGACCCCATAGGACCCCCGTTTGTTGCCATGGACCCCTTGGGACCCCGAGTCGGACCCCGGACCCCCAGGACCCCACAGGACCCCGGGACCCCGGACCCCATCGATCGGACCCCGGACCCCGGACCCCTACGGACCCCGCTAAGGGACCCCGGACCCCGTGGACCCCCACGGACCCCTAGGACCCCGGGACCCCGGACCCCGGACCCCAGTTGGACCCCCCTTAGGACCCCAGGACCCCACATGAGGACCCCGGGACCCCGGGACCCCTGGACCCCTGGACCCC

GGACCCCGG

AAACCTTCGGACCTTCGTACCTTCGACCTTCGTGAAACCTTCGGACCTTCGTGACCTTCGGGTGACCTTCGGGTACTCACCTTCGGCACCTTCGTATACCTTCGAGACCTTCGCTGGACTGTAAACCTTCGATCGTTATTTTCTACCTTCGCACACCTTCGACCTTCGAGCCGAACCTTCGGACTGGCCTCACCTTCGATCACCTTCGACCTTCGAACCTTCGACCTTCGACCATACCTTCGACCTTCGGACCTTCGGACCTTCGCACCTTCGAACCTTCGCCACCTTCGACCTTCGAACCTTCGGTGGGCCACCTTCGTTACCTTCGAACCTTCGGGCACCTTCGACCTTCGGCTTACACCTTCGGACCTTCGATGACTGACCTTCGAAACCTTCGACCTTCGACCTTCGATTCCACACCTTCGACCTTCGACACCTTCGACCTTCGCGTTCACCTTCGACCTTCGACCTTCGGGCCAATACCTTCGACCTTCGGCCTAATACCAACCTTCGAACCTTCGAACCTTCGGGTGACCTTCGGACCTTCGCGCACACCTTCGCACTTACCTTCGACCTTCGGTGACCTTCGACCTTCGTGCACCTTCGCGTCTCACCTTCGCACCTTCGTTACCTTCGAACCTTCGACCTTCGACAACCTTCGTACCTTCGCTACCTTCGGGAACCTTCGATGTCCTACCTTCGACCTTCGCGACCTTCGACCTTCGTAACCTTCGAACCTTCGACCACCTTCGAACCTTCGAAGGAAGGACCTTCGAATCATTACCTTCGTTACCTTCGTAGTACCTTCGCCGACCTTCGGACCTTAACCTTCGTGA

ACCTTCGAC

AAGTGGTATTCCCAAAAAGATATTCCCGTTATTCCCGCTATAAATATTCCCAACTGTTTATTCCCACTTATTCCCGTATTCCCGTATTCCCCGGGCCTTTTATTCCCCATTCTATTCCCGTATTCCCATATTCCCACTTATTCCCCTATTCCCATATTCCCTATTCCCTATTCCCGTAAGTATTCCCCCCCTTATTCCCCCATTTGTCCATATTCCCTCAATATTCCCTATTCCCCTATTCCCGCTTCCCCAGTCTATCCCAAATTGTATTCCCAATATTCCCCTTTATTCCCTTTATTCCCGTATTCCCTATTCCCGTGGCTATTCCCCCGACGTTATTCCCTCTATATTCCCGTCTCAAAGTATTCCCTTATTCCCAGTATTCCCATTCGCCTGGTATTCCCGTATTCCCGGTATTCCCGTTTATTCCCATCTATTCCCATTATTCCCTGGTATTCCCGCTATTCCCGCAACGGTATTCCCTATTCCCTATTCCCGTATTCCCAATATTCCCTTCCTATTCCCGTATTCCCTAGCTATTCCCATATTCCCGATATTCCCTATTCCCGTATTCCCTATTCCCTTATTCCCTATTCCCTAGTTATTCCCAAACAGGCTATTCCCTATTCCCTATTCCCATATTCCCTTATTCCCTACGTATTCCCTCTTAAAAATATTATTCCCCTATTCCCTATTCCCGTATTCCCCGAAAGTATTCCCCCATATTCCCAATTATATTCCCGGTATTCCCCCGTATTCCCTATTCCCTTATTCCCATATTCCCATATTCCCAAAGTTTATTCCCGTATTCCCCACTATTCCCATAGATTATTCCCTATTCCCTATTCCCTTATTCCCCTATTCCCTATTCCCTTATTCCCGTATTCCCAGTTTATTCCC

TATTCCCTA

GAAAGGATAAAGGATGTCAAAGGATGCAGATATAAAGGATAAAGGATACAAAGGATTAAGTATGTCCGAAAAAGGATAAAGGATTAAAAGGATCGCTGAACCTTACACAAAGGATAGTGAACAAAGGATTCAATAAAAGGATAAAGGATCAAAGGATCAAAAGGATACAAAGGATGCCGATGAAAGGATAAAGGATCTAAAGGATAAAGGATGAGAAAGGATTGGAAAGGATTATGAGATCAAAGGATAAAGGATGGTGTAAAAGCTAAAGGATTCGGCAAAGGATAAAGGATTAAAGGATAAAAAGGATAAAGGATGAAAGGATCCGCAGGGACCAGCAAAAGGATAAAGGATCGAATGGGTAAGAAAGGATCAAAGGATGAAAAAGGATTACTAAAGGATAAAGGATCCTGAAAGGATTACAAAGGATCTTAAAAGGATTCGGAAAGGATCCATAGGAAAAGGATAAAGGATCGCGAAAGGATAAAGGATTAAAGGATAAAGGATAAAAGGATTCAAAGGATAAAGGATAGACAAGGGAGAAAGGATGCAAAGGATGAAAGGATAAAGGATAAAGGATAAAGGATTAGGTTAAAGGATCGAAAGGATCCAAAGAGAGCGAAAAGGATGAGGACGAAAGGATCAAAGGATCCAAAGGATCAAAGGATAAAGGATTGCGATGGAAAGGATGTCAAAGGATCACCAAAGGATAAAGGATAATAAAGGATAAAGGATAAAGGATCAAAGGATTTGAAAAGGATAAAGGATCACGAAAGGATAAAAGGATTGGCAAAGGATGAAAAGGATGAAAGGATAAGCCCTCCAAAGGAT

AAAGGATAA

ATTCAAATACTTCAAATAGTCAAATATCAAATATGCTTCAAATATCAAATATCGGTCAAATAAGTCAAATAAGTCAAATATCAAATATCAAATACTCAAATATCAAATAGTGTCAAATACGAATTGGGTCAAATATCAAATATGTGTTCAAATATTCTTCAAATACTGGACACTCAAATAGGAGTCAAATATCAAATATCAAATACGTGAAGTGCTTGTCAAATATTTCAAATACTTTCAAATATTCAAATACTCAAATATGTTCAAATATCAAATATCAAATATTTCAAATATCAAATATCAAATAGTCCTCAAATAGCAAACCAGTCAAATATCAAATATGTCAAATATGCTCACGGCAACCTCAAATATCAAATATCAAATAGATGTCAAATAAGTTTCAAATATCAAATATCAAATATCAAATATCAAATATTCAAATATCAAATAGCGGGCTCAAATAAGCTCAAATACGTCAAATAGGGGGTCAAATAGATCAAATAGTCAAATATTCAAATACATCAAATAGTCAAATACAAGAACCACCGAGATCTCAAATAATCAAATATCAAATATGGTCAAATAAATATTCAAATAGGTCAAATAAGATCAAATATCAAATAAGTCGTCCATCAAATAGTCAAATAGCTCAAATAACCTCAAATATCAAATATTCAAATACCGGTCATCAAATAGGACAAATCAAATAGTCAAATAAGATCCTCTCAAATAATTCAAATAGCTGTTCAAATACTCAAATATCAAATAGTCAAATATCAAATATTCAAATATCAAATACTTCAAATATGTTCAAATAATCAAATACTCAAATATCAAATAATCAAATAATACTTCAAATATACCAAACGCTCAAATATTAGTTGGATCAAATATCTTCAAATATCAAATAA

TCAAATATC

```

# Answer
Let's start by doing some explaining about why your code doesn't work.

```
with open('input_2.txt', 'r') as input_2:
    nt = input_2.readline() #This is reading the first line of the file (which is 10)
    nt = int(nt) # this is storing "10" as an int

    for i in range(nt): # This is the same as "for i in range(10):" which just returns nums 0-9, does not access lines of your file
        s = input_2.readline() # this is reading the second line which is blank
        print(s) # this will print a blank line 10 blank lines
        t = input_2.readline() # this is reading the rest of the file
        print(t)
        BruteForce(s, t) # this is calling your function on 10 blank lines and the restof your file, not each seq/motif

```

I made some modifications that will help you achieve your goal I think:

```
def BruteForce(s, t):
    occurrences = []
    for i in range(len(s)-len(t)+1): # loop over alignment
        match = True
        for j in range(len(t)): # loop over characters
            if s[i+j] != t[j]:  # compare characters
                match = False   # mismatch
                break
        if match:   # allchars matched
            occurrences.append(i)

    print(occurrences)

def group(lst, n):
    """ This function will be used to group your sequences and motifs later """
    for i in range(0, len(lst), n):
        val = lst[i:i+n]
        if len(val) == n:
            yield tuple(val)

with open('input_2.txt') as input_2: #standard open in read mode,
    lines = input_2.readlines()[2:] # skip 10 and the blank space bc we don't care about them
    split_line = [line.strip('\n') for line in lines if line!='\n'] #this will store lines that aren't blank into a list and strip new_line chars '\n'
    seq_tuple = list(group(split_line, 2)) # this will use group func above to group the list into tuples of 2's which are sequence and motif.

for x in seq_tuple: # this accesses the list of tuples (seq_tuple)
    BruteForce(x[0], x[1]) # x[0] in each tuple is seq, or BruteForce(s), x[1] in each tuple is motif, or BruteForce(t)

```

Running this will output Each motifs start position in each sequence:

```
[2, 9, 48, 152, 175, 205, 212, 261, 277, 284, 300, 307, 370, 410, 607, 623, 717, 735,
751, 758, 792, 844]
[14, 33, 48, 55, 116, 192, 217, 253, 260, 267, 312, 370, 405, 412, 445, 452, 511,
518, 561, 613, 657, 672, 679, 735, 769, 776, 847]
[37, 46, 64, 103, 120, 226, 243, 277, 284, 296, 303, 318, 345, 352, 359, 366, 381, 418, 439, 488, 495, 510, 590, 608, 628, 648, 697, 704, 748, 764, 771]
[29, 80, 95, 105, 128, 145, 191, 219, 328, 423, 489, 642, 713, 758, 808, 831]
[28, 57, 77, 102, 129, 148, 258, 265, 404, 419, 449, 553, 569, 584, 629, 664, 713, 728]
[2, 87, 94, 109, 204, 254, 332, 350, 487, 511, 518, 577, 603, 611]
[18, 154, 201, 216, 223, 235, 283, 339, 390, 397, 418, 425, 434, 453, 460, 481, 565, 582, 637, 644, 695, 711, 735]
[154, 161, 221, 303, 476, 483, 526, 554, 569, 584, 591, 617, 624, 647, 685, 754, 829, 836, 859]
[1, 33, 70, 135, 182, 198, 241, 279, 294, 303, 340, 395, 460, 478, 493, 500, 517, 559, 566, 573, 661, 704, 711, 721, 728, 754, 772, 808]
[19, 37, 44, 73, 80, 95, 128, 184, 191, 263, 270, 286, 293, 327, 364, 371, 400, 407, 414, 421, 436, 561, 609, 659, 754, 769, 784, 826, 886, 896]

```
> 5 votes

# Answer
This looks hard work. In Perl the function "index" and "rindex" as well as "substr" will perform this without the incremental loop. There must be an equivalent in Python?

I've not tried this, but at a guess it would be something like,

```
my @array;
my $target = 'AGATAGAGAGATAGAGAGATAGAGAGATAGAGAGATAGAG';
my $query = 'AGAGA';
while(length($target) > 4){
    my $location = index ($query, $target);
    push @array, $location;
    $target = substr ($location, length($target), $target);
}
print join("\s", @array), "\n";

```
> 3 votes

---
Tags: python, phylogenetics, sequence-alignment, phylogeny, motifs
---