title: Converting .bw files to .fq (fastq)

text:
# Question
Title: Converting .bw files to .fq (fastq)
I am able to convert .bw files to .fq(fastq) manually in Shell, but I would like to automate the process coz I have hundreds of .bw file that I need to convert. Till now, I could think of following script, but unable to solve the problem. The .bw files are in H3K27ac folder.

```
import subprocess

path_bw = "~/Documents/H3K27ac"
import os

files = os.path.expanduser(path_bw)
for file in files:
    if file.endswith('.gz'):
        files.remove(file)
files = [x[:-3] if x.endswith('.bw') else x for x in files]

# for .bw files

for file in files:
    bigwigtowig = "bigWigToWig " + path_bw + file + ".bw " + path_bw + file +".wig"
    subprocess.call(bigwigtowig, shell=True)

    wigtobed = "wig2bed < " + path_bw + file + ".wig > " + path_bw + file + ".bed"
    subprocess.call(wigtobed, shell=True)

    getfasta = "bedtools getfasta -fi " + path_bw + "hg19.fa -bed " + path_bw + file + ".bed -fo " + path_bw + file + ".fa"
    subprocess.call(getfasta, shell=True)

    perl = "perl " + path_bw + "fasta_to_fastq.pl" + path_bw + file + ".fa > " + path_bw + file + ".fq"
    #os.system(perl)
    subprocess.run(perl, shell=True)

```

Manually, I used the following to convert the .bw files:

1. `# bigWigToWig file1.bw fileout.wig`
2. `# wig2bed <fileout.wig> fileout.bed`
3. `# bedtools getfasta -fi hg19.fa -bed fileout.bed -fo fileout.fa`
4. `# perl fasta_to_fastq.pl fileout.fa > file1.fq`

All the above packages were preinstalled via shell and works perfectly. hg19 is the genom that i am using. fasta\_to\_fastq.pl is openly available in github and helps to convert .fa to .fq.

MAin thing i require is to convert all the files from .bw to .fq and this is the only way in my research that i could find. I need all the names to be the same as output. For example file1.bw shoud be file1.fq , file2.bw --\> file2.fq and so on.

Do you have any ideas if this is possible or are there easier ways to convert it? It would be awesome to know your view point.

# Answer
> 4 votes
What you are doing makes **no** sense at all. A bigwig file contains intervals and then indicates the coverage of that interval. That coverage again comes from the alignments covering that region. There is no direct connection between the interval and the fastq file other than that the fastq file is the basis for the coverage calculation. Converting the interval to a fasta and then forcing it into a fastq file is technically possible, but has absolutely no practical meaning. You cannot recover sequencing reads from a bigwig, that information is lost once you convert your alignments (usually a BAM file) into a bigwig. You should open a new question describing the problem that this is actually based on.

# Answer
> 1 votes
The code is verbose and I ain't at all sure some of it will work, at least in two places.

The first point is `os.path.expanduser` ... admittedly I prefer `pathlib`, however that is definitely not how I understand that command. My understanding of `expanduser` is if `~` is being used. You definitely want to look at `globb`-ing for the opening part of the code.

```
from pathlib import Path
import os

mypath = Path(os.path.expanduser('~/Documents/H3K27ac'))
files = mypath.glob('*.bw') # or use ('*' + '.bw')

```

I'm pretty certain that is what you want and its a lot easier than lines 1 to 10 in your code. The `files` list will then go straight into the loop.

Also, please check

```
files = [x[:-3] if x.endswith('.bw') else x for x in files]

```

Comprehensions are never logical in Python but convention is always to place the `if` statement after the list of interest.

Thus the following is much closer:

```
files = [x for x in files if not x.endswith('.bw') else x[:-3]]

```

**Note** not tested, so this may throw a bug, but its better. However, using `glob` bing you don't need this statement.

The other issue is the sub.process call ... Normally they are written:

```
os.chdir(path_bw)
file1 = 'file1.bw'
file2 = 'fileout.wig'
subprocess.call(["bigWigToWig",  file1, file2], shell=True)

```

The way you've written is:

* Really hard to read
* I don't think it would work because of the absence of the `[]`, certainly if you are using `popen` it will not work

Also check the following variant:

```
subprocess.call(["bigWigToWig",  file1, file2], text=True)

```

I know its seems a bit weird, but it very likely work (don't ask why)

I agree `subprocess.call` is a good idea in a Python pipeline. Issuing a hard path within `subprocess call` I'm not a fan hence using `os` to get Python to move into the `cwd` is a better approach.

A further point is the way you've used 'file' throughout your pipeline line, hmmmm... okay, but its difficult to read. Its not generally how a pipeline is written, but just saying.

`nextflow` is an increasingly popular pipeline language and you may consider learning that rather than going down a Python route (Python is very good for data science), of course this particularly question could be very easily solved in `bash` (but for some reason it ain't trendy).

---
Tags: python, fasta, gene, fastq, hg19
---


# Question
Title: How to know if the DNA sequence has been assembled and why is it important to know how it was assembled?
I have downloaded my FASTA format files, that have the DNA sequences of the coding region of the genes and the DNA sequence of the complete genome, from NCBI. How can I recognize if these sequences have been assembled or not ?

In particular, I read this paper Genome Sequencing of a Genetically Tractable Pyrococcus furiosus Strain Reveals a Highly Dynamic Genome and I downloaded the sequences that they indicate in **Nucleotide sequence accession number** where they say:

> The fully assembled DNA sequence of COM1 has been deposited in GenBank (accession number CP003685).

I clicked on the link CP003685 and I downloaded the DNA sequences that I mentioned at the beginning by clicking on "Send to", then "Complete Record",then "File" (in Choose Destination), then FASTA (in Format). For the cds file I made the same steps but clicking on "Coding Sequences" and I could do the same to have the sequences of the genes clicking on "Gene Features". But are these sequences assembled? By reading the extract of the article I would say yes... but is there a way to know this from the FASTA files that I downloaded or the NCBI website ?

Moreover, why is it important to know how a sequence has been assembled (so which software has been used) ?

I apologize for my very basic questions and I thank you in advance.

# Answer
> 2 votes
Genomes are always assemblies. Assembly basically means making a consensus sequence from many partially overlapping readouts. As there is currently no sequencing technology that allows to sequence an entire genome or even an entire chromosome with a single super high fidelity read that is without errors we have to combine several reads to form the consensus which is then the assembly.

> Moreover, why is it important to know how a sequence has been assembled (so which software has been used) ?

Tbh, I never cared on how a genome I used was assembled. To be fair I only work with mouse or human, and for people who work with genomes that are not as well curated as these it might indeed matter, but others would need to comment on that, not really my field.

# Answer
> 2 votes
A nucleotide sequence may be a part of a genome assembly. If so, just click the assembly link from the right panel and you'll be redirected to the assembly page.  In this assembly page, it says the assembly only contains one chromosome.  Assembly method is also available for some assemblies.

---
Tags: fasta, genome, assembly, sequence-analysis, phylogenetics
---


# Question
Title: Way to get genomic sequences at given coordinates without downloading fasta files of whole chromosomes/genomes first?
So I have a list of start and stop positions along chromosomes in different species, and I'd like to get the corresponding DNA sequence for each set of coordinates. In the past, I've just download the genome as a fasta file and then use pyfaidx to extract the sequences at the given positions. But now that I'm working with several species at once, I was wondering if there's any kind of tool in Python or R that can fetch your sequences of interest without downloading a bunch of large files. Thanks

# Answer
using a http request.

if there is a DAS server, you can always use this protocol to download the xml -\> fasta. see https://www.biostars.org/p/56/

```
$ curl -s "http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chrM:100,200"  | xmllint --xpath '/DASDNA/SEQUENCE/DNA/text()' - | tr -d '\n'
ggagccggagcaccctatgtcgcagtatctgtctttgattcctgcctcattctattatttatcgcacctacgttcaatattacaggcgaacatacctacta

```

or use the UCSC utility **twoBitToFa** which works with remote files.

If you have an NCBI accession number you can use the E-utilities with seq\_start and seq\_end.

```
$ wget -q  -O -  "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=AE014134.1&seq_start=100&seq_stop=300&rettype=fasta"
>AE014134.1:100-300 Drosophila melanogaster chromosome 2L complete sequence
TGCCAACATATTGTGCTCTTTGATTTTTTGGCAACCCAAAATGGTGGCGGATGAACGAGATGATAATATA
TTCAAGTTGCCGCTAATCAGAAATAAATTCATTGCAACGTTAAATACAGCACAATATATGATCGCGTATG
CGAGAGTAGTGCCAACATATTGTGCTAATGAGTGCCTCTCGTTCTCTGTCTTATATTACCG

```
> 8 votes

# Answer
*I am not aware of any tools that offer this functionality out-of-the-box, although that doesn't necessarily mean it doesn't exist.*

The tool you describe would have to be some kind of web service, where you submit a request with something like the following:

* species
* assembly version
* list of intervals, each with:
    + chromosome ID
    + start position
    + end position

The response would be one or more subsequences in Fasta format that you could save to a file.

This approach pre-supposes that there will be a relatively small number of subsequences, that those subsequences won't be too large, and you won't be requesting them very frequently. Otherwise, this approach will offer very little practical benefit over your current approach.
> 1 votes

# Answer
If you are comfortable with Linux command line, try Entrez Direct and GNU parallel.

First, I recommend acquiring a NCBI API key to increase the maximum requests per second.

After that, if you already have the accession numbers that you want to download, just create a tab delimited file (eg arg.tsv) contains 3 columns (no header), including accession number, sequence start and sequence stop.

To fetch the desired sequence region and save it in separate file in parallel:

```
parallel -a arg.tsv -C '\t' -j8 "efetch -db nuccore -id {1} -seq_start {2} -seq_stop {3} -format fasta >{1}.fa"

```

Or you can use Biopython `from Bio import Entrez`. But I have not tried yet.
> 0 votes

---
Tags: r, python, genome, phylogenetics, sequence-analysis
---


# Question
Title: Convert sequence in MS doc to fasta or genbank
How can I convert a sequence provided in a Microsoft doc file into a fasta or genbank format?

# Answer
The CrazyDoc Python package can convert Genbank/Fasta/Snapgene/MSDoc sequences into Biopython records, which can be saved as Genbank or Fasta.

A web interface / demo is provided at EGF CUBA (Collection of Useful Biological Apps): Convert Sequence Files.

*Disclaimer: I'm the current maintainer of CrazyDoc.*
> 2 votes

# Answer
A fasta file is just in plain text format. Pandoc can convert `.docx` format to plain text format:

```
pandoc -i input.docx -t plain >output.fasta

```
> 0 votes

---
Tags: fasta, format-conversion, genbank, egf
---


# Question
Title: What do these characters mean in a PDB file?
The 33rd atom information in ***3nir.pdb*** is written as follows:

`ATOM 33 N ATHR A 2 4.935 -11.632 15.046 0.74 2.95 N`

How should I interpret `N A` in the above line?

I mean, what does it mean? Is it Sodium or Nitrogen?

If it is Nitrogen, what does `A` mean then?

# Answer
from https://www.wwpdb.org/documentation/file-format

Protein Data Bank Contents Guide: Atomic Coordinate Entry Format Description Version 3.30 Document Published by the wwPDB This format complies with the PDB Exchange Dictionary (PDBx) http://mmcif.pdb.org/dictionaries/mmcif\_pdbx.dic/Index/index.html. I believe PDB format is kind of outdated, but not sure..

At page 177 it says about altLoc:

* AltLoc is the place holder to indicate alternate conformation. The alternate conformation can be in the entire polymer chain, or several residues or partial residue (several atoms within one residue). If an atom is provided in more than one position, then a non-blank alternate location indicator must be used for each of the atomic positions. Within a residue, all atoms that are associated with each other in a given conformation are assigned the same alternate position indicator. There are two ways of representing alternate conformation- either at atom level or at residue level (see examples).
* For atoms that are in alternate sites indicated by the alternate site indicator, sorting of atoms in the ATOM/HETATM list uses the following general rules: • In the simple case that involves a few atoms or a few residues with alternate sites, the coordinates occur one after the other in the entry. • In the case of a large heterogen groups which are disordered, the atoms for each conformer are listed together.

your question could be more useful if you post the entire description for residue 2:

```
ATOM     33  N  ATHR A   2       4.935 -11.632  15.046  0.74  2.95           N
ATOM     34  N  BTHR A   2       5.078 -11.406  15.180  0.31  2.78           N
ATOM     35  CA ATHR A   2       5.757 -11.521  13.850  0.81  3.02           C
ATOM     36  CA BTHR A   2       5.773 -11.153  13.921  0.20  2.67           C
ATOM     37  C  ATHR A   2       7.070 -10.839  14.210  0.74  2.82           C
ATOM     38  C  BTHR A   2       7.155 -10.559  14.193  0.29  1.80           C
ATOM     39  O  ATHR A   2       7.152  -9.941  15.050  0.80  3.31           O
ATOM     40  O  BTHR A   2       7.214  -9.641  15.012  0.25  2.41           O
ATOM     41  CB ATHR A   2       4.976 -10.693  12.813  0.87  5.53           C
ATOM     42  CB BTHR A   2       4.896 -10.354  12.941  0.25 12.07           C
ATOM     43  OG1ATHR A   2       4.611  -9.432  13.388  1.00  6.88           O
ATOM     44  OG1BTHR A   2       3.743 -11.083  12.501  0.25  9.57           O
ATOM     45  CG2ATHR A   2       3.858 -11.584  12.293  0.75 10.03           C
ATOM     46  CG2BTHR A   2       5.683  -9.885  11.726  0.27  5.90           C
ATOM     47  H  ATHR A   2       4.547 -10.814  15.527  0.75  3.44           H
ATOM     48  H  BTHR A   2       5.510 -10.211  15.754  0.25  2.90           H
ATOM     49  HA ATHR A   2       5.962 -12.339  13.548  0.75  3.32           H
ATOM     50  HA BTHR A   2       4.036  -9.929  13.477  0.25  2.86           H
ATOM     51  HB ATHR A   2       5.648 -10.589  11.938  0.75  5.43           H
ATOM     52  HB BTHR A   2       4.644  -9.326  13.574  0.25  5.67           H
ATOM     53  HG1ATHR A   2       5.030  -9.344  14.216  0.75  8.74           H
ATOM     54  HG1BTHR A   2       3.236 -11.198  13.399  0.25 10.21           H
ATOM     55 HG21ATHR A   2       4.096 -12.441  11.924  0.75 10.92           H
ATOM     56 HG21BTHR A   2       6.542  -9.278  12.024  0.25  9.66           H
ATOM     57 HG22ATHR A   2       3.222 -10.974  11.650  0.75 10.92           H
ATOM     58 HG22BTHR A   2       5.039  -9.142  11.179  0.25  9.66           H
ATOM     59 HG23ATHR A   2       3.163 -11.738  13.200  0.75 10.92           H
ATOM     60 HG23BTHR A   2       5.904 -10.639  11.169  0.25  9.66           H

```
> 2 votes

---
Tags: pdb
---


# Question
Title: Compressing read data (converting fastq to fastq.gz) on Windows
I downloaded some of my data in fastq format instead of fastq.gz. I want to upload my data on a server now.

I was wondering if there is a way to convert my fastq files to fastq.gz on windows (I read some stuff about python, but I'm not quite sure) before uploading them on the server, since it takes a lot of time uploading the fastq files.

# Answer
> 7 votes
`7zip` is file archiver capable of running `gzip`-style compression, and runs on Windows: https://www.7-zip.org/

# Answer
> 2 votes
If your Windows 10 build number is 17063 or later, you can use `tar` command to compress file in Windows PowerShell:

```
tar -vczf output.fq.gz input.fastq

```

The `-z` flag means to compress archive with gzip.

# Answer
> -1 votes
Bioinformatics is 99% done in Unix-like environments, and trying to come up with Windows solutions is unnecessary effort. The Windows Subsystem Linux (WSL2) provides you with a native Linux distribution within Windows that easily allow you to do any Unix-ish thing such as running gzip. It is straight-forward to install, see its online documentation.

---
Tags: fasta, fastq
---


# Question
Title: deseq2 single factor design output interpretation
This is regarding the single factor design For example if I have Age or other continuous numerical variable how to provide that into the design formula.

For this post do i need to 'You could dichotomise your continuous variables into meaningful groups' or it can go without it?

```
My simple design

dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered, colData=coldata, design= ~ Age)

dds <- estimateSizeFactors(dds)

dds <- estimateDispersions(dds)

dds  <- DESeq(dds, parallel = TRUE)

```

Here in case of metadata/coldata I m giving a single numerical value which is Blast percentage. In case of sample I have like 5 sub-types from M0 to M5.

Now for the interpretation part How do I interpret the result?

Would it be as such

the expression differences between the my subtypes due to 'Age'

I'm bit confused since in my coldata I'm not providing any information regarding my subtypes. So when I run this

```
resultsNames(dds)  I get this

[1] "Intercept"     "Age"

```

So if I would like to know if there is a difference between which subtype due to this Age variable how do i know or my question is violating the statistical design and rationale here?

To know the differences in subtype I have do which is proving the `FAB` which are basically my sub-types information where I have tested pairwise.

I would like to know if I give any numerical variable to my design how do I interpret output?

Any suggestion or help would be really appreciated.

# Answer
> 2 votes
Your design is `~ Age` so you can only test how Age affects each gene and it is assumed that the effect is linear. If you would like to test the effect of sub type on gene expression you'll have to introduce it to the design e.g. `~ Age + Sub.Type` (or whatever the column name is). This will not test how the effect of Age is different between between the subtypes, if you would like to test this you will have to include the interaction term in the design formula: `~ Age + Sub.Type + Age:Sub.Type`, this design have a lot of parameters and you will need a lot of samples to have a reliable test.

Hope it helped

---
Tags: deseq2
---


# Question
Title: How can I parse alternative atom information?
I am trying to parse PDB files.

Say, a PDB file has the following data:

```
ATOM     33  N  ATHR A   2       4.935 -11.632  15.046  0.74  2.95           N
ATOM     34  N  BTHR A   2       5.078 -11.406  15.180  0.31  2.78           N
ATOM     35  CA ATHR A   2       5.757 -11.521  13.850  0.81  3.02           C
ATOM     36  CA BTHR A   2       5.773 -11.153  13.921  0.20  2.67           C
ATOM     37  C  ATHR A   2       7.070 -10.839  14.210  0.74  2.82           C
ATOM     38  C  BTHR A   2       7.155 -10.559  14.193  0.29  1.80           C
ATOM     39  O  ATHR A   2       7.152  -9.941  15.050  0.80  3.31           O
ATOM     40  O  BTHR A   2       7.214  -9.641  15.012  0.25  2.41           O
ATOM     41  CB ATHR A   2       4.976 -10.693  12.813  0.87  5.53           C
ATOM     42  CB BTHR A   2       4.896 -10.354  12.941  0.25 12.07           C
ATOM     43  OG1ATHR A   2       4.611  -9.432  13.388  1.00  6.88           O
ATOM     44  OG1BTHR A   2       3.743 -11.083  12.501  0.25  9.57           O
ATOM     45  CG2ATHR A   2       3.858 -11.584  12.293  0.75 10.03           C
ATOM     46  CG2BTHR A   2       5.683  -9.885  11.726  0.27  5.90           C
ATOM     47  H  ATHR A   2       4.547 -10.814  15.527  0.75  3.44           H
ATOM     48  H  BTHR A   2       5.510 -10.211  15.754  0.25  2.90           H
ATOM     49  HA ATHR A   2       5.962 -12.339  13.548  0.75  3.32           H
ATOM     50  HA BTHR A   2       4.036  -9.929  13.477  0.25  2.86           H
ATOM     51  HB ATHR A   2       5.648 -10.589  11.938  0.75  5.43           H
ATOM     52  HB BTHR A   2       4.644  -9.326  13.574  0.25  5.67           H
ATOM     53  HG1ATHR A   2       5.030  -9.344  14.216  0.75  8.74           H
ATOM     54  HG1BTHR A   2       3.236 -11.198  13.399  0.25 10.21           H
ATOM     55 HG21ATHR A   2       4.096 -12.441  11.924  0.75 10.92           H
ATOM     56 HG21BTHR A   2       6.542  -9.278  12.024  0.25  9.66           H
ATOM     57 HG22ATHR A   2       3.222 -10.974  11.650  0.75 10.92           H
ATOM     58 HG22BTHR A   2       5.039  -9.142  11.179  0.25  9.66           H
ATOM     59 HG23ATHR A   2       3.163 -11.738  13.200  0.75 10.92           H
ATOM     60 HG23BTHR A   2       5.904 -10.639  11.169  0.25  9.66           H

```

We see that there are many alternative atoms in the 2nd residue.

How should I choose the atoms?

Should I just randomly choose one atom from the alternatives? How can my parser learn to differentiate between those two atoms (say, 33 and 34)? Coz, from the parsing point of view, there is no indication that they are related.

How can I parse alternative atom information in a PDB file?

# Answer
> 4 votes
Generally alternate location are used in structures with partial occupancy.

In https://www.wwpdb.org/documentation/file-format-content/format33/sect9.html you will see that the fourth numerical column is occupancy. So you could take the higher occupancy version. The tool PDBfixer can do it for you: github.com/openmm/pdbfixer

* Do not assume the sum of occupancy is 1 —it should but errors do creep up especially for close alternate location that overlap.
* The anisotropic and isotropic B-factors will be lower than were the structure refined with a single conformation (like normal atoms).
* Same letters in different residues should not be assumed to be paired up. Models as seen for NMR structures are a solution for that.

Lastly... Your structure has protons, which is unusual (exceptions exist). It is fine if it is a neutron diffraction structure or NMR. If you are not expecting protons (X-ray structure), skip it out of safety.

---
Tags: proteins, pdb
---


# Question
Title: How to plot a 2D RNA structure in python using dot-bracket input?
I need some help with forgi library for visualisation of RNA secondary structure.

I have an RNA sequence and dot-bracket notation of the secondary structure and I would like to plot 2D graph like that you can find here.

The input file is a `.fx` file, or a fasta file where the sequence is single-line and the third line is the dot-bracket notation of the RNA structure.

```
>1y26
CGCUUCAUAUAAUCCUAAUGAUAUGGUUUGGGAGUUUCUACCAAGAGCCUUAAACUCUUGAUUAUGAAGUG
((((((((((..(((((((.......)))))))......).((((((.......))))))..)))))))))

```

If there any other possibilities to plot it in real-time (for instance, in jupyter notebook avoiding any crossing external scripts), I will be very thankful to learn them!

# Answer
> 1 votes
I wrote a function to do that a while ago, maybe it helps. It only works with dot-bracket notation, and it turns it into ordered pairs which then you can easily plot.

```
import pandas as pd
import seaborn as sns

def dot_bracket_to_pairs(ss_string):
    '''
    Takes dot and bracket string, returns dataframe
    with paired bases.
    If any invalid characters are in the structure, it
    will interpret them as dots, as it only reads parentheses.
    '''
    index_list = []
    pairs = {}

    for index, char in enumerate(ss_string):
        if char == '(':
            index_list.append(index)
        if char == ')':
            try:
                # pair to last item in the list in dictionary
                pairs[index_list.pop()] = index
            except IndexError:
                print(f'Invalid structure, found extra \')\' in position {index}')

    if len(index_list) != 0:
        for item in index_list:
            print(f'Invalid structure, found extra in \'(\' in position {item}')

    df_pairs = pd.DataFrame(pairs.items(), columns=['Y', 'X'])

    return df_pairs

db_string = '((((((((((..(((((((.......)))))))......).((((((.......))))))..)))))))))'
df = dot_bracket_to_pairs(db_string)

sns.scatterplot(data=df, x="X", y="Y")

```

# Answer
> 1 votes
One way simply is to use VARNA, here is a link ... https://varna.lri.fr. This is well known for its power to view secondary structure (e.g. its in Jalview). It is generally considered the best viewer app for secondary structure.

However, the solution by @nquinoneso whilst simple, is preferable. The output is here:

The y-axis is the 3' and the x-axis is the 5', the dots represent hairpins. Whilst its simple, the power of this approach is the automation by direct access to the (pandas) dataframe.

*I do recommend upvoting the post.*

Complex secondary structure would not work because it will confuse the 5' stems with the immediate 3' stems which don't hybridise and a more advanced solution using nucleotide code would be needed. Thats the bit for further development, using the index both for the dot-bracket and the nucleotide sequence to confirm the hairpin. VARNA does perform this calculation (because it identifies canonical mutations)... but you can't access the underlying array/dataframe which is certainly a limitation.

---
Tags: python, visualization, secondary-structure
---


title: Caching and parallelization

text:
# Question
Title: Caching and parallelization
Currently algorithmically developing a calculation for "co-dependent" mutation (quotation appropriate). The algorithm uses chunked biopython alignment objects which are parallelized.

**Aim** One optization is caching a database within parallelization, but I ain't sure whether the parallelization duplicates the cache.

**Method** Each alignment chunk (biopython object) is passed onto a separate CPU and refers back to the original database to obtain its metadata in the final stage of the calculation (dragging metadata around this algorithm makes no sense).

The idea is to cache the database to prevent pulling it in from the harddrive:

```
from multiprocessor import Pool
from myOOPlibrary import mystuff

```

Facade design, the parallelised index to the alignment chunks enter into the myOOPlibrary:

```
from cachetools import cached, TTLCache

cache = TTLCache(maxsize=1, ttl=500)

Class mystuff (stuff2):
... # init omitted

    def mymethod(self, ..):
    # to be updated

Class stuff2 ():
... # init omitted

    @cached(cache)
    def getDB(self, ...):
.... #  pull it in (harddrive) and parse to a dictionary
    return metadataDB

```

**Question** I'm wondering if caching is the right approach in parallelization, specifically whether this will end up with multiplying the database held in RAM...

**Alternatives** The alternativees to caching are:

* IO to JSON (not cool)
* IO via pickle (binary IO): a JSON will readily`pickle`

The possible advantage of the alternativees are the database is written once and accessed by all cores for certain.

**Underlying question**

* When caching is the database cached for each core? If so the RAM overhead becomes an issue because the database is duplicated `total CPUs - 1` times
* Alternatively is the cache shared between CPUs centrally and each CPU refers to this location?

# Answer
***Summarising***

1. @Jesse's comment I think is the right answer, open the database outside parallelisation. What I didn't thoroughly describe (and what was essentially requested) is the architecture, the alignments entering into an entire OOP library and the whole library is parallelised, i.e. every new alignment entering library is on its own CPU. The answer is obvious: instantiate the OOP library with the entire database and thats a key purpose of the design.
2. The technical answer is `joblib`. This is built to cache across parallelised archetechture as the `joblib` team point out **here**. The criticism is that there's parallelisation within parallelisation. The discussion\* is as much in the comments as in the response.
3. @user172818 response is correct too, the algorithm would benefit from C. Its a lot of code to rewrite in `C`, but sidestepping via Cython is a future development (its a hybrid language connecting C and Python).

\*, https://stackoverflow.com/questions/25033631/multiple-processes-sharing-a-single-joblib-cache

---

There was a comment hinting this is Stack Overflow, that's difficult. Ultimately the technical answer (which is Stack Overflow) isn't the correct answer as far as I understand OOP.

---

***Feedback***

Just to wrap this, point 1. worked fine, thank you. It means a database sits on each processor in parallelisation, but is imported once. In truth, it worked nicely and provided cleaner code. Re-writing to instantiate once would affect other stuff so would be more costly than a cost-saving.

Python can be very fast (e.g. numpy), but also really slow (hence it has critiques), extensively parallelising helps bridge the gap. Python parallelisation has critiques which I strongly, strongly believe are unfounded: I parallelise loads via `multiprocessor` library no problems.

Caching here, i.e. parallelising across a parallelisation could stretch things too far.

***Last point***

Point 3, re-writing in Cython/`C`, is a greater priority than the testing required for parallelised caching: can't parallelise all data and that code is better as `C`.
> 1 votes

---
Tags: python, covid-19, algorithms, mutations
---


# Question
Title: What tool I can use to map short-reads sequences to reference genome and get specific mapped size
I have about 300 90-bp sequences which I would like to map to a reference genome to make it longer to a 300bp sequence, wherein the 90bp is at the middle.

Anyone knows what bioinformatics tool I can use for this? Thank you!

*Actual scenario:* The SNP sequences are from a SNP panel data (90bp sequence, SNP is at the middle). We would like to extend the sequence to 300bp, retaining the SNP at the middle so we can make a sequence capture panel of the SNP.

# Answer
If you have a SNP panel then you have the coordinates on the genome. You can use tools such as bedtools slop to expand the coordinates to the desired length and then extract those regions from the target genome using bedtools getfasta.
> 3 votes

---
Tags: sequence-alignment, snp, read-mapping
---


# Question
Title: Find enriched processes in modules obtained by WGCNA
I am doing a gene enrichment of different clusters obtained by WGCNA, but there are some of them in which I don't get enriched processes, is there any way to find processes to these modules, for example to associate them to orthologues of other organisms. If there are different ways, could you send me links on how to do it?

# Answer
Well, not all modules will have significant enrichment in gene sets such as GO terms or KEGG pathways. In my experience, the smaller modules are often not significantly enriched. They could be the result of some technical variation, they could contain genes localized to a particular genomic region (do you also calculate enrichment in genomic regions?) or they could be a particular gene family. WGCNA modules often are enriched in known pathways and processes, but that does not mean that all will be.
> 2 votes

---
Tags: r, rna-seq, wgcna
---


# Question
Title: Which atoms are not found in protein PDB files?
I am developing an educational Bioinformatics framework, I need to know

* Which atoms are absolutely not found in any PDB files?

in the following list -

# Answer
## PDB file readability

It is a question of implementation. For example:

```
from rdkit import Chem

mol = Chem.MolFromSmiles('[Og]')
pdbblock = Chem.MolToPDBBlock(mol)

import pymol2

with pymol2.PyMOL() as pymol:
    pymol.cmd.read_pdbstr(pdbblock, 'Oganesson')

```

shows that both RDKit and PyMol can deal with oganesson. However, this is just handling the atoms in coordinate space. But as soon as we touch it in molecular mechanics way things blow up:

```
from rdkit.Chem import AllChem

AllChem.EmbedMolecule(mol)  # UFFTYPER: Unrecognized atom type: Og (0)

```

But generally most molecular mechanics tools struggle with coordinated iron or similar. Take the mangled case of https://www.rcsb.org/ligand/SIR which ought to have been planar...

## PDB DB

In the PDB polymeric units and ligands are called chemical components. They have a 3 letter code. The list is easily parsable. The problem is that some lists are incomplete: https://ftp.ebi.ac.uk/pub/databases/pdb/data/component-models/complete/chem\_comp\_model.cif.gz when read

```
from pdbecif.mmcif_tools import MMCIF2Dict
models: dict = MMCIF2Dict().parse('chem_comp_model.cif')

```

has a large number of entries missing because of \<reasons\>.

So using a more crude approach and simply querying component by component:

```
with open('chem_comp.list') as fh:  # this is from somewhere else in the ftp
    chem_resns = *map(str.strip, fh),

import requests
import re
elements = set()
for i in range(len(chem_resns)//100):
    response = requests.post('https://www.ebi.ac.uk/pdbe/api/pdb/compound/summary/', ','.join(list(mia)[100*i:100*(i+1)]))
    response.raise_for_status()
    for resn in response.json().values():
        for entry in resn:
            if entry['formula'] is None: # the unknown ligand/ions etc.
                continue
            elements.update(re.findall('[^\d ]+', entry['formula']))

```

Gives the following elements that are present:

```
{'Ag',
 'Al',
 'Am',
 'Ar',
 'As',
 'Au',
 'B',
 'Ba',
 'Be',
 'Bi',
 'Br',
 'C',
 'Ca',
 'Cd',
 'Ce',
 'Cf',
 'Cl',
 'Cm',
 'Co',
 'Cr',
 'Cs',
 'Cu',
 'D',
 'Dy',
 'Er',
 'Eu',
 'F',
 'Fe',
 'Ga',
 'Gd',
 'H',
 'Hf',
 'Hg',
 'Ho',
 'I',
 'In',
 'Ir',
 'K',
 'Kr',
 'La',
 'Li',
 'Lu',
 'Mg',
 'Mn',
 'Mo',
 'N',
 'Na',
 'Nd',
 'Ni',
 'O',
 'Os',
 'P',
 'Pb',
 'Pd',
 'Pr',
 'Pt',
 'Pu',
 'Rb',
 'Re',
 'Rh',
 'Ru',
 'S',
 'Sb',
 'Sc',
 'Se',
 'Si',
 'Sm',
 'Sn',
 'Sr',
 'Ta',
 'Tb',
 'Te',
 'Th',
 'Ti',
 'Tl',
 'U',
 'V',
 'W',
 'X',
 'Xe',
 'Y',
 'Yb',
 'Zn',
 'Zr'}

```

While those that are absent are:

```
pt = Chem.GetPeriodicTable()
full = set(map(pt.GetElementSymbol, range(1, 118+1)))
missing = full - elements
for m in sorted(missing, key=pt.GetAtomicNumber):
    print(pt.GetAtomicNumber(m), m)

```

```
2 He
10 Ne
32 Ge
41 Nb
43 Tc
61 Pm
69 Tm
84 Po
85 At
86 Rn
87 Fr
88 Ra
89 Ac
91 Pa
93 Np
97 Bk
99 Es
100 Fm
101 Md
102 No
103 Lr
104 Rf
105 Db
106 Sg
107 Bh
108 Hs
109 Mt
110 Ds
111 Rg
112 Cn
113 Nh
114 Fl
115 Mc
116 Lv
117 Ts
118 Og

```
> 4 votes

---
Tags: pdb
---


# Question
Title: BiocManager::install("DiffBind") mergeOne.c:341:3: note: use option -std=c99 or -std=gnu99 to compile your code
I am using

```
BiocManager::install("DiffBind")

```

to install the package "DiffBind" on CentOS\_7, I got the error:

```
gcc -I"/opt/common/CentOS_7/R/R-4.0.5/lib64/R/include" -DNDEBUG -D_FILE_OFFSET_BITS=64 -I'/lila/home/lid/R/x86_64-pc-linux-gnu-library/4.0/Rhtslib/include' -I'/admin/opt/common/CentOS_7/R/R-4.0.5/lib64/R/library/Rcpp/include' -I/usr/local/include   -fpic  -g -O2  -c mergeOne.c -o mergeOne.o
mergeOne.c: In function ‘mungeTarget’:
mergeOne.c:341:3: error: ‘for’ loop initial declarations are only allowed in C99 mode
   for (int i=l;i<l+r;i++) {
   ^
mergeOne.c:341:3: note: use option -std=c99 or -std=gnu99 to compile your code
make: *** [mergeOne.o] Error 1
ERROR: compilation failed for package ‘DiffBind’

```

Can anyone tell me what should I do? Thanks

# Answer
> 1 votes
You can set a specific compiler in `~/.R/Makevars` file, create the folder and write it down in Makevers:

```
CC = gcc -std=c99

```

or:

```
CC = gcc -std=gnu99

```

# Answer
> 3 votes
The previous answer is correct, in that the compiler wants to be able to use `C99` to compile your library.

However, the solution of putting `CC = gcc -std=c99` into your `~/.R/Makevars` is a little bit blunt, since **everything** that you now compile not directly from source is going to use the `-std=c99` flag. This is a recipe for some confusing and hard to track down error messages. It's a bit like installing everything in your base conda environment and waiting for everything to break.

It's a safer option to use `withr`, which allows you to specify compiler flags for individual installations, keeping everything nice and neat. It works a bit like the `with` (or `within`) function in R which constructs an environment to execute some command in.

```
library(withr)

with_makevars(
    c(PKG_CFLAGS = "-std=c99"),
    BiocManager::install("DiffBind"),
    assignment = "+="
)

```

It can also be performed directly with `BiocManager`:

```
BiocManager::install(
    "DiffBind",
    configure.vars = c("CC" = "gcc", "-std" = "c99")
)

```

---
Tags: r, bioconductor
---


# Question
Title: How good does the assembly of an NCBI prokaryotic genome have to be in order to argue gene loss?
NCBI has several labels for assembly completeness - Complete, Scaffold, Chromosome and Contig. Complete would be a circularized genome (or linear, rarely)

For a Complete genome it's fairly straightforward - if not present, a gene can safely assumed to be missing (or possibly present on plasmids, if those aren't available in the database).

But what about the others? With Contig-level genomes, I'd feel uncertain about making any statements (especially since some of these projects are split into thousands of contigs). But what about Scaffold, or Chromosome? Are they expected to be complete in terms of genes, but just missing repetitive sequences?

# Answer
I would suggest that in the absence of "Complete" no very confident assertion can be made about gene presence/absence in a genome assembly. You need to actually have the chromosome sequence to make strong statements.

I think that most people (me included) would reluctantly accept that genomes with very good BUSCO or CheckM numbers (depending on your organism's taxonomy) are probably more or less accurate in terms of presence/absence. For metagenomics, I would extend that to the MIMAG criteria, which are pretty stringent.

These criteria are of course independent of the NCBI categories, which are pretty coarse-grained.

That said, I think that you can make stronger arguments about gene presence/absence from the sequencing reads from which the genome assembly is derived. Find all the reads for your genome, map them to the genome assembly and also to the closest related ortholog of the gene in question. How well is the gene covered? How do the alignments look? Harder for metagenomics, but there are still probably tricks one could use.
> 2 votes

---
Tags: gene, assembly, phylogenetics
---


# Question
Title: Why do molecular generation models maximize “penalized logP” as a measure of drug-likeliness?
I found that Lipinski's rule of five states that Log P (octanol-water partition coefficient, lipophilicity measure) usually should not exceed 5.

Many papers about drug discovery machine learning models tell about maximization of "penalized logP", but the paper they refer in the end does not contain any information on it. Let me show that.

Let's check , here we read about penalized log P:

> The second task is to produce novel molecules with desired properties. Following (Kusner et al., 2017), our target chemical property y(·)is octanol-water partition coefficients (logP) penalized by the synthetic accessibility (SA)score and number of long cycles
>
> ```
> y(m) = logP(m)−SA(m)−cycle(m)
>
> ```
>
> where cycle(m)counts the number of rings that have more than six atoms.

Then we go into Kusner et al., 2017:

> For the second optimization problem, we follow (Gómez-Bombarelli et al., 2016b) and optimize the drug properties of molecules. Our goal is to maximize the water-octanol partition coefficient (logP), an important metric in drug design that characterizes the drug-likeness of a molecule.As in Gómez-Bombarelli et al. (2016b) we consider a penalized logP score that takes into account other molecular properties such as ring size and synthetic accessibility (Ertl & Schuffenhauer, 2009).

Here is Gómez-Bombarelli, 2016 and it seems to contain nothing on the topic (I can be mistaken, please point me then), except

> The objective we chose to optimize was 5×QED−SAS, where QED is the Quantitative Estimation of Drug-likeness (QED),37 and SAS is the Synthetic Accessibility score.36Thisobjective represents a rough estimate of finding the most drug-like molecule that is also easy to synthesize.

Other examples of the papers also leading to the same papers are , , etc.

What does penalized logP mean for drug-likeliness and why it should be maximized?

# Answer
> 2 votes
**A positive logP is good (hydrophobic)**, but high results in ADMET problems (cytochromes break overly hydrophobic compounds) and aggregation (gunk, tar, crud are words happily thrown around). However, in most cases the logP is not greater than +3 because good binding is result of hydrogen bonds etc —some bad algorithms will suggest a polyaromatic hydrocarbon that makes lots of pi-pi interactions but these are few.

Confusingly, **a high synthetic accessibility score is bad (hard)**. 3 is the average and 10 is basically something impossible. This score is how hard it is to make (i.e. pricier), which is one of the major drivers.

The number of cycles is number of rings. Even if you get pi-pi, cation-pi and sulfur-pi bonds you also get aggregation due to stacking (a crud). So is somewhat bad.

So basically that score is

```
high=good minus high=bad minus high=bad

```

If difference in relative Gibbs energy of binding is thrown it, negative ∆∆G is good, so that would be subtracted. As its atom count dependent ligand efficiency is better, but that ranges from 0 to -1.2 kcal/mol. So would need scaling.

# Answer
> 1 votes
You were trapped by arXiv versioning convention that shows you the **latest** version of the preprint as default. It seems that this experiment was removed from the final publication.

If you have a look at the v1 (https://arxiv.org/abs/1610.02415v1) instead of the latest version (v3), you find the explanation:

> As a simple example, we first attempt to maximize the water-octanol partition coefficient $logP$, as estimated by RDkit. \[43\] $logP$ is an important element in characterizing the druglikeness of a molecule, and is of interest in drug design. To ensure that the resulting molecules to be easy to synthesize in practice, we also incorporate the synthetic accessibility \[44\] (SA) score into our objective. Our initial experiments, optimizing only the $logP$ and SA scores, produced novel molecules, but ones having unrealistically large rings of carbon atoms. To avoid this problem, we added a penalty for having carbon rings of size larger than 6 to our objective. Thus our preferred objective is, for a given molecule m, given by: $J(m) = logP(m) − SA(m) − ring-penalty(m)$, (1)
>
> where the scores $logP(m)$, $SA(m)$, and $ring-penalty(m)$ are normalized to have zero mean.

---
Tags: drugs, small-molecules
---


# Question
Title: Generating contigs from illumina Paired End Reads (Metagenome)
I'm a wet-lab biologist desperately trying to get to grips with some Bioinformatics.

I have some shotgun Illumina paired-end data from a gut microbiome project I'm working on. I'm considering how to go about generating contigs from reads. I've been recommended to use SPAdes which I'm reading up on at the moment.

Originally I was thinking of using FLASh however I was told it might not be suitable for metagenomic data.

**My question simply is, what makes certain programs suitable for merging paired-end reads from a genome but not a metagenome?** I find myself a little spoilt for choice (eg. FLASh/SPAdes/BBmap... ) with bioinformatics tools I'm always concerned I'm picking the wrong one.

Should I ideally be running my data through a couple of tools and comparing them?

# Answer
NB: "merging PE reads" (such as e.g. FLASH does) is not the same as "sequence assembly from PE reads" (what MetaSPAdes does). The latter generates contigs, the former generates (some) merged shotgun library inserts. I think you want sequence assembly, i.e. a tool like MetaSPAdes.

**Don't be tricked by the loose language of "merging" that is sometimes applied to the sequence assembly problem. It is much more complicated than that.**

If you want to learn more about sequence assembly, I suggest you read the sequence assembly wikipedia page.

I would probably recommend MegaHit over metaSPAdes, because it requires much less memory / runtime, and lower frequency of errors.

But generally speaking, both programs are building different flavors of de Bruijn Graphs from metagenomic reads to generate contigs/scaffolds. This is the accepted method by this point for Illumina metagenomic data, because it is the method that is capable of handling the amount of small pieces of data in a tractable fashion. Alternatives such as overlap-layout-consensus graphs or string graphs just don't work as well for the metagenomic problem.

You are probably looking for methods that talk about complicated, computationally complex assembly algorithms and graphical methods for summarizing lots of short reads. You are not looking for simple read processing/alignment tools like BBMap/BBTools or FLASH. **To reiterate, simple read merging tools will not give you assembled contigs.**
> 2 votes

---
Tags: metagenome
---


# Question
Title: How can I obtain FTP links to studies in ENA given a BioProject ID?
How can I programmatically obtain ftp links to RNA seq fastq files in ENA? Here's an example of a link that I would be interested in obtaining:

`ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_1.fastq.gz`

In particular, is there some tool that, given the BioProject ID (here, PRJNA506829), would be able give me all ftp links for the runs in the project, or would I need to write a web scraper to do it?

# Answer
Here are some command lines I used for that purpose in bash. Simply prepare a text file containing each accession number (SRR/ERR) you want and create a for loop. Here I used prozilla to speed up downloads but you may use `wget` either.

for index in $(cat list_of_accessions) ; do

if [ ${#index} -eq 9 ]; then
            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index}/${index}_1.fastq.gz
            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index}/${index}_2.fastq.gz
            else
            if [ ${#index} -eq 10 ]; then
                    proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/00${index:9:9}/${index}/${index}_1.fastq.gz
                            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/00${index:9:9}/${index}/${index}_2.fastq.gz
                    else
                    if [ ${#index} -eq 11 ]; then
                     proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/0${index:9:10}/${index}/${index}_1.fastq.gz
                                        proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/0${index:9:10}/${index}/${index}_2.fastq.gz
                    else
                        proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index:9:11}/${index}/${index}_1.fastq.gz
                                                proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index:9:11}/${index}/${index}_2.fastq.gz
                        fi
                fi

            fi
done
```
</code></pre>
> 2 votes

# Answer
They have an API you can interact with.

If you need to get files for only a few different projects: Search for your accession ID in the browser (leading to https://www.ebi.ac.uk/ena/browser/view/PRJNA506829). Filter the show selected columns to only fastq\_ftp, click download tsv to get the list of ftp links.
> 1 votes

# Answer
pysradb can fetch ENA/SRA fastq/bam links (if available):

```
$ pysradb metadata SRR8240860 --detailed
run_accession study_accession experiment_accession experiment_title                                               experiment_desc                                                organism_taxid  organism_name           library_strategy library_source  library_selection library_layout sample_accession sample_title instrument   total_spots total_size  run_total_spots run_total_bases run_alias      sra_url_alt1                                    sra_url_alt2                                    sra_url                                                                                experiment_alias source_name strain/genotype developmental stage ena_fastq_http ena_fastq_http_1                                                                 ena_fastq_http_2                                                                 ena_fastq_ftp ena_fastq_ftp_1                                                                     ena_fastq_ftp_2
SRR8240860    SRP170618       SRX5059122           GSM3487689: fer-15 Day1 rep1; Caenorhabditis elegans; RNA-Seq  GSM3487689: fer-15 Day1 rep1; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              PAIRED         SRS4075216       N/A          HiSeq X Ten  40494415    4755049874  40494415        12148324500     GSM3487689_r1  gs://sra-pub-src-3/SRR8240860/RRA0719_R2.fq.gz  s3://sra-pub-src-3/SRR8240860/RRA0719_R2.fq.gz  https://sra-downloadb.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR8240860/SRR8240860.1  GSM3487689       whole worm  fer-15(b26ts)   Adult day 1         N/A            http://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_1.fastq.gz  http://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_2.fastq.gz  N/A           era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR824/000/SRR8240860/SRR8240860_1.fastq.gz  era-fasp@fasp.sra.ebi.ac.uk:vol1/fastq/SRR824/000/SRR8240860/SRR8240860_2.fastq.gz

```
> 0 votes

# Answer
You can do this in one line using NCBI's EntrezDirect e-utils:

`esearch -db sra -query PRJNA506829 | efetch -format runinfo`

The output will be in csv, but you can format it to tab using csvkit's `csvformat -T`
> 0 votes

---
Tags: rna-seq, public-databases, data-retrieval
---


title: Converting genetic matrix to STRUCTURE

text:
# Question
Title: Converting genetic matrix to STRUCTURE
I want to convert a genetic matrix (rows are samples, columns are alleles) into a STRUCTURE file (output .str) within R. Is there a preferred way to achieve this? I've found genind\_to\_structure in the graph4lg R package, but I wanted to see if there are other approaches.

An example of what I mean by a "genetic matrix":

```
# Build a matrix of alleles
temp <- lapply(1:30, function(i) sample(1:9, 4, replace = TRUE))
temp <- sapply(temp, paste, collapse="")
temp <- matrix(temp, nrow = 10, dimnames = list(paste("ind",1:10), paste("loc",1:3)))
temp

# Convert matrix to genind format (for adegenet package)
test <- df2genind(temp, ploidy = 4, sep = "")
test

# How to convert to STRUCTURE?

```

If anyone has input, please chime in! Thanks

# Answer
I ultimately found this function, which I used for the purpose of converting a genind to a STRUCTURE file. The graph4lg genind\_2\_structure command would regularly error on me (present working directory couldn't be passed as the `output` argument).
> 1 votes

---
Tags: r, file-formats, phylogenetics
---


# Question
Title: Extracting clusters, hierarchical clustering - bulk RNASeq
I have been attempting to extract particular clusters, given by hierarchical clustering outputs from the `pheatmap` function (in R).

Please find the `dput` output of the dataset in question below:

```
structure(list(`72v24.logFC` = c(-1.47511280138092, -1.10386310802134,
-1.14476608274903, -1.63563574551239, -0.895377907265352, -1.35216488439027,
-2.04949381135359, -1.0344565728341, -0.967008696515716, -1.11233305998901,
-2.02663681865178, -1.74756562789736, -1.02076838668592, -1.27981009782201,
-1.13181804725663, -1.19537930856639, -2.23476188503608, -1.67500131867632,
-0.950059974858215, -2.3124598902615, -1.93757457440581, -3.24602437386808,
-0.839869104053795, -1.40329424193006, -1.63484798229567, -0.912148120310492,
-1.33467209346367, -1.27031027753064, -0.925498608058671, -0.802956638081314,
-0.85457796397113, -1.0591348681193, -0.825278459132589, -0.671570119524199,
-1.8439768319437, -0.947033608766851, -1.63849733567277, -1.40584263305038,
-1.1317958487115, -1.40724695083764, -1.21842141154588, -2.29248800050516,
-1.41699787526721, -0.879218050664484, -2.25064962332145, 1.4685855647974,
-2.02553848518381, 1.32808748053321, -2.46062484389627, -2.45830217362141,
-0.880536960021976, -0.689060222340137, -0.863262774594, -1.64961312646758,
1.10572346456398, -0.891760767302047, -1.22716047163938, -0.680697741299578,
-1.24551344420046, -0.891918251648935, -1.64722787848732, -0.832845262108387,
-0.937129577660045, -0.891684122691785, 1.96452911491088, -1.59420664438654,
-0.921295293242016, -1.8304480291305, -0.632859040349798, 1.69982819847455,
-0.982668795391973, -1.61816842714933, -1.23602172818196, -2.15356245206668,
-1.67961932563564, 1.15478114880933, 1.32665444151868, -1.04590896017973,
-2.29770026795487, -1.27498150009975, -1.36310217581123, -1.09169550738617,
1.5349413586602, -1.66842984591719, -0.760122936706606, -1.10054393155051,
-1.00344072779106, -0.993259123487476, -1.43136824963825, -1.07214189069438,
-0.990813303020613, -1.51324812740189, -1.23465012540084, -1.00098231964124,
-1.92195764952845, -1.00139179559519, -0.940118366042393, -2.35969198434574,
-0.574626090826994, -0.884870428965888, -0.87601460107501, 1.61336157755029,
-0.777670642457121, -0.7312214522498, -0.812926777785119, -1.13443421373561,
-1.56498420350444, -0.69126972981765, -1.44351026281231, -1.21962275860236,
-2.21155555407497, -1.02015250027788, -0.604620179871691, -1.53013572266942,
-0.790050747397837, -1.17950736026735, -0.781075150161975, -0.888554144538286,
3.59172865862026, -0.8665071699621, -2.97701792242157, -1.44424941630738,
-1.50528882788757, -0.772465957539335, -1.57664185213987, -1.02371443439054,
-0.931807409493221, 1.77728969821973, -1.05455101137715, -1.09765136526264,
-0.881436416213335, -0.985605985458853, -1.31644334991692, 2.26917914828263,
-1.43844640540614, -1.75106843701475, -1.22483991752173, -1.65032906954133,
-0.686319914242903, 2.45473055753775, -1.15594379837188, 1.46397122502826,
-0.91260428248088, -0.726540148411507, -1.18227838019011, -0.940722690371935,
-1.28074043595918, -0.629033607437977, 1.01053893402133, -0.735590704487842,
-2.0620695625823, -0.692251638781327, -0.717818968339447, -0.564023101160995,
1.51541391124355, -1.89444914881727, -0.737931396419173, -1.17763728235316,
-0.887605568116903, -2.8392997898765, -1.48018352596355, 1.60698486226092,
-1.96804216410227, -0.881437183709851, 3.08730764755417, -0.600282670155725,
-0.905719394102218, -1.05553935946213, -1.13508709782814, -0.613381276596874,
-0.642647139688194, -0.861929974132332, -1.30652597213552, -1.73647794867498,
-0.619909931049252, -0.880145154869759, -0.60084263914673, -1.16825573255754,
-0.942919502394457, -1.3709930239358, -0.80595477526712, -0.628115199442632,
-1.11502998682131, -1.01274335471793, -1.16595378001606, -3.36098944183559,
-1.19617345920841, -2.10811718229285, 3.33752122936, -1.67696701892834,
-1.90639199489612, -0.574737209544065, -1.66853541929234, -0.595093777485911,
-0.7964813653225, -0.791102813817634, -0.601987407016077, -1.52559881980829,
1.15480324851026, -0.830167728660421), `72v0.logFC` = c(-5.10937036648522,
-8.41670756910048, -8.18671783299995, -2.79415692883273, -3.62090016917187,
-8.37523504570539, -3.79100494829738, -7.20659380156151, -5.46135634576091,
3.43493215839382, -1.67650075643956, -2.99391647274475, -5.77357941089076,
-5.22675557853878, 3.26436094268851, -5.98693257837598, -5.85864442290848,
-2.39995052889286, -5.19221784614079, -4.72634343497625, -7.46407809296449,
-3.04588366694082, -3.85014181238163, 4.44756662720351, -4.42031795110716,
-2.58478993868624, -5.46121416816988, -2.09364835701952, -5.42397023596482,
-5.26192471384563, -4.48979100280513, -5.14171942073048, -5.5267515390142,
3.88762033289334, -5.88489034628499, -2.85490179670763, -5.31586479048071,
-2.64877712735579, -6.24766069294439, -5.35826818217055, -4.97204642950896,
-2.82989963758335, -4.75016577988789, -3.51684070536132, -3.20931950649022,
-3.62394584217955, -2.81449967053025, -5.68308664428684, -3.11500402887324,
-4.08965188145047, -5.43176189260205, -4.30737780091769, 2.52666133143524,
-3.88213803962076, -5.517506523237, -3.27502948978619, -5.37724717693793,
2.77090845859761, -4.22554364178978, -5.81944792042417, -2.23700806354102,
-1.81804660477942, -3.45613077773548, -4.88444075082709, -5.48877199958173,
-4.90077849970695, -3.79304730634383, -0.953597409660203, -5.573920008042,
-1.83250444732307, -1.73248469438834, -2.35460087734472, -4.43029236446579,
-3.37654864236982, -3.63410171518379, -2.78276966423079, -2.01664452557169,
-3.54834187597395, -2.04395822264753, -2.15693095907293, -4.04759665569609,
-4.87928960250415, 1.1358557257212, -4.89657506999093, -2.42698466771808,
-4.55703295523896, -2.53230009478764, -1.74571568381108, 1.93875826054123,
-2.29178704384166, -1.88631492984573, -4.17811684330152, -4.55155213716622,
1.93093495783903, -3.66429213603731, 2.73817157445051, -3.96977406502995,
-4.18487497756956, 2.17512996773018, -4.41168757499049, -2.96267804428878,
-4.60683474159913, -3.0573109426718, -1.59500391193586, -4.65882130908176,
-2.68903554919404, -4.44676809434802, 1.60275224879078, -3.83760729167497,
-3.09907199540008, -2.37905633686461, -4.29650784478684, -1.88615561094275,
-4.44026540480975, -2.01996858091239, -2.46842498723876, -4.24053886261736,
-4.32561109033911, -2.96097578889787, -3.76708395833968, -2.45229520135946,
-2.33856203921674, -4.23241140779399, -3.48627475173498, 2.25888018977625,
-3.85271202182208, -4.61317704601435, -1.95186634348091, 1.3205477461731,
-2.51561846705872, -1.7031812569795, -4.44717736673796, -3.88478215160905,
-3.27547652595564, -2.81690169904308, -4.83981492357284, -2.4486942125238,
-4.31360736969916, -4.22122930700826, -3.49045260391997, -2.87862947524065,
-4.14268228866617, -2.62252967420329, -4.23902905919912, -4.14777621540367,
-2.86056694414781, -1.70668584273452, 1.3624660569909, -2.59153842674144,
-4.12393693326418, -2.12959469550212, 1.46923996462969, -1.98633421258236,
3.22176343046404, -3.52110060349997, -2.69900299957347, -4.20126962458439,
-4.68116989520068, 1.57592717464553, -3.93999603065316, -2.79446235029201,
1.53128139579193, 2.65531461899762, -1.33692502254506, -1.40309746138702,
7.1503230589606, -3.03232696595877, -3.99799016693196, 1.36106448869886,
2.48321747034381, -1.81760324077147, -2.93792828398095, -4.24318581706171,
-5.02935968459404, -2.8067837982494, -3.18182441332852, -3.17383228514224,
-3.75210155230786, -3.87059449392306, -3.38923922246314, -3.8492188722049,
-4.54207754798088, -3.65594980956956, -4.083692310451, 1.30846845334021,
-2.14067263149232, 0.87594188820449, -2.36549598411797, -1.6026649015159,
3.06122694909985, -2.80125769402892, -1.21061533381979, 2.45621751938726,
-1.29479960958983, -1.32215076834374, -3.40665434275434, -3.87581533996246,
2.25817917663687, -4.47531796954113, -3.0627832955887), `24v0.logFC` = c(-5.63038870666508,
-6.44381761773614, -9.71855520117086, 2.83617729007977, -3.39220401105853,
-9.03206578896647, 4.72316302034493, 2.69600384046815, -3.70005627557764,
1.57004091545159, -2.69837465613229, -1.96398920017231, -6.78881951127235,
2.88395305274553, 2.2815237548239, 4.26299665984366, 1.80759581466942,
-6.13159865127454, 2.31984691801986, 2.42017005000243, 1.90522291464817,
-6.08799743507621, 1.93712248733409, -2.85900752558212, -2.53314737897312,
-5.39825301043452, -2.67723981805717, 1.99371094854688, -6.30018858053764,
-5.51269257077614, 1.22241656195939, 2.21258569930309, -3.34360626488873,
2.84622416535537, 1.62066860166637, -4.78372642162402, -5.70058663575352,
-5.38484766744183, -5.09622306588101, 1.4720608195385, 2.81999896316429,
-5.00082275569591, -4.22688753149231, 3.13575564914666, -3.65108053066376,
-4.96357882349085, 1.41296031244116, 2.02777799934303, -2.69828900301196,
3.3280993785528, -4.91524305899842, 1.33559456145784, -3.20704073441948,
3.4589509748052, -4.76220190139647, -5.41080867153898, -2.52368474702804,
2.35227625126478, -1.66843249760588, -5.98409816921998, 1.97783827997787,
-5.05713668502354, -5.66050145489872, -3.23578154397873, -5.47590590732768,
1.79478248294125, 1.53715323981994, 1.17924158724962, -3.38135082340585,
-6.0533528170201, 3.64133939415158, -5.18911327526881, 1.12836935335145,
-5.50962223997537, -5.13018380364927, 2.43347308480167, 3.09821960202745,
1.08371028206663, 1.15073045150599, 1.66618928654645, -4.57926216582183,
-4.39710926752864, 2.29767598663356, 1.11369483753166, 1.90655785334723,
-4.57402537033236, -1.64795888694238, -4.90903507231585, -4.69225607476014,
-3.78658812613927, 2.71158710397443, 2.0014280958994, 1.42976609883019,
1.69469179645011, 2.02542067939079, 1.29620189704961, -2.11262483065018,
-4.45917113273047, -5.00209626672923, -2.22313144605634, 2.28015100958449,
2.99018492359164, -4.55487718731015, -4.4548605761304, -2.74564764773372,
1.50453692418261, -1.80916071906589, -3.78464149018, 5.44303553967597,
1.63876298449708, -2.55262897048338, -4.52516489294599, 1.86420921266273,
1.92394285952202, -1.87942876800687, -4.72261410600152, -2.57805004195914,
-4.47016439150427, -2.79823793519691, 1.49664827810044, -5.12377787039451,
-1.64302437544027, -2.32977998370269, -1.54481646746192, 1.51575792656483,
-1.66914425228449, -5.67499677983564, -2.41648326328024, 1.4178331810756,
1.34477707145884, 1.6168542644666, 1.51101526992765, 2.85584943327062,
-4.19149110174371, 1.69988651757158, 2.42852927065545, 1.66468253437414,
3.21304236161755, 2.52923647014318, -3.32575501787469, 1.32217968706286,
-4.3404208961354, 1.60065358035513, 1.11559852847827, 2.24622871386089,
-2.07936338453629, -4.096641542765, -1.81072034410725, -1.71317603678557,
3.48120770906848, -3.87801910383696, -4.08020770103266, 1.11478378387699,
-3.46079055745071, 1.7322630648363, 1.54861321913459, 1.96850013686628,
2.50282552968501, -1.81077760463368, 2.02487013804019, 1.95047564703855,
-2.99930567848563, -1.79193606152703, -2.54724389366544, -4.15112678166129,
1.13069494609896, -1.43241741260756, 2.36383913233529, 1.65786717131242,
-2.29507370488367, 1.5786351022915, -3.31042999789919, 2.582031624397,
-3.62438687496853, -4.20056844288066, 1.59787844454997, -4.39077226372811,
-4.50758458577604, -3.5310616440773, 1.34960913912344, 1.56097908933761,
1.80910781056029, -4.1248594098749, 1.5318048796277, -3.97753680837585,
-2.80501778688795, -3.20129716412254, -1.98305552020422, -4.25045728757901,
-4.08742040090535, -4.31190449886363, 1.39895836419766, 1.67102855271528,
-5.13520809951127, -4.31701572736756, 0.832778953035371, 2.30899115735356,
-3.18620970041199, -3.80052092618969, 1.59131320240714)), row.names = c("ENSMUSG00000097417.1",
"ENSMUSG00000084950.7", "ENSMUSG00000021824.12", "ENSMUSG00000085505.1",
"ENSMUSG00000018381.15", "ENSMUSG00000098716.1", "ENSMUSG00000092134.1",
"ENSMUSG00000097272.1", "ENSMUSG00000039206.13", "ENSMUSG00000097860.1",
"ENSMUSG00000047880.14", "ENSMUSG00000087700.1", "ENSMUSG00000094127.1",
"ENSMUSG00000101162.2", "ENSMUSG00000079179.9", "ENSMUSG00000107317.1",
"ENSMUSG00000102451.1", "ENSMUSG00000098609.1", "ENSMUSG00000109669.1",
"ENSMUSG00000086836.1", "ENSMUSG00000092545.7", "ENSMUSG00000106093.1",
"ENSMUSG00000090625.1", "ENSMUSG00000009585.17", "ENSMUSG00000047988.1",
"ENSMUSG00000027950.7", "ENSMUSG00000086430.1", "ENSMUSG00000105941.1",
"ENSMUSG00000004071.12", "ENSMUSG00000104208.1", "ENSMUSG00000022223.8",
"ENSMUSG00000031145.15", "ENSMUSG00000047369.15", "ENSMUSG00000097347.1",
"ENSMUSG00000086680.2", "ENSMUSG00000097660.1", "ENSMUSG00000002109.14",
"ENSMUSG00000108524.1", "ENSMUSG00000067786.16", "ENSMUSG00000021485.13",
"ENSMUSG00000097768.1", "ENSMUSG00000085437.8", "ENSMUSG00000086220.1",
"ENSMUSG00000075576.12", "ENSMUSG00000102299.1", "ENSMUSG00000079293.11",
"ENSMUSG00000090952.2", "ENSMUSG00000040723.14", "ENSMUSG00000027329.9",
"ENSMUSG00000106005.1", "ENSMUSG00000085425.1", "ENSMUSG00000046027.17",
"ENSMUSG00000006517.5", "ENSMUSG00000004221.16", "ENSMUSG00000049775.16",
"ENSMUSG00000109881.1", "ENSMUSG00000023988.8", "ENSMUSG00000092371.1",
"ENSMUSG00000087611.1", "ENSMUSG00000102349.1", "ENSMUSG00000092103.2",
"ENSMUSG00000089636.1", "ENSMUSG00000086360.1", "ENSMUSG00000036138.16",
"ENSMUSG00000045382.6", "ENSMUSG00000090081.1", "ENSMUSG00000107999.1",
"ENSMUSG00000041774.14", "ENSMUSG00000030720.16", "ENSMUSG00000027398.13",
"ENSMUSG00000090246.1", "ENSMUSG00000104823.1", "ENSMUSG00000032528.4",
"ENSMUSG00000100783.1", "ENSMUSG00000082848.1", "ENSMUSG00000021250.13",
"ENSMUSG00000009292.17", "ENSMUSG00000109491.1", "ENSMUSG00000090252.1",
"ENSMUSG00000087443.1", "ENSMUSG00000110148.1", "ENSMUSG00000085754.7",
"ENSMUSG00000064339.1", "ENSMUSG00000091160.1", "ENSMUSG00000108332.1",
"ENSMUSG00000108321.1", "ENSMUSG00000048721.2", "ENSMUSG00000032363.15",
"ENSMUSG00000003868.14", "ENSMUSG00000107075.1", "ENSMUSG00000102961.1",
"ENSMUSG00000096957.1", "ENSMUSG00000093006.2", "ENSMUSG00000109754.1",
"ENSMUSG00000027774.16", "ENSMUSG00000090016.1", "ENSMUSG00000012126.16",
"ENSMUSG00000074136.9", "ENSMUSG00000089875.1", "ENSMUSG00000092190.1",
"ENSMUSG00000074213.2", "ENSMUSG00000086564.9", "ENSMUSG00000055725.11",
"ENSMUSG00000028451.12", "ENSMUSG00000078786.9", "ENSMUSG00000092655.1",
"ENSMUSG00000098482.1", "ENSMUSG00000087442.1", "ENSMUSG00000108574.1",
"ENSMUSG00000087669.1", "ENSMUSG00000053442.4", "ENSMUSG00000101674.6",
"ENSMUSG00000078653.4", "ENSMUSG00000108353.1", "ENSMUSG00000097748.1",
"ENSMUSG00000050010.8", "ENSMUSG00000029623.14", "ENSMUSG00000097060.1",
"ENSMUSG00000050627.13", "ENSMUSG00000092386.1", "ENSMUSG00000085471.1",
"ENSMUSG00000087607.1", "ENSMUSG00000061607.14", "ENSMUSG00000026729.9",
"ENSMUSG00000097854.1", "ENSMUSG00000064330.9", "ENSMUSG00000039753.16",
"ENSMUSG00000027706.8", "ENSMUSG00000029360.3", "ENSMUSG00000103864.1",
"ENSMUSG00000085208.3", "ENSMUSG00000073628.5", "ENSMUSG00000102776.1",
"ENSMUSG00000032487.8", "ENSMUSG00000086998.1", "ENSMUSG00000006728.12",
"ENSMUSG00000021273.9", "ENSMUSG00000062591.5", "ENSMUSG00000040936.14",
"ENSMUSG00000002985.16", "ENSMUSG00000046145.7", "ENSMUSG00000038274.11",
"ENSMUSG00000087366.7", "ENSMUSG00000039452.5", "ENSMUSG00000087660.1",
"ENSMUSG00000087206.1", "ENSMUSG00000079355.3", "ENSMUSG00000086693.1",
"ENSMUSG00000030579.10", "ENSMUSG00000084878.1", "ENSMUSG00000108235.1",
"ENSMUSG00000097077.1", "ENSMUSG00000101438.1", "ENSMUSG00000087367.1",
"ENSMUSG00000064351.1", "ENSMUSG00000086312.1", "ENSMUSG00000097221.1",
"ENSMUSG00000107090.1", "ENSMUSG00000037982.16", "ENSMUSG00000095280.1",
"ENSMUSG00000017837.12", "ENSMUSG00000029422.14", "ENSMUSG00000102564.1",
"ENSMUSG00000026749.11", "ENSMUSG00000067288.12", "ENSMUSG00000025353.4",
"ENSMUSG00000109251.1", "ENSMUSG00000086600.8", "ENSMUSG00000052214.9",
"ENSMUSG00000090330.2", "ENSMUSG00000090009.1", "ENSMUSG00000059791.14",
"ENSMUSG00000101695.1", "ENSMUSG00000103698.1", "ENSMUSG00000097472.1",
"ENSMUSG00000084128.10", "ENSMUSG00000026156.8", "ENSMUSG00000075297.10",
"ENSMUSG00000055692.18", "ENSMUSG00000038502.16", "ENSMUSG00000015942.9",
"ENSMUSG00000106478.1", "ENSMUSG00000039307.16", "ENSMUSG00000054556.6",
"ENSMUSG00000045576.16", "ENSMUSG00000037936.15", "ENSMUSG00000087486.1",
"ENSMUSG00000082180.1", "ENSMUSG00000014470.4", "ENSMUSG00000056078.5",
"ENSMUSG00000049928.15", "ENSMUSG00000013997.10", "ENSMUSG00000105198.1",
"ENSMUSG00000086714.1", "ENSMUSG00000109274.1", "ENSMUSG00000100548.1",
"ENSMUSG00000053706.2", "ENSMUSG00000109854.1", "ENSMUSG00000060802.8",
"ENSMUSG00000093598.1"), class = "data.frame")

```

Please note, I have checked out many posts regarding this, most notably: https://www.biostars.org/p/287512/#287518, which highlights how to perform this procedure. Regardless, there are some unanswered questions.

The problem is the following: After running the code suggested in the post, I do indeed obtain x number of clusters, depending on the height I decide to cut the tree at. However, they are in an incomprehensible order, making it extremely difficult to extract the clusters as they are observed in the heatmap. Please note, the dataframe I have uploaded is only a snippet of the true data - and is called logFC\_mat, in the code below:

```
pheatmap(logFC_mat,
         show_rownames=F, cluster_cols=F, cluster_rows=T, scale="row",
         cex=1, clustering_distance_rows="euclidean", cex=1,
         clustering_distance_cols="euclidean",
         clustering_method="complete", border_color=FALSE) -> heatof

#make clusters based on h=1.5
sort(cutree(heatof$tree_row, h=1.5)) -> hi_clusters

#check if clusters are in same order as in heatmap
hi_clusters[hi_clusters==1] %>% names() -> clust1
hi_clusters[hi_clusters==2] %>% names() -> clust2

length(clust1) #56
length(clust2) #78
all(clust1 %in%
      rownames(logFC_mat[heatof$tree_row$order[1:56],])) #TRUE

all(clust2 %in%
      rownames(logFC_mat[heatof$tree_row$order[57:134],])) #FALSE

```

Now, from visual inspection of the heatmap, it becomes clear that cluster 2 is in fact the last 78 genes rather than the middle cluster. The problem is, why is it labeled as cluster 2 then !? In addition, this is a snippet of the data. The real data contains in fact thousands of rows, and so such visual deductions can simply be impossible.

Now for the question: How to extract the clusters, such that they are in the SAME ORDER, as what is seen in the heatmap (Top to bottom)?

Thanks.

# Answer
> 1 votes
I would recommend not to use `h` to cut the tree but to use `k`, as `k` instructs the number of clusters to obtain, and then leaves it to the function to find the correct height for that. From there on you can easily extract the genes per cluster, for example with this code obtaining the first three clusters:

heatof <-
pheatmap(logFC_mat,
         show_rownames=F, cluster_cols=F, cluster_rows=T, scale="row",
         cex=1, clustering_distance_rows="euclidean", cex=1,
         clustering_distance_cols="euclidean",
         clustering_method="complete", border_color=FALSE)

#/ extract the top three clusters
k <- 3
sort(cutree(heatof$tree_row, k=3)) -> hi_clusters

genes_per_cluster <- do.call(rbind, lapply(unique(hi_clusters), function(x){

  data.frame(cluster=x, gene=names(hi_clusters[hi_clusters %in% x]))

}))

genes_per_cluster
#>     cluster                  gene
#> 1         1  ENSMUSG00000097417.1
#> 2         1  ENSMUSG00000084950.7
#> 3         1 ENSMUSG00000021824.12
#> 4         1 ENSMUSG00000018381.15
#> 5         1  ENSMUSG00000098716.1
#> 6         1 ENSMUSG00000039206.13
#> 7         1 ENSMUSG00000047880.14
#> 8         1  ENSMUSG00000094127.1
#> 9         1  ENSMUSG00000098609.1
#> 10        1  ENSMUSG00000106093.1
#> 11        1 ENSMUSG00000009585.17
#> 12        1  ENSMUSG00000027950.7
#> 13        1 ENSMUSG00000004071.12
#> 14        1  ENSMUSG00000104208.1
#> 15        1 ENSMUSG00000047369.15
#> 16        1  ENSMUSG00000097660.1
#> 17        1 ENSMUSG00000002109.14
#> 18        1  ENSMUSG00000108524.1
#> 19        1 ENSMUSG00000067786.16
#> 20        1  ENSMUSG00000085437.8
#> 21        1  ENSMUSG00000086220.1
#> 22        1  ENSMUSG00000102299.1
#> 23        1 ENSMUSG00000079293.11
#> 24        1  ENSMUSG00000085425.1
#> 25        1  ENSMUSG00000006517.5
#> 26        1 ENSMUSG00000049775.16
#> 27        1  ENSMUSG00000109881.1
#> 28        1  ENSMUSG00000102349.1
#> 29        1  ENSMUSG00000089636.1
#> 30        1  ENSMUSG00000086360.1
#> 31        1 ENSMUSG00000036138.16
#> 32        1  ENSMUSG00000045382.6
#> 33        1 ENSMUSG00000030720.16
#> 34        1 ENSMUSG00000027398.13
#> 35        1  ENSMUSG00000104823.1
#> 36        1  ENSMUSG00000100783.1
#> 37        1  ENSMUSG00000082848.1
#> 38        1  ENSMUSG00000110148.1
#> 39        1  ENSMUSG00000085754.7
#> 40        1  ENSMUSG00000108321.1
#> 41        1 ENSMUSG00000032363.15
#> 42        1 ENSMUSG00000003868.14
#> 43        1  ENSMUSG00000107075.1
#> 44        1  ENSMUSG00000074136.9
#> 45        1  ENSMUSG00000089875.1
#> 46        1 ENSMUSG00000055725.11
#> 47        1 ENSMUSG00000028451.12
#> 48        1  ENSMUSG00000078786.9
#> 49        1  ENSMUSG00000087442.1
#> 50        1  ENSMUSG00000053442.4
#> 51        1  ENSMUSG00000101674.6
#> 52        1  ENSMUSG00000097748.1
#> 53        1  ENSMUSG00000050010.8
#> 54        1 ENSMUSG00000029623.14
#> 55        1  ENSMUSG00000097060.1
#> 56        1 ENSMUSG00000050627.13
#> 57        1  ENSMUSG00000085471.1
#> 58        1 ENSMUSG00000039753.16
#> 59        1  ENSMUSG00000027706.8
#> 60        1  ENSMUSG00000032487.8
#> 61        1 ENSMUSG00000002985.16
#> 62        1 ENSMUSG00000038274.11
#> 63        1  ENSMUSG00000087206.1
#> 64        1  ENSMUSG00000079355.3
#> 65        1  ENSMUSG00000086693.1
#> 66        1 ENSMUSG00000030579.10
#> 67        1  ENSMUSG00000108235.1
#> 68        1  ENSMUSG00000097077.1
#> 69        1  ENSMUSG00000087367.1
#> 70        1 ENSMUSG00000037982.16
#> 71        1 ENSMUSG00000029422.14
#> 72        1 ENSMUSG00000026749.11
#> 73        1 ENSMUSG00000067288.12
#> 74        1  ENSMUSG00000090330.2
#> 75        1 ENSMUSG00000059791.14
#> 76        1  ENSMUSG00000103698.1
#> 77        1  ENSMUSG00000097472.1
#> 78        1  ENSMUSG00000026156.8
#> 79        1 ENSMUSG00000075297.10
#> 80        1 ENSMUSG00000055692.18
#> 81        1 ENSMUSG00000039307.16
#> 82        1 ENSMUSG00000045576.16
#> 83        1  ENSMUSG00000087486.1
#> 84        1  ENSMUSG00000014470.4
#> 85        1  ENSMUSG00000056078.5
#> 86        1 ENSMUSG00000049928.15
#> 87        1  ENSMUSG00000086714.1
#> 88        1  ENSMUSG00000109274.1
#> 89        1  ENSMUSG00000109854.1
#> 90        1  ENSMUSG00000060802.8
#> 91        2  ENSMUSG00000085505.1
#> 92        2  ENSMUSG00000092134.1
#> 93        2  ENSMUSG00000097272.1
#> 94        2  ENSMUSG00000087700.1
#> 95        2  ENSMUSG00000101162.2
#> 96        2  ENSMUSG00000107317.1
#> 97        2  ENSMUSG00000102451.1
#> 98        2  ENSMUSG00000109669.1
#> 99        2  ENSMUSG00000086836.1
#> 100       2  ENSMUSG00000092545.7
#> 101       2  ENSMUSG00000090625.1
#> 102       2  ENSMUSG00000047988.1
#> 103       2  ENSMUSG00000086430.1
#> 104       2  ENSMUSG00000105941.1
#> 105       2  ENSMUSG00000022223.8
#> 106       2 ENSMUSG00000031145.15
#> 107       2  ENSMUSG00000086680.2
#> 108       2 ENSMUSG00000021485.13
#> 109       2  ENSMUSG00000097768.1
#> 110       2 ENSMUSG00000075576.12
#> 111       2  ENSMUSG00000090952.2
#> 112       2 ENSMUSG00000040723.14
#> 113       2  ENSMUSG00000027329.9
#> 114       2  ENSMUSG00000106005.1
#> 115       2 ENSMUSG00000046027.17
#> 116       2 ENSMUSG00000004221.16
#> 117       2  ENSMUSG00000023988.8
#> 118       2  ENSMUSG00000087611.1
#> 119       2  ENSMUSG00000092103.2
#> 120       2  ENSMUSG00000090081.1
#> 121       2  ENSMUSG00000107999.1
#> 122       2 ENSMUSG00000041774.14
#> 123       2  ENSMUSG00000090246.1
#> 124       2  ENSMUSG00000032528.4
#> 125       2 ENSMUSG00000021250.13
#> 126       2 ENSMUSG00000009292.17
#> 127       2  ENSMUSG00000109491.1
#> 128       2  ENSMUSG00000090252.1
#> 129       2  ENSMUSG00000087443.1
#> 130       2  ENSMUSG00000064339.1
#> 131       2  ENSMUSG00000091160.1
#> 132       2  ENSMUSG00000108332.1
#> 133       2  ENSMUSG00000048721.2
#> 134       2  ENSMUSG00000102961.1
#> 135       2  ENSMUSG00000096957.1
#> 136       2  ENSMUSG00000093006.2
#> 137       2 ENSMUSG00000027774.16
#> 138       2 ENSMUSG00000012126.16
#> 139       2  ENSMUSG00000092190.1
#> 140       2  ENSMUSG00000074213.2
#> 141       2  ENSMUSG00000086564.9
#> 142       2  ENSMUSG00000092655.1
#> 143       2  ENSMUSG00000098482.1
#> 144       2  ENSMUSG00000108574.1
#> 145       2  ENSMUSG00000087669.1
#> 146       2  ENSMUSG00000078653.4
#> 147       2  ENSMUSG00000108353.1
#> 148       2  ENSMUSG00000092386.1
#> 149       2  ENSMUSG00000087607.1
#> 150       2 ENSMUSG00000061607.14
#> 151       2  ENSMUSG00000026729.9
#> 152       2  ENSMUSG00000064330.9
#> 153       2  ENSMUSG00000103864.1
#> 154       2  ENSMUSG00000085208.3
#> 155       2  ENSMUSG00000073628.5
#> 156       2  ENSMUSG00000102776.1
#> 157       2  ENSMUSG00000086998.1
#> 158       2 ENSMUSG00000006728.12
#> 159       2  ENSMUSG00000021273.9
#> 160       2  ENSMUSG00000062591.5
#> 161       2 ENSMUSG00000040936.14
#> 162       2  ENSMUSG00000046145.7
#> 163       2  ENSMUSG00000087366.7
#> 164       2  ENSMUSG00000039452.5
#> 165       2  ENSMUSG00000087660.1
#> 166       2  ENSMUSG00000084878.1
#> 167       2  ENSMUSG00000101438.1
#> 168       2  ENSMUSG00000064351.1
#> 169       2  ENSMUSG00000086312.1
#> 170       2  ENSMUSG00000097221.1
#> 171       2  ENSMUSG00000107090.1
#> 172       2  ENSMUSG00000095280.1
#> 173       2 ENSMUSG00000017837.12
#> 174       2  ENSMUSG00000109251.1
#> 175       2  ENSMUSG00000086600.8
#> 176       2  ENSMUSG00000090009.1
#> 177       2  ENSMUSG00000101695.1
#> 178       2 ENSMUSG00000084128.10
#> 179       2 ENSMUSG00000038502.16
#> 180       2  ENSMUSG00000015942.9
#> 181       2  ENSMUSG00000106478.1
#> 182       2  ENSMUSG00000054556.6
#> 183       2  ENSMUSG00000082180.1
#> 184       2 ENSMUSG00000013997.10
#> 185       2  ENSMUSG00000100548.1
#> 186       2  ENSMUSG00000053706.2
#> 187       2  ENSMUSG00000093598.1
#> 188       3  ENSMUSG00000097860.1
#> 189       3  ENSMUSG00000079179.9
#> 190       3  ENSMUSG00000097347.1
#> 191       3  ENSMUSG00000092371.1
#> 192       3  ENSMUSG00000109754.1
#> 193       3  ENSMUSG00000090016.1
#> 194       3  ENSMUSG00000097854.1
#> 195       3  ENSMUSG00000029360.3
#> 196       3  ENSMUSG00000102564.1
#> 197       3  ENSMUSG00000025353.4
#> 198       3  ENSMUSG00000052214.9
#> 199       3 ENSMUSG00000037936.15
#> 200       3  ENSMUSG00000105198.1

Created on 2022-07-07 by the reprex package (v2.0.1)

```

---
Tags: r, rna-seq, clustering, unsupervised-learning
---


# Question
Title: How to retrieve dataset names from python statsmodels?
The Python statsmodels library grants access to many R-datasets which can be downloaded e.g. with:

```
import statsmodels.api as sm
df = sm.datasets.get_rdataset(dataname='iris', package='datasets').data

```

But for this you have to know the id or the name of the dataset.

```
sm.datasets.get_rdataset??
Signature: sm.datasets.get_rdataset(dataname, package='datasets', cache=False)
Source:
def get_rdataset(dataname, package="datasets", cache=False):
    """download and return R dataset

    Parameters
    ----------
    dataname : str
        The name of the dataset you want to download

```

There is a list of available datasets here, but this list does not contain the identifiers. Also not the detailed view.

Is there an overview of all the available datasets that contains the ids e.g. "iris"?

There is an index that is supposed to list all the datasets, but e.g. 'rice' as a name does not work.

# Answer
> 1 votes
Here is the website of the creator of rdatasets I guess with the info you are looking for:

http://vincentarelbundock.github.io/Rdatasets/datasets.html

Here is a short look at it:

```
Package     Item    Title   Rows    Cols    n_binary    n_character     n_factor    n_logical   n_numeric   CSV     Doc
boot    acme    Monthly Excess Returns  60  3   0   1   0   0   2   CSV     DOC
boot    aids    Delay in AIDS Reporting in England and Wales    570     6   1   0   0   0   6   CSV     DOC
boot    aircondit   Failures of Air-conditioning Equipment  12  1   0   0   0   0   1   CSV     DOC
boot    aircondit7  Failures of Air-conditioning Equipment  24  1   0   0   0   0   1   CSV     DOC
boot    amis    Car Speeding and Warning Signs  8437    4   1   0   0   0   4   CSV     DOC
boot    aml     Remission Times for Acute Myelogenous Leukaemia     23  3   2   0   0   0   3   CSV     DOC
...

```

I assume the Item column are the identifiers.

# Answer
> 1 votes
You can view all the available datasets at this URL https://github.com/vincentarelbundock/Rdatasets/blob/master/datasets.csv

The available ones in statsmodels are those with the tag `datasets` which looks like line 372 - 455.

I didn't test them all, but everyone I tried from a small random sample seemed to work fine.

# Answer
> 1 votes
Please, try this code to get all modules in datasets:

```
import statsmodels.api as sm
import inspect

def get_list_of_modules_in_datasets():
  # get list of strings for the objects in datasets
  list_str_in_datasets = sm.datasets.__all__
  print(list_str_in_datasets)
  # get list of objects in datasets
  list_objs_in_datasets = [eval('sm.datasets.' + dsstr) for dsstr in list_str_in_datasets]
  print(list_objs_in_datasets)
  # get list of modules in datasets
  list_modules_in_datasets = []
  for dso in list_objs_in_datasets:
    if (inspect.ismodule(dso)):
      list_modules_in_datasets.append(dso)
  print(list_modules_in_datasets)
  print("Number of Modules in datasets = %d" %(len(list_modules_in_datasets)))
  return (list_modules_in_datasets)

if __name__ == '__main__':

  list_modules_in_datasets = get_list_of_modules_in_datasets()

  # print data for all dataframes from the datasets
  for ds in list_modules_in_datasets:
    df = ds.load_pandas().data
    print(list(df.columns))
    print(df.shape)
    print(df)

```

---
Tags: python, statsmodels
---


# Question
Title: Bootstrap analysis - why it is called bootstrap?
*Bootstrap analysis*, *bootstrapping* etc are quite common jargons of bioinformatics and phylogenetics.

However, it is not very much clear to us, what exactly being meant by *"a boot's straps".*

Does it means a *"comparison"*? (Such as we hold the two straps of a boot (shoe) in "two hands"; may be an analogy to comparing 2 things (say left-'hand' side and right-'hand' side of an algebraic proof-task)) or it is something about *sequence alignment* (since we 'tie' a knot in boot (shoe) so that the 2 ends of the strap remain attached to each-other).

Some of my classmates even assumed it may be a process during computer booting (but I could not agree with that).

I've searched internet, but could not found anything very helpful. This wikipedia page (permalink) about bootstrapping, tells,

>  Tall boots may have a tab, loop or handle at the top known as a bootstrap, allowing one to use fingers or a boot hook tool to help pulling the boots on. The saying "to pull oneself up by one's bootstraps" was already in use during the 19th century as an example of an impossible task. The idiom dates at least to 1834, when it appeared in the Workingman's Advocate: "It is conjectured that Mr. Murphee will now be enabled to hand himself over the Cumberland river or a barn yard fence by the straps of his boots." In 1860 it appeared in a comment on philosophy of mind: ***"The attempt of the mind to analyze itself an effort analogous to one who would lift himself by his own bootstraps"***. **Bootstrap as a metaphor, meaning to better oneself by one's own unaided efforts**, was in use in 1922. This metaphor spawned additional metaphors for a series of self-sustaining processes that proceed without external help.
>
> The term is sometimes attributed to a story in Rudolf Erich Raspe's The Surprising Adventures of Baron Munchausen, but in that story Baron Munchausen pulls himself (and his horse) out of a swamp by his hair (specifically, his pigtail), not by his bootstraps – and no explicit reference to bootstraps has been found elsewhere in the various versions of the Munchausen tales.

(In the quotation, highlighted portion seemed to tried to tell a meaning), but it wasn't also very much helpful. Could not correlate with bioinformatics. (Btw the same wikipedia article mentions computer booting).

Wikipedia article about Bootstrapping (statistics) (permalink) tells:

> The bootstrap was published by Bradley Efron in "Bootstrap methods: another look at the jackknife" (1979), inspired by earlier work on the jackknife.Improved estimates of the variance were developed later. A Bayesian extension was developed in 1981. The bias-corrected and accelerated (BCa) bootstrap was developed by Efron in 1987, and the ABC procedure in 1992.

Still that is not very helpful.

However, the textbook, *Introduction to Bioinformatics/ Arthur M. Lesk/ 3rd Edition*/ Oxford / (Low Price edition) ; in its *chapter 5* (Alignments and phylogenetic trees), section The problem of varying rates of evolution -\> *computational consideration* (page 296); a clearer definition has been given:

> 3. Formal statistical tests, involving rerunning the calculation on subsets of the original data, are known as **jackknifing** and **bootstrapping**.
>     * **jackknifing** is calculation with data sets sampled randomly from the original data. For phylogeny calculations from from multiple sequence alignments, select different subsets of the position in the alignment, and rerun the calculation. Finding that each subset gives the same phylogenetic tree lends it credibility. If each subset gives a different tree, none of them are trustworthy.
>     * **Bootstrapping** is similar to jackknifing except that the position chosen at random may include multiple copies of the same position, to form data sets of the same size as original, to preserve statistical properties of data sampling.

More clear definition of "Boot-strapping"; but it does not explain, where the relation exists with the boot (shoe) and its straps.

# Answer
The method itself has nothing really to do with boots or straps, just as the jack-knife method also has nothing to do with knives. In the case of bootstrapping, the goal is to determine the accuracy of an estimate from random subsets of a population. Normally estimating something like the variance of the mean requires multiple independent samples, but bootstrapping allows you to perform estimates from a single population. This is essentially allowing the estimate to "pull itself up by its bootstraps", which is the only reasonable source for the term. See also the more general wikipedia article on bootstrapping outside of statistics.
> 9 votes

# Answer
In general, branch support values indicate how much of the total data support that branch in a proposed tree. But what do you do if you used all of the available data to infer the tree in the first place?

With Bootstrap supports, you *generate* new alignments (almost) *ex nihilo*. You calculate a trees support values, using its own data.

This is the connection to "pulling yourself up by your bootstraps": you don't "pull" on additionaly collected data (a rope or something in the metaphor), but rather on the very thing that your observation is based on (your boots).
> 6 votes

# Answer
**Origin of name** Bootstrapping is a term of derision in many ways, prior widespread use in phylogenetics. The response here is also the story that I recall the name meaning 'pulling yourself up by your bootstraps'. If I remember the story of it origin correctly, it means its impossible to 'pull-yourself-up' with you bootstraps because there's nothing external to leverage.

**Historic usage** Thus traditionally, it was a last resort approach in statistics, or else a crude data mining approach, because a lot statistics at that time worked by comparing a result against a theoretical reference, i.e. critical probability. So for example, T-test, Chi-squared, ANOVA etc ... a theoretical null distribution is used or else created against which the observation is compared. The idea of not using any theoretical external distribution was not deemed cool in the days of critical probability. Bootstrapping does weird things to the distribution that few really understand (spent time investigating).

**Method** Technically bootstrapping is simply resampling with replacement. So if you've a data set of 100 nucleotides replicates of a resampled data set are constructed each of 100 nucleotides by randomly selecting 100 nucleotides (with replacement) from the original data set. Performing tree estimate calculations over multiple replicates permits a robustness estimate of each 'node' in the tree (hypothetical common ancestor).

**History in phylogenetics** With onset of phylogenetics however bootstrapping became a statistical norm for several different reasons, which are quite complex. For example, attempting to perform critical probability across ever branch of a tree could run into statistical problems, what is a 'null distribution' in evolutionary biology (that question has since been reasonably well answered), in the days of the rows about the parsimony method - bootstrapping was well under radar. Hence bootstrapping is now mainstream, with some caveats.

I assume bootstrapping was originally supported by simulation studies, however it so established I don't even know it theoretical history.

**Jack-knifing** 'Jack-knifing' is named because it was deemed a rough and ready tool (a 'pen-knife') to assess a data set and a forerunner of the bootstrap. Its so old now I don't even remember the resampling strategy it uses, chopping the data set in half rings a bell - I honestly forget.

**Today** These days there is no real reason to know the method, other than using it and the thresholds commonly accepted. Very little statistical development is based around bootstrapping. The underlying tenant is that each site is independent of the others and that is a widespread assumption across phylogenetics. That assumptions has been assessed (obviously) but the solutions were not implemented. Continually reassessing this assumption remains very interesting and could be productive longterm, but it doesn't negate the basic approach to treat each alignment position as an independent event to every other alignment position because its a baseline used throughout molecular phylogenetics.
> 0 votes

---
Tags: terminology
---


# Question
Title: AttributeError [in fastStructure's chooseK.py]: module 'vars' has no attribute 'insum'
I've been working with the fastStructure program and am on the step of analyzing model complexity using the provided chooseK.py script. I have been running into the following error:

```
Traceback (most recent call last):
  File "chooseK.py", line 109, in <module>
    bestKs = parse_varQs(files)
  File "chooseK.py", line 54, in parse_varQs
    Q = Q/utils.insum(Q,[1])
AttributeError: module 'vars' has no attribute 'insum'

```

I have run into 'AttributeError' messages in the past and have attempted the solution of renaming files in my working directory that may be shadowing the real module (similar to: https://stackoverflow.com/questions/36530726/using-numpy-module-object-has-no-attribute-array), but this hasn't fixed anything. Reinstalling Numpy/Scipy/fastStructure has also been a dead end.

Any suggestions are much appreciated! Happy to provide more info as needed.

-A

# Answer
Could you try loading `fastStructure` via conda?

```
conda create --name faststruc python=3.9 # there can be issues with 3.10
conda activate faststruc
conda install -c bioconda faststructure

```

It will then ask you whether you want to install `numpy` and `scipy` automatically. Furthermore all dependencies are installed and connectivity will have already been tested. You can easily delete everything afterwards by deleting the environment `faststruc` so keeps everything clean.

---

The code in question is here:

https://programtalk.com/python-examples/vars.utils.insum/

Its a `numpy` calculation, what they are doing is navigating Python by leveraging the power of `numpy`, which is very fast.

Code ...

```
bestKs = []
for file in files:
    handle = open(file,'r')
    Q = np.array([map(float,line.strip().split()) for line in handle])
    Q = Q/utils.insum(Q,[1])
    handle.close()
    N = Q.shape[0]
    C = np.cuemsum(np.sort(Q.sum(0))[::-1])
    bestKs.append(np.sum(C<N-1)+1)

return bestKs

```

The code looks good, they know their stuff.
> 0 votes

---
Tags: python, scipy
---


title: Different line length in fasta file

text:
# Question
Title: Different line length in fasta file
I am currently using VEP for variant annotation.

I am facing an error as below:

```
[E::fai_build_core] Different line length in sequence 'Pn9'

```

I understand there is an issue with the difference in line length of Pn9 in the fasta file. However, the sequences of fasta file all have different line length. I don't get why there is an error specifically on Pn9.

I have tried using both snpEff and Annovar but doesn't work. Any thoughts on this would be really appreciated.

I am attaching the length of sequence here :

Thanks in advance.

# Answer
This could be a bug in the fasta file reader component of fai\_build\_core. The fasta file format is very relaxed, and allows some odd constructions, for example:

```
>cb3c6a5f-17dd-4cd6-a9ca-02124d2a32de
AAGAAATTGACCTCCGTAACCCCAAAGGAGACAAAAAGAATGCTGGCTTTACTCCC
CAGGAGCGCCAGGTTTTGGGGAACTGCCTACAAGCTGTACCATTGGCTGACAGCTT
TTAGGGCATCTACCCCCAACACTGCTTACACTTACACTTTCTGGACTTACATG
ATGAATGCCCGCTCTAAGAATGTTGGTTGGCGCCTTGATTACTTTTTGCTTTCC
CACTCTCTTTTTACCTGCATTGTGACAGCAAGATCCGGTCAAGGCTCTTG
GCAGTACACTGTCCCATCACCGCTACCTGTGATACCCCTCCCTGCAGTAGCTT
CCTGCCTGGGAGATGGCTCTCTCTGCAGAAGTCTGGTGTTTTAGCC
TTCAGGTGTTTGGTTTGTATGTGCTCCCTCATTTAAACATTAAACCAAACTT
CTGGTTTCCTTTAGACAATCCAAGAGAAATAAAAGGCCCTAC
TTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

```

A few programs will complain about fasta files that look like this, because it usually means that data has been corrupted. These programs only accept fasta files that have sequences with identical line lengths for all but the last line:

```
>cb3c6a5f-17dd-4cd6-a9ca-02124d2a32de
AAGAAATTGACCTCCGTAACCCCAAAGGAGACAAAAAGAATGCTGGCTTTACTCCC
CAGGAGCGCCAGGTTTTGGGGAACTGCCTACAAGCTGTACCATTGGCTGACAGCTT
TTAGGGCATCTACCCCCAACACTGCTTACACTTACACTTTCTGGACTTACATGATG
AATGCCCGCTCTAAGAATGTTGGTTGGCGCCTTGATTACTTTTTGCTTTCCCACTC
TCTTTTTACCTGCATTGTGACAGCAAGATCCGGTCAAGGCTCTTGGCAGTACACTG
TCCCATCACCGCTACCTGTGATACCCCTCCCTGCAGTAGCTTCCTGCCTGGGAGAT
GGCTCTCTCTGCAGAAGTCTGGTGTTTTAGCCTTCAGGTGTTTGGTTTGTATGTGC
TCCCTCATTTAAACATTAAACCAAACTTCTGGTTTCCTTTAGACAATCCAAGAGAA
ATAAAAGGCCCTACTTTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

```

One way to work around this input error problem is to get rid of *all* line breaks in the sequence of a fasta file, but then you run the risk of programs not properly being able to read all the characters in a line. A better approach is to feed the fasta file through a less restrictive data cleaning program and convert the lines back to all have a consistent length.

There's a good chance this is *not* the problem you are experiencing, but without more information about your input data, all we have to go on are the error messages and the sequence lengths \[not line lengths\] you have provided.
> 1 votes

---
Tags: fasta, sequence-annotation, samtools
---


# Question
Title: Get number of reads with a single, (almost) exact match to the full length of a reference sequence
I can't find an answer to this in previous questions, so hoping someone can help me now.

We have nanopore sequenced a PCR product, and I've filtered our reads to +/- 10 bp of the expected product length (~4,400bp), leaving ~55,000 reads.

We now want to work out how many of the reads correspond to each of 6 different potential arrangements. The section we're interested in is in the middle of the 4,400 bp PCR product, and each different potential arrangement is ~1,500 bp. We have reference sequences for each of the 1,500 bp arrangements. Basically, we want to know what portion of our starting cells had each of the different arrangements of genes in this operon (e.g arrangement 1=30%, arrangement 2=15%, etc.).

This seems like it should be a really easy thing to do using read mapping to count alignments for each reference separately, but I just can't get my head around it. No matter what options I've tried with minimap2 and samtools, every read maps every time on every different arrangement. I think the problem is that the genes in our reference sequences *are* all present in every read, just in different arrangements. So we need to **only** count alignments which align to the whole 1,500 bp reference sequence, not partial hits or hits with huge gaps in them.

I have tried the minimap2 options -r, -G and --secondary=no to restrict the output to single alignments with no large gaps, and combined that with samtools view -m 1000 to output only alignments longer than 1,000 bp (this is far stricter than it needs to be to only output alignments to any of the different potential arrangements). But still, every time, every read is mapping to every reference, usually multiple times.

What am I doing wrong?! Am I totally over-complicating this, or doing something obviously stupid? Almost feel like using blastn might be better at this point...

**TL;DR: how do I use minimap2 and/or samtools (or any other tool better suited to this) to count how many of our ~4,400 bp reads contain *only full-length, almost exact matches* to our 1,500 bp reference sequence, excluding any partial hits?**

My latest attempt (broken into individual steps, though I'd usually pipe the samtools steps together, also not bothered with samtools sort at this point):

```
minimap2 -ax map-ont -r10 -N5 --secondary=no --no-long-join reference.fasta reads.fastq > reads.minimap.sam

samtools view -m 1000 -b -o reads.minimap.bam reads.minimap.sam

samtools -c reads.minimap.bam

```

# Answer
Because minimap2 is a minimiser-based mapper, it's less useful for determining mapping accuracy to single-base precision.

Where precision matters, you can use LAST, bearing in mind that it calculates *local* alignments, and can report multiple hits per query/reference pair. The usual approach I use for LAST mapping is as follows:

1. `lastdb -uNEAR -R01 reference.fa reference.fa`
2. `last-train -Q 1 reference.fa <(cat reads.fastq | head -n 40000) > train.mat`
3. `lastal -p train.mat reference.fa reads.fastq > reads_vs_ref.maf`

That final MAF file can be inspected for mapping. The default output produced from LAST is not very human-readable, so I've written a perl script that converts MAF results into a CSV format with reference-relative statistics that I have found useful:

`maf2csv.pl reads_vs_ref.maf > reads_vs_ref.csv`

This can be inspected in a spreadsheet viewer (e.g. LibreOffice Calc, or Excel) to look at the mapped locations, mapped fraction, identities, etc..

Alternatively, the MAF results can be converted to SAM/BAM format using `maf-convert` (bearing in mind the previously-mentioned caveats about multiple hits):

`maf-convert sam reads_vs_ref.maf | samtools view -b -T reference.fa | samtools sort > reads_vs_ref.bam`
> 1 votes

---
Tags: sam, samtools, nanopore, read-mapping, minimap2
---