# Question
Title: How can I make this Biopython program (to correct erroneous barcodes) run faster, and is there any alternative method?
*This question has also been asked on Biostars*

I am looking forward to getting a valuable suggestion for a bioinformatic problem.

**Background:**

Currently, I am performing a de novo whole genome assembly. At the stage of barcode correction, I lost nearly half of all the reads due to erroneous barcodes.

Barcodes, in this context, are 18-nucleotide-long DNA sequences that were added to the DNA insert during NGS library preparation. Now, I have a file that contains all the 18-base barcodes used in my projects. These barcodes are important because they are unique for reads that came from a DNA molecule and therefore helpful in generating long-range information.

While library construction, 18-base molecular barcodes were used for the linked-read generation, and the barcodes have a high incidence of "N" in positions between 12 and 18, as you may see below.

There is a very large occurrence of Ns in positions 12 to 18 of the barcodes and for this reason, a large pool of erroneous barcodes were filtered.

Now, I have been asked by the company to do one of two things, to recover the lost reads:

1. Either replace the barcodes(barcode fastq file from sequencing core) with the correct ones from a pool of barcodes (This is a text file I received from the company that contains the pool of all 18 base barcode sequences)
2. Since most "N"s are seen in positions 12 to 18 of the 18-base barcode sequences (barcode fastq file from the sequencing core), replace the 11-mers in the barcode file, and replace them with the 11-mers from the correct barcode file (text file).

I have decided to use Python code to execute the first method.

```
#Define the path where fastq file is stored.
import os
from Bio.SeqIO.QualityIO import FastqGeneralIterator
path=os.path.join(os.path.expanduser("~"), "Desktop", "SVA1_S1_R2_001.fastq")

#create a new file handle to write the file
path_new=os.path.join(os.path.expanduser("~"), "Desktop", "trimmedfastq.fastq")
path_ust_barcode=os.path.join(os.path.expanduser("~"), "Desktop", "UST40MBarcodes.txt")
write_file=open(path_new, "w")

#open the text file containg the correct pool of barcodes

with open(path_ust_barcode, "r") as fh:
   barcodes = fh.readlines()

handle=open(path, "r")
output = []

#===================================================
# Return the Hamming distance between string1 and string2.
def hamming_distance(string1, string2):
    # Start with a distance of zero, and count up
    distance = 0
    # Loop over the indices of the string
    L = len(string1)
    for i in range(L):
        # Add 1 to the distance if these two characters are not equal
        if string1[i] != string2[i]:
            distance += 1
    # Return the final count of differences
    return distance
#===================================================

#set correct barcode to "None"
corr_barcode=None
#set the minimum hamming distance to "None"
min_distance=19
#open the fastq file

for title, sequence, qual in FastqGeneralIterator(handle):

    # if a perfect match exists, just output as is
    if sequence in barcodes:
        output.append("@%s\n%s+\n%s\n" % (title, corr_barcode, qual))
        continue

    for read in barcodes:

        h_distance=hamming_distance(sequence, read)

        # since if there was a perfect match we won't be here,
        # if distance is 1, we'll never do better
        if h_distance == 1:
            corr_barcode = read
            break

        if  h_distance <= min_distance:
            min_distance=h_distance
            corr_barcode=read

    output.append("@%s\n%s+\n%s\n" % (title, corr_barcode, qual))

with open(path_new, "a") as wf:
    for line in output:
        wf.write(line)

handle.close()
print("Done")

```

**Problem**

I would like to know if this is the only Pythonic way to tackle this problem. This method is dead slow. My erroneous barcode fastq file is 30+ Gb, and the UST10x barcode file (the text file that contains all the correct barcode pool) is of size 750Mb. I have split this text file into 1000 smaller files, each about 750kb and having ~40,000 lines per file. It's been the third day and the program is still running without outputting anything. I am looking forward to getting valuable suggestions here, with respect to the code, especially how I can make the process faster.

NB: I've tried replacing

```
with open(path_ust_barcode, "r") as fh:
   barcodes = fh.readlines()

```

with

```
barcodes = list()
with open(path_ust_barcode, "r") as fh:
   for b in fh:
      barcodes.append(b.rstrip())

```

Still, the code runs extremely slowly.

I am using an IBM server with a 197GM RAM, and an octa-core processor. In my initial code, it was iterating through all the correct barcodes of the text file (UST10x file), for each erroneous barcode in the fastq file, at the end of every iteration, it opened a text file in append mode and wrote the correct version of fastq reads. Since my server has high memory capacity, I rewrote the code in such a way that it stored all the correct barcodes in the memory, and finally wrote it in a text file, by bringing the write function outside of the loop.

Python is not absolutely necessary. I'm not familiar with C or C++, but that being said, I would really love to get this process done. I'm using an IBM server with Ubuntu (Jammy), and 197GB Ram, and the current local storage has only 150GM free space. if you can provide me suggestions as to how to execute it in C or C++, I can do it.

# Answer
Try profiling (benchmarking) your code to work out where it is slow.

My other suggestions are:

* replace your simple Pure Python hamming\_distance function with something optimised like RapidFuzz - see https://github.com/maxbachmann/RapidFuzz
* cache the corrections you are making - you have tonnes of RAM so a building a dictionary mapping the fragments with N to the corrected fragment ought to be possible. And then for each new sequence, check if you already solved it and if so reuse the same choice.
> 2 votes

---
Tags: sequence-analysis, biopython
---


# Question
Title: regex: samtools command to "refine" PacBio IsoSeq data?
I am working with preprocessed data (IsoSeq PacBio) and I cannot understand one of the steps:

```
samtools view -b -e [rq] >= 0.9 flnc.bam

```

What exactly was filtered from the bam file?

To obtain flnc they used the standard IsoSeq pipeline (lima -\> isoseq3 refine).

# Answer
> 2 votes
The `samtools view` command you've provided just uses a filter expression to discard records that do not meet the specified expression criteria. In this case, the expression is looking for alignments where the value of the 'rq' (read quality) tag is greater or equal to 0.9. The isoseq tags documentation says that the *read quality* of the alignment is really the *predicted accuracy of the polished isoform*. Note that the command produces BAM output, so be sure to redirect stdout or write to a file for example using:

```
samtools view -b -e '[rq] >= 0.9' -o out.bam flnc.bam

```

---
Tags: samtools, pacbio
---


# Question
Title: Low fraction reads in cell metric in snRNA-seq data
For some of our snRNA-seq samples we are finding low fraction of reads detected in single-nuclei rna-seq samples from cellranger, while the other metrics are perfect.

While I understand this could be caused by library preparation and abundance of ambient RNA. Just wondering if this behavior is expected for single-nuclei datasets. My confusion is specifically coming from the fact that 10X website also has such a sample that has low number of fraction reads in cells.

https://cf.10xgenomics.com/samples/cell-exp/7.0.0/5k\_mouse\_brain\_CNIK\_3pv3/5k\_mouse\_brain\_CNIK\_3pv3\_web\_summary.html

Would love to know if this behavior is expected or seen often who are more experienced in snRNA-seq. The metric documentation page is here

Here is snapshot from 10X report

# Answer
> 1 votes
Yes, "low fraction reads in cells" is more common with snRNA-seq from frozen samples, I wouldn't worry about it by itself (you're just throwing away some sequencing data). You should still check though that the other QC metrics are good. In particular, sometimes such samples have an unclear seperation of cells from empty droplets on the UMI rank plot and the consequent estimated cell count is off, requiring alternative approaches to distinguish cellular droplets from ambient RNA etc.

---
Tags: scrnaseq, transcriptome, cellranger
---


# Question
Title: What statistical test to apply for DE after CibersortX deconvolution?
*This question was also asked on Biostars*

I am running CibersortX in high-resolution mode (which yields estimates of gene values per sample).

After that, I want to perform DE between two conditions on the resulting gene estimates.

* What statistical test would one need to apply to perform DE after high-resolution mode de-convolution?
* If I use the output of DE for ranking purposes only (in pre-ranked GSEA), does it make sense to use a simple t-test?

The data is in TPM (CibersortX received TPMs as its input). I transform to log2 scale before t-test. There are only 7 samples available (4 control, 3 treatment)

Here are some genes which get the smallest p-value after t-test (but do not pass FDR)

```
ETFBKMT 0.08727244  0.2549699   0.1277664   0.000000    0.597810    0.4984627   0.4672444
AFF3    3.28551359  2.9969676   2.8793547   3.013827    2.483018    2.5070328   2.3119873
PDLIM1  4.26118568  4.0109785   3.9223976   3.856899    3.388990    3.5155023   3.2931998
MYOF    6.17754486  5.9286025   5.7452668   5.845756    5.276884    5.3981318   5.1251924
LINC02210 3.22922278    3.1785117   3.5162067   3.307364    2.736807    2.9091521   2.8052582
CD274   5.20308488  5.0063801   4.8715090   4.922224    4.440247    4.6192961   4.5082217

```

And these are the values before log

```
ETFBKMT 1.062360    1.193311    1.092601    1.000000    1.513417    1.412707    1.382466
AFF3    9.750753    7.983202    7.358209    8.077041    5.590659    5.684498    4.965666
PDLIM1  19.175412   16.122220   15.162099   14.489126   10.475814   11.435934   9.802840
MYOF    72.381287   60.909803   53.641094   57.510614   38.770421   42.169610   34.900902
LINC02210 9.377626  9.053727    11.441519   9.899558    6.665934    7.511766    6.989834
CD274   36.837031   32.141830   29.273209   30.320557   21.709387   24.578008   22.756735

```

And that is the a distribution of the expression values (after deconvolution) before log2, in one of the samples

```
 Min.   :    1.00
 1st Qu.:    5.01
 Median :   18.03
 Mean   :   89.08
 3rd Qu.:   47.82
 Max.   :34406.87

```

# Answer
> 0 votes
To specifically answer the question - a T-test is a parametric test which requires a normal distribution. So it's not a good idea in context to this data. However, the non-parametric equivalent of the T-test is the Wilcoxon and thats absolutely perfect and described below together with the entire rationale for parametrics versus non-parametrics. Basically, I agree with your approach, except I would simply shift to non-parametric statistics.

---

I noted the response on Biostars from the link recommends **DSeq2** and **EdgeR**. Both tests along with the additional test described below, along with the T-test, are used to identify whether each gene in the analysis is significantly up-, or down-regulated and can then be ranked. The underlying issue is the validity of the tests associated *P-values*.

Both DSeq2 and EdgeR have been notably criticised by Li *et al* 2022, in a paper entitled

> Li *et al* 2022 *Exaggerated false positives by popular differential expression methods when analyzing human population samples*" Genome Biology: 79

Li *et al* (2022) needs to be read carefully before making a choice. Li *et al* 2022 is recommending using a non-parametric test, via Wilcoxon, which is implemented in the method **NOIseq** from the R-Cran Bioc library (described below). Genome Biology has a current impact factor of 18 (exceptional). I was told that one member of the DSeq2 team produced a later paper advocating non-parametrics methods.

The underlying issue is the dependency on a parametric distribution for both DSeq2 and EdgeR. Basically, the central issue is the validity of parametric statistics (assumed distribution) versus non-parametrics statics (no assumptions about the distribution).

In parametric statistics outside DE, i.e. *in general*, when these are used numerous tests are performed to assess the fit of the data set to the parametric distribution assumed. For example, for ANOVA it's the Fmax test. For a T-test there's a number of tests for the normal distribution. If the data set fails the test, a multitude of transformations are possible usually aimed at normalising the variance.

From this perspective the concern is no comparative tests were devised for DSeq2 or EdgeR given the age of the papers DSeq2 is from 2014 Love *et al* (2014). If there's no test how could the data use alternative transformations to fit the assumptions? Again this is unusual given the age of papers because such as tests would be expected.

**Non-parameteric DEG**

The non-parametric method mentioned, **NOISeq**, is from the R/ **Bioc** package by Tarazona *et al* 2015, Nucleic Acids:43. This was the approach use by Li *et al* 2015 and compared against both DSeq2 and EdgeR. There are more recent examples, again I believe one of the DSeq2 team produced a non-parametric DEG as well.

The GitHub page for Li *et al* 2022, is https://github.com/xihuimeijing/DEGs\_Analysis\_FDR

---

I understand the criticism recently stated about Wilcoxon (non-parametrics in general) is their reduced power (response on Biostars). What happens is you get fewer significant 'hits' that are under or over-expressed genes. This is the basic trade-off between parametrics vs. non-parametrics. The key issue is if you've trapped your key results within a non-parametric test, *that is a very solid result*. If the equivalent results are trapped using EdgeR or DSeq2 there is always going to be question mark whether it's a false positive.

*When to use EdgeR/DSeq2*: if you are doing alternative statistical tests on the expression data which are directly comparable with under- and over-expression (e.g. ML and DL analysis) thats okay because we know Edge/DSeq2 gives false positives. When I've done this there's a 'centre ground' which can give conflicting result, its the central quartile type area, i.e. 25 - 75% region (to be honest it can be a bit outside that).

If you are singly dependent on EdgeR, DSeq2 or NOISeq it is better to be conservative in my opinion, i.e. fewer hits which you are more certain about.

**Final point** The final issue is simply how close to the DSeq2 and EdgeR data sets is the data you are working with, i.e. same species? What I would caution is these distributions have been applied broadly across species. The extent of the false positive rate must vary in my opinion.

---
Tags: rna-seq, differential-expression, deconvolution, t-test, cibersortx
---


# Question
Title: maxquant error: Unhandled Exception: System.IndexOutOfRangeException: Index was outside the bounds of the array
An error occurred while I was running the program maxquant (2.0.3):

Reading search engine results (SP)

Unhandled Exception: System.Exception: Exception during execution of external process: 1228404 Unhandled Exception: System.IndexOutOfRangeException: Index was outside the bounds of the array. at MsLib.Search.AndromedaPeptide.FromString (System.String line, BaseLibS.Ms.Data.Protein.ProteinSet proteinSet) \[0x000a9\] in :0 at MsLib.Search.AndromedaQuerySp.Add (System.Int32 count, System.String line, BaseLibS.Ms.Data.Protein.ProteinSet proteinSet, System.Int32 minLength) \[0x00001\] in :0

\[ERROR\] FATAL UNHANDLED EXCEPTION: System.IndexOutOfRangeException: Index was outside the bounds of the array. at MsLib.Search.AndromedaPeptide.FromString (System.String line, BaseLibS.Ms.Data.Protein.ProteinSet proteinSet) \[0x000a9\] in :0

how can i fix it, thank your very much!

# Answer
\****this was the xml file that i configured to maxquant. i changed modification and added the fastq and raw file. Other parameters were default. however, 'Unhandled Exception: System.IndexOutOfRangeException: Index was outside the bounds of the array' was reported when I run the command '*maquant mq.xml*'. but, When I only kept the N modification, the program ran successfully, but a small amount of peptides were identified.
> 1 votes

---
Tags: maxquant
---


title: SRA toolkit and SRA- prefix

text:
# Question
Title: SRA toolkit and SRA- prefix
I am trying to download raw sequence data from the SRA database using the SRA toolkit (v.3.0.2). Unfortunately, the authors of the article I am reading gave SRA numbers (i.e. with the SRA prefix, not the SRR). This does not seem to work with the prefetch command:

`prefetch --max-size 100000000 SRAxxxxxx -O downloads_SRAxxxxxx/`

I ended up going on NCBI and manually finding the relevant SRR number for each sample; however, I was wondering if I could have done it programmatically. I appreciate your help.

Edit: Thank you for your comments, here are a couple of examples:

SRA030738 ; SRA030736 ;

Searching for them on SRA I found the SRR numbers, respectively:

SRR135604 ; SRR135602.

The paper is Hasan et al, 2012: DOI: 10.1073/pnas.1207359109

# Answer
> 2 votes
Using Entrez Direct E-utils, run the following command for your two examples:

```
esearch -db SRA -query 'SRA030738 OR SRA030736' | efetch -format docsum | xtract -pattern DocumentSummary -element Runs | xtract -pattern Run -element Run@acc

```

(I know, it seems strange to run `xtract` twice, but the formatting on the plain xml output is dodgy. The first pass parsed out the unwanted `&gt` or `&lt` or `&quot` characters. The second pass worked on the re-formatted `stdout`)

# Answer
> 3 votes
@Max\_IT and @user3479780 Thanks for bringing this issue to our attention. It has been fixed in EDirect 19.5, which is now available on our public ftp site (See installation instructions here). Please update your installation and try again.

The problem was that using capitalized `-db "SRA"` bypassed normalization of the record. This was easily solved within EDirect by lower-casing the `-db` argument value to make it case-insensitive.

Normalization unescapes the XML-encoded ExpXml and Runs field contents to integrate them as regular XML within the docsum in situ, before passing them to the user. Thus:

```
efetch -db SRA -format docsum -id SRA030738,SRA030736 |
xtract -pattern DocumentSummary -element Run@acc

```

will now work properly, without the need for chaining two `xtract` commands together.

When creating structured records, NCBI internal software enforces conformation to the ASN.1 specification or XML document type definition. This is normally beneficial, blocking syntactically invalid records (and preventing new software bugs from silently corrupting production pipelines), but it limits the ability to add experimental structured data that is not (yet) in the spec.

Since `xtract` does not use a spec or DTD, it can find any element within the newly-unescaped contents by object name. (It still checks for properly-balanced opening and closing tags.)

Note also that you can skip the `esearch` step if you know the identifiers, and can instead pass them directly to `efetch` in an `-id` argument. For larger sets of identifiers, you can now pipe them from a file, one ID per line, without the need for an `epost` step.

EDirect has logic to split large lists into smaller chunks, sized appropriately for each database, and will internally loop through the fetch requests until all of the IDs have been processed, automatically retrying in the event of transient network or server failure.

---
Tags: sratoolkit
---


# Question
Title: Tidying MEGAN Taxonomy for Python/R Analysis
I am analysing some WGS data in MEGAN and would like to do some additional analysis in Python/R

I am having trouble Tidying the Taxonomic data in a format which would be conducive to this. Originally I tried exporting Taxon Path to count and using the Split function in Pandas to break up taxonomy into new columns as this worked well for KEGG Pathways. Due to varying specificity of taxonomic assignments this didn't work.

I've had more luck Exporting Data as a BIOM file format and converting it to TSV resulting in a file with the following column containing all taxonomic info (See below)

Ideally , I'd like to generate columns for; Domain, Phylum , Class , Order, Family and Species by accessing the features contained in Taxonomy. I'm unaware of any functions which would let me do that.

an idea of what I'm trying to achived is

| Domain | Phylum | Class | Order | Family | Species |
| --- | --- | --- | --- | --- | --- |
| Bacteria |  |  |  |  | uncultured bacterium |
| Bacteria | Bacteroidetes | Bacteroidia | Bacteroidaceae | Bacteroides | Bacteroides acidifaciens |

From

| Taxonomy |
| --- |
| d\__Bacteria; s\__uncultured bacterium |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides acidifaciens |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides caccae |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides cellulosilyticus |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides faecis |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides finegoldii |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides fragilis |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides intestinalis |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides salyersiae |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides stercoris |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides thetaiotaomicron |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides uniformis |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides xylanisolvens |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides cellulosilyticus CAG:158 |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__uncultured Bacteroides sp. |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides sp. 2\_2\_4 |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides sp. 3\_1\_23 |
| d\__Bacteria; p\__Bacteroidetes; c\__Bacteroidia; o\__Bacteroidales; f\__Bacteroidaceae; g\__Bacteroides; s\__Bacteroides sp. D2 |

# Answer
Save the following as a python script (e.g., named `splitty.py`). Make the script executable by doing `chmod +x splitty.py`. Then, if your `MEGAN` data are in a file called `taxonfile.csv`, run `splitty.py` by doing `./splitty.py taxonfile.csv`. Get help on the script by doing `splitty.py -h`.

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

"""
Splitty

Use this to split up a taxonomy table
"""

def argsparser():
    """
    argsparser
    """
    import argparse
    parser = argparse.ArgumentParser(prog="splitty",
                                    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
                                    description="Uset splitty to get taxon table.")
    parser.add_argument("table", help="Taxon table in semicolon separated format.")
    return parser

def parseline(line_in):
    """
    Parse each line_in as return a dict
    """
    transl = {"d": "Domain",
             "p": "Phylum",
             "c": "Class",
             "o": "Order",
             "f": "Family",
             "g": "Genus",
             "s": "Species"}
    return {transl[val[0]]: val[1] for val in [i.split("__") for i in line_in]}

def main(args):
    """
    The main routine.
    """
    import pandas as pd
    df_list = []
    with open(args.table, "r") as input_handle:
        for line in input_handle.readlines():
            linelist = parseline(list(filter(None, [i.strip() for i in line.rstrip().split(";")]))[1:])
            df_list.append(pd.DataFrame(linelist, index=["Taxon"]))
    print(pd.concat(df_list).to_csv(sep="\t", index=False))

if __name__ == "__main__":
    main(argsparser().parse_args())

```
> 3 votes

---
Tags: r, python, pandas, dataframe
---


title: Protein fold topology search web-server

text:
# Question
Title: Protein fold topology search web-server
I'm finding a web server (or standalone program) that allows the upload of a custom PDB for protein fold topology search.

Do you know any of them?

More details:

Actually I started off with Foldseek but have a question in mind:

how do I know if I have a novel fold? (because similarity search will often return partially similar hits).

Then, I tried to think of topology, as an illustration:

which simplifies structures to spatial arrangements of secondary structures (or topology).

I want to know if there's any database of protein fold topologies and tools for topology similarity search.

ADD COMMENT • link • edit •

# Answer
There are a few, but the main candidates for this would be :

1. CATH - http://www.cathdb.info/
2. SCOP - https://scop.mrc-lmb.cam.ac.uk/
3. Ptgl - http://ptgl.uni-frankfurt.de/

CATH and Ptgl have search services on their website, but not sure about SCOP.
> 1 votes

---
Tags: protein-structure
---


# Question
Title: WES variant calling with DNBSEQ-T7: technical quality assessment
I recently received whole exome sequencing samples that were sequenced on an MGI sequencing instrument, DNBSEQ-T7. I am interested in somatic variant-calling on the paired tumor-normal samples. As this isn't the canonical Illumina sequencer, an issue that I am concerned with is that the technical artifacts arising from the sequencer are not yet well characterised and this may lead to spurious artifacts in the variants.

I asked the company that did the sequencing for a panel of normals (PON) for the instrument to run Mutect2, but they simply replied that they use the standard GATK PON for their somatic variant calls.

How should I go about assessing whether the variants are "biased" in any way, as a result of the instrument's technical properties? Is there a way to quantify the magnitude of these technical errors from the VCF files?

# Answer
> 1 votes
Okay it seems the company has given you the VCF (maybe), you want to re-run it via Mutect2. Mutect2 is part of the GATK suite, it's an upgrade, and they are refusing to run the appropriate PON. I think this is more personnel management than a major technical hurdle.

I would reply to the tech and CC the technical director of the company explaining Mutect2 vs. regular GATK PON. You simply use a standard process of escalation within technical justification. Maybe send them the commend to demonstrate how simple it is:

```
     gatk Mutect2 \
     -R reference.fa \
     -I tumor.bam \
     -I normal.bam \
     -normal normal_sample_name \
     --germline-resource af-only-gnomad.vcf.gz \
     --panel-of-normals pon.vcf.gz \
     -O somatic.vcf.gz

```

The company will have to justify their position regarding PONs or accommodate your requirements either now or the future. There may simply be a technical misunderstanding between you and the tech.

You've paid for the service - thats the standard you expect. It's easier for them to assist you than all the agro of a complaints procedure and I'm pretty sure they'll help.

---
Tags: variant-calling, gatk, wes
---


# Question
Title: How do I force Cromwell / WDL to put input files in the same directory?
I have a cromwell task that, among other variables, takes these three variables:

```
   File Ref             # Reference genome fasta
   File RefFai          # Reference Genome index
   File RefDict         # Reference Dictionary

```

While all three files make it into the task's "inputs" directory, they end up in different sub-directories of inputs:

```
ME $ find inputs
inputs

inputs/-181926642
inputs/-181926642/Sample.vcf.bgz.tbi
inputs/-181926642/Sample.vcf.bgz
inputs/-1548908437
inputs/-1548908437/GCA_000001405.15_GRCh38_no_alt_short_headers_nonACTG_to_N.dict
inputs/44588251
inputs/44588251/GCA_000001405.15_GRCh38_no_alt_short_headers_nonACTG_to_N.fa
inputs/44588251/GCA_000001405.15_GRCh38_no_alt_short_headers_nonACTG_to_N.fa.fai

```

Which means that when I call GATK and pass it the path to the .fa file, it gets rather cranky that the .dict file isn't in the same directory.

How do I force cromwell to put all three files in the same directory?

# Answer
> 0 votes
My hacktastic fix for this is to use basename to get the name of each file, then do

```
refName=`basename "${Ref}"`
ln -s "${Ref}" "$refName"

```

for each file. This puts aliases to each file in the executions dir, and then I pass "refName" to GATK rather than "${Ref}".

Not a good solution, but a solution

# Answer
> 1 votes
as for wdl version 1.0 (other versions untested)

you can pack your reference and related files into a struct:

```
struct Reference {
    File fasta
    File dict
    File fai
}

```

You will pass the `Reference` to your tasks instead of `File`. Inside your workflow/config define:

```
Ref example = {
    "fasta" : "reference.fasta",
    "fai" : "reference.fasta.fai",
    "dict" : "reference.dict"
}

```

Inside the task you just use the file want use by `example.fasta`.

---
Tags: wdl
---


# Question
Title: Generating 3D cartesian coordinates of a peptide and computing van der Waal interactions (1-4 atoms) for a given set of dihedral angles?
I wanted to generate 3D Cartesian coordinates of a peptide sequence for fixed bond lengths and bond angles for a given set of backbone dihedral angles and side dihedral angles. I also wanted to compute non-bonded van der Waal interactions (1-4 atoms) for each of the peptide generated from the dihedral angles. I initially tried frag builder module in python although it works fine. There are few issues,

1. I could not fix bond length and bond angle of my own choice.
2. It generates PDB file and save it inside the hard disk. So I need to read it again and parse it to compute non-bonded van der Waal interactions (1-4 atom). It consumes lot of time
3. Also I find fragbuilder is not that fast in generating peptides.
4. Reading , parsing and computing the interactions takes around 30 mins to complete the task for 200,000 3D ensemble of structures of given peptide sequence of length 18.

Kindly suggest me a better module in python which does the job faster and also give reference to it.

# Answer
The PDB module in Biopython is worth checking out for obtaining angles, distances, torsion angles. It is described from p208 at http://biopython.org/DIST/docs/tutorial/Tutorial.pdf. It doesn't do van der Waal interactions.

For example,

```
from Bio.PDB.PDBParser import PDBParser
from Bio.PDB.Chain import Chain
from Bio.PDB.internal_coords import *
from Bio.PDB.PICIO import write_PIC, read_PIC, read_PIC_seq
from Bio.PDB.ic_rebuild import write_PDB, IC_duplicate, structure_rebuild_test
from Bio.PDB.SCADIO import write_SCAD
from Bio.PDB.PDBIO import PDBIO

parser = PDBParser(PERMISSIVE=1)
myProtein = parser.get_structure("1a8o", "1A8O.pdb")
myChain = myProtein[0]["A"]
IC_Chain.MaxPeptideBond = 4.0
myChain.internal_coord = None
myChain.atom_to_internal_coordinates(verbose=True)
r2 = myChain.child_list[1]
r2ic = r2.internal_coord
r2chi1 = r2ic.get_angle("chi1")

```
> 1 votes

---
Tags: python, protein-structure, pdb, molecular-dynamics
---


# Question
Title: When it is stated in your consumer genetics testing report that you have x% (lets say 1.3 percent...) Neanderthal DNA, what is actually meant?
*This question has also been asked on Biology SE*

When a raw number like 1.3 percent is stated, it is confusing without context,and can be misinterpreted. So does 1.3 percent mean that 1.3 percent of the letters/nucleotides/bases of my genome consists of established (likely/high confidence) Neanderthal markers or SNPS? Does it mean 1.3 percent of my entire genome consists of these Neanderthal markers, i.e 1.3 percent of 3.1 billion or 6.2 billion base pairs, or is it 1.3 percent of protein coding regions, or exons?\**

How do they obtain these markers for Neanderthal ancestry to begin with? I know the Neanderthal genome has been sequenced from bone fragments. Do they compare SNPs and variable regions in the genomes between Neanderthal and Modern Human population genomes to see the relative frequencies of these variations/SNPs. Then likely Neanderthal markers are flagged by statistical models and algorithms..for example combinations of SNPs that appear infrequently or hardly at all in the reference Sapien population genome but occur frequently in the Neanderthal reference genome..that marks these SNPs/variable regions as of likely Neanderthal origin? Could you please describe some of these statistical methods and algorithms? Can anything else explain the "markers", or SNPs other than interbreeding with Neanderthal (such as mutations)? Please illustrate using some hypothetical alignments like AAGGGT..vs GAGGGT..Please give a clear, concise and authoritative answer of high confidence. It is okay to be detailed and technical, although building up to complexity so that it can be understood with relatively little familiarity with population genetics and Mathematics.

# Answer
I ran a basic search of 23andme and it looks like they use known variant loci from the Neanderthal genome. Also, they expect a 2-2.5% overlap. See:

1. https://customercare.23andme.com/hc/en-us/articles/212873707-Neanderthal-Ancestry-Report-Basics
2. https://blog.23andme.com/articles/new-neanderthal-report (Section "The Science behind the report")
3. https://permalinks.23andme.com/pdf/23-05\_neanderthal\_ancestry\_inference.pdf (Eric Durand's white paper on the technical details)
> 2 votes

# Answer
Based on reading a similar post it seems like if the report says I am 1.3 percent Neanderthal then it does NOT mean 1.3% of my base pairs or SNPs of my genome (or exons or protein coding region or total number of my SNPs) has a "Neanderthal" label inferred, rather some kind of statistical inference is made as a whole.

The author suggests when it comes to 23AndMe that my genome is segmented into chunks of a chromosome (containing a set of SNPs, 23AndMe does genotyping only reads parts of chromosome where bases are known to vary.)

Technical details could be referred to the following sites: https://academic.oup.com/bib/article/19/5/765/3038466

https://www.cd-genomics.com/the-principles-and-workflow-of-snp-microarray.html

https://www.illumina.com/techniques/popular-applications/genotyping/snp-snv-genotyping.html

Since 23andMe is using v5 genotyping chip, the fraction of Neanderthal variant could be calculated. Moreover, the basis of v5 genotyping chip is related to the concept that DNA-based polymorphisms (SNPs or SNVs) which differ substantially in frequency among human populations can be employed to infer ancestry. With 23andMe, an average of 236 Neanderthal SNPs/SNVs could be obtained with the analysis based on a selected group of consented customers.

And these chunks are fed into a Machine learning classification algorithm like Support vector machine (SVM). The Machine learning algorithm was trained chunks of genome of reference population sets with chunks of known ancestry (like "Scandinavian" or "South European" or "Sub Saharan African").

Then the classification algorithm makes a multilabel prediction... lets say the algorithm predicts that a particular segment is 95 percent likely to be human and 5 percent Neanderthal..and maybe the next segment has a 2% probability of being human and 98 percent Neanderthal..Then the algorithm looks for segments of high confidence (\>90 percent or 70 percent whatever threshold you set) and then add them up.

For example if in a toy example genome has 9 segments of the same size has a greater than 90 percent probability of being estimated "human", and 1 segment with greater than 90 percent of being estimated Neanderthal by the classification algorithm than I would be considered to be 10 percent Neanderthal...

The author suggests this is the method used to infer ancestry like "Norwegian" and "South European" and "Sub Saharan African" by 23And Me..but will it be similar for inferring Neanderthal ancestry?

Calculating ancestry percentage is not seem as simple as counting up number of "Neanderthal" SNPs/Markers and dividing by your total SNPs or base pairs in your genome/exon...Rather an inference of ancestry mix needs to be made based on reading SNPs on segments of chromosomes and comparing it to reference data sets via classification algorithms.
> 2 votes

# Answer
> How do they obtain these markers for Neanderthal ancestry to begin with? I know the Neanderthal genome has been sequenced from bone fragments

1. Since ancient nuclear DNA is scarce, enough preserved molecules are needed.
2. The level of scarcity about ancient nuclear DNA : often, 1% or fewer of the total DNA fragments from an extraction among soil microbes and other sources
3. After extraction and sequencing, distinguish the genuine DNA fragments from those of modern contaminants and from those of microorganisms that are co-extracted from the sample
4. Since ancient DNA has characteristic damage patterns, e.g. a high degree of cytosine deaminations on the five prime ends, and at the three prime end, there is a high frequency of Gs changing to As.
5. Match these ancient DNA fragments to a map of the complete human genome reference sequence; this allowed the book of the target's genome to be assembled from the millions of random sequence fragments recovered from the bones
6. The more times we had a read confirmed by multiple fragments, the more confident the sequencing could be in distinguishing the true DNA sequence, it is the sequencing depth issue
7. Then genotype samples to discover, filter, and analyze variants.
> 1 votes

---
Tags: snp, variation
---


# Question
Title: Unable to install an R Bioconductor package on Apple M1/M2
I have an Apple latest arm64 machine. Running the latest R 4.3.0 version. However, I'm unable to install an R package - `S4Arrays`.

Doesn't look like there is an arm64 build for this package. Anybody know how to resolve it?

# Answer
> 2 votes
The easiest solution is to install the x86 version of R in your M1/M2 Mac here. It will be then emulated by the Rosetta. After that you will have no problems with Bioconductor packages

---
Tags: r, bioconductor
---


# Question
Title: How to get the nucleotide sequence for protein NCBI accession code, programatically?
It is a lot easier in my work to do blast searches against protein databases when I can, since I work with distantly related organisms. But sometimes, I need to work with the actual nucleotide sequences. Until now, I've often done it manually, but it is getting unfeasible due to increasing data amounts.

The NCBI website is structured such that I suspect automating the process (e.g. input a list of protein accession numbers, get out a fasta file with the corresponding nucleotide seqs) should be possible. At the same time, I am not quite sure what tools to use (biopython?), and if a tool already exists, I would rather not reinvent the wheel. Would you happen to know of anything like this? Thank you for your time!

**Update and WIP:**

I have tried Entrez Utils. In part it works, but I am having some issues scaling it up. For example,

`elink -db protein -target nuccore -id WP_244350345.1 | efetch -format fasta_cds_na | awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);} END {printf("\n");}' | grep WP_244350345.1 -A 1`
works. (elinks to connect the protein accession WP\_244350345.1 with the nucleotide database | efetch to grab the nucleotide sequence, which is a bit of a pain because it grabs the entire genome rather than the cds I want | so I need a bit of cleaning, including using awk to turn the resulting multi-line fasta to single-line fasta | and grep to actually get the cds sequence associated with WP\_244350345.1. As I said, it works, and it results in:

```
>lcl|NZ_JAFIRA010000020.1_cds_WP_244350345.1_35 [gene=bchM] [locus_tag=JX360_RS09120] [protein=magnesium protoporphyrin IX methyltransferase] [protein_id=WP_244350345.1] [location=35695..36360] [gbkey=CDS]
ATGCAAGAGAAAAGCATCGTCCGCGACTATTTCAACTCCATTGGCTTCGACCGTTGGCGGCGCATCTACGGAGAGGATGAGGTCAATTTTGTCCAAAAAGACATTCGCTCCGGCCATGCCCGTACCGTAGCCACTGTGTTGGAATGGCTGGGGGATCCCTCAGGTTTGCAGATTTGTGATGCGGGGTGTGGGGTCGGTAGTCTTAGTTTGCCGTTGGCTGCCAAAGGGGCCCGCGTCTTCGCCAGCGACATTTCCGAACAAATGGTGAACGAAGCCCGCCGTCGCCAGCAAATCCAGTTGGGATCCACCGAGAACCCCCAATTTCGGGTCTCGGACTTAGAGGAACTAAGCGGAGAATATGATGCCGTCATTTGCCTGGATGTGATGATCCACTACCCAGAAGCCGATGCCCTGCGCATGCTGGAACATCTAACCCGCCTGGCCCGCTCTCGCTTGATCTTCAGTTTTGCCCCCAAAACCCCTCTGTTGACGTTGCTAAAAAAAGTGGGAGAGTTCTTCCCCGGCCCCAGCAAAGCCACCCGTGCCTACCAACACCGAGAAGCGATTTTGGTGAGCAAATTAACCGAATTGGGCTGGAACGTTCAGCAACGACAGACCATTCGTAGCCGCTTTTACTTTGCCTGGATTCTGGATTTGGCCCGGTGA

```

*However*, trying to read multiple accession numbers from a file, in a loop, is not working. My current code (essentially the same as before, just replacing WP\_244350345.1 with ${LINE}:

```
while IFS='' read -r LINE || [ -n "${LINE}" ]; do
echo "${LINE}"
elink -db protein -target nuccore -id "${LINE}" | efetch -format fasta_cds_na | awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);}  END {printf("\n");}' | grep "${LINE}" -A 1
done < "results.txt" | more

```

The `echo "${LINE}"` part works fine, the rest does not. It has also worked perfectly fine with other parts of Entrez Utils (see below), so I am confused and I would appreciate some help.

```
while IFS='' read -r LINE || [ -n "${LINE}" ]; do
organism=$(elink -db protein -target taxonomy -id "${LINE}" | efetch | tr --delete '\n' | sed 's/    species, /_/g'| sed 's/^.*\. //g' | sed 's/ /_/g')
echo "${LINE}"
echo "$organism"
done < results.txt | more

```

results in:

```
WP_244350345.1
Thermostichus_vulcanus_cyanobacteria
WP_073547903.1
Chroogloeocystis_siderophila_cyanobacteria
WP_088429095.1
Halomicronema_hongdechloris_cyanobacteria
...

```

# Answer
Yes its easy using Biopython.

What I'm requesting is the sequence must be \>9500 otherwise it's discarded. In your case this value will be MUCH bigger, because these are bacterial genomes.

```
from Bio import Entrez
from Bio import SeqIO, AlignIO

def getgenbank (smallGenDict):
    bacteriadb = []
    for bacteria in smallGenDict.keys():
        Entrez.email = "name@email.com"
        with Entrez.efetch(
                db="nucleotide", rettype="gb", retmode="text", id=bacteria
                ) as handle:
            seq_record = SeqIO.read(handle, "gb")
            # record = get_CDS(seq_record, bacteria, smallGenDict.keys)
            record = seq_record # better to remove the above function because its outside the question
            if len(record.seq) > 9500:
                bacteriadb.append(record)
    return bacteriadb

def get_CDS(record, bacteria, smallGenDict):
    for feature in record.features:
        # custom manipulation of input sequences
        # better to forget this bit of code
        # this whole part of the code has been skipped
        # its capturing the coding sequence (not relevant)
return record

def main(smallGenDict):
    db = getgenbank(smallGenDict)
    SeqIO.write(db, '/Users/username/Desktop/myoutput.fasta', "fasta")

if __name__ == '__main__':
    smallGenDict = {'ABxxxxxx':'AGACTGATA', 'ABxxxxx2:'GAGATGA'} # better to automate this bit (different question)
    main(smallGenDict)

```

Thats it.

**Notes** In the code below, the file smallGenDict is a dictionary of Genbank accession No.s and their sequence. `{'ABxxxxxx':'AGACTGATA', 'ABxxxxx2:'GAGATGA'}` This better written here as a list `['ABxxxxxx','ABxxxxx2' ...]` say call it `smallGenList`.

Where it says `smallGenDict.keys():` just replace it with `smallGenList` and remove the method `.keys()`, otherwise the code will through a bug. The function `get_CDS` is processing the sequence file to remove the non-coding regions. I've not included it because its specific stuff.
> 1 votes

# Answer
You can use R Rentrez, in the link you can find a tutorial for the package usage.
> 0 votes

---
Tags: proteins, gene, ncbi
---


# Question
Title: CNV Kit ` from . import commands ImportError: cannot import name 'commands' from '__main__'`
I am trying to run some code for my colleague in bash

```
/path/to/cnvkit.py batch /path/to/my/folder/with/bams/*.bam \
  --normal --targets ${bed_file} \
  --fasta path/to/my/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta \
  --output-reference /path/to/my/CD_BATCH1_reference.cnn \
  --output-dir /path/to/my/Group_1 --scatter

```

however, I keep getting this truly peculiar error

```
Traceback (most recent call last):
  File "/path/to/cnvkit.py", line 4, in <module>
    from . import commands
ImportError: cannot import name 'commands' from '__main__' (/path/to/cnvkit.py)

```

**Attempt 1**

I checked if the path was correct and multiple times I found that it was correct. It did not change the error.

**Attempt 2**

Perhaps I need to `chmod +x *` everything in the CNV kit directory and that will help. It did not change the error.

Why is this error occurring? How do I solve this issue?

# Answer
> 1 votes
Could you download it from `conda` and try again? Its here, https://docs.conda.io/en/main/miniconda.html

Then,

```
conda create --name myenv
conda activate myenv
conda install -c bioconda cnvkit
cnvkit batch /path/to/my/folder/with/bams/*.bam --normal --targets ${bed_file} --fasta path/to/my/resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta --output-reference /path/to/my/CD_BATCH1_reference.cnn --output-dir /path/to/my/Group_1 --scatter

```

The path to the file is not used, the command is simply `cnvkit` This avoids all the privilege stuff and ensures you're running a virtual environment to protect your system Python.

---

To address the comments. "conda" is slow and doesn't ever load.

There is an alternative to conda called **mamba**, this is fast and designed to address this exact issue. I have never needed to use it so the following is a rough guide only.

Download and install via

```
conda create --name myenv
conda activate myenv
conda install -c conda-forge mamba

```

The documentation is here, https://mamba.readthedocs.io/en/latest/user\_guide/mamba.html

I've never used it but the idea is ...

```
mamba create -n mambaenv <list of packages>
mamba activate mambaenv

```

The list of packages could be left blank for now, then,

```
mamba install "Python>=3.10.0" numpy pandas pathlib
mamba install cnvkit # SEE NOTE

```

Can also do

```
mamba install -n mambaenv cnvkit

```

**NOTE** cnvkit exists strictly in -c bioconda. You do not seem to search and install via channels in mamba

You search for a package and its dependencies via

```
mamba repoquery search cnvkit
mamba repoquery depends -a cnvkit

```

It is important to place a virtual environment between the system Python and package Python (don't want a system corruption). More to the question, this would solve all the permission problems without changing your system in anyway (kinda important).

---
Tags: bash, cnv
---


# Question
Title: NCBI Blast databases - cannot find some databases
I'm working with a system that does local BLAST searches using some preconfigured BLAST databases. From what I see, these were updated within the last year, but I cannot find any current databases available that have those same names. One specific set is a non-redundant set of human protein or nucleotide sequences. The databases are named human\_nr and human\_nt (with various extensions as per the BLAST db format). There are others names up\_sp\_human (Uniprot-Swissprot human?) and up\_tbl\_human (Uniprot-Tremble human?), but a cannot find any similar databases on the NCBI site or other sites.

Where else should I look for such preconfigured databases?

# Answer
I don't think there is a human\_nr or human\_nt database for NCBI.

An equivalent database is accessible via the standard `nr/nt` database and simply restrict the organism parameter to `human` or to taxa id `9606` . This essentially will do the same thing.

Or maybe this sort of thing?,

```
blastdbcmd -db nr -entry all -outfmt "%g %T" | awk ' { if ($2 == 9606) { print $1 } } ' > human_flat-list.db

```
> 3 votes

---
Tags: phylogenetics, phylogeny, blast, public-databases, ncbi
---


# Question
Title: Error while running computeMatrix command in Deeptools
I am trying to get computeMatrix for bigwig file in specific genomic region using deeptools.
Below is the code I am using

```
computeMatrix reference-point --referencePoint TSS \
-b 1000 -a 1000 \
-R ~/Desktop/ATAC/ATAC/Inducible_elements_Greenberg.bed \
-S ~/Desktop/ATAC/Control.mRp.clN.bigWig \
--skipZeros \
-o ~/Desktop/ATAC/matrix_controlmRPATACBasalcomputematrix.gz \

```

But I get the following error

```
File "/Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/site-packages/numpy/__init__.py", line 305, in __getattr__
    raise AttributeError(__former_attrs__[attr])
AttributeError: module 'numpy' has no attribute 'float'.
`np.float` was a deprecated alias for the builtin `float`. To avoid this error in existing code, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:
    https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations. Did you mean: 'cfloat'?

```

These are the packages in my Conda environment

```
# packages in environment at /Users/abhinavsoni/anaconda3/envs/ATAC_Seq:
#
# Name                    Version                   Build  Channel
bzip2                     1.0.8                h3422bc3_4    conda-forge
ca-certificates           2023.5.7             hf0a4a13_0    conda-forge
libblas                   3.9.0           16_osxarm64_openblas    conda-forge
libcblas                  3.9.0           16_osxarm64_openblas    conda-forge
libcxx                    16.0.3               h4653b0c_0    conda-forge
libffi                    3.4.2                h3422bc3_5    conda-forge
libgfortran               5.0.0           12_2_0_hd922786_31    conda-forge
libgfortran5              12.2.0              h0eea778_31    conda-forge
liblapack                 3.9.0           16_osxarm64_openblas    conda-forge
libopenblas               0.3.21          openmp_hc731615_3    conda-forge
libsqlite                 3.41.2               hb31c410_1    conda-forge
libzlib                   1.2.13               h03a7124_4    conda-forge
llvm-openmp               16.0.3               h1c12783_0    conda-forge
ncurses                   6.3                  h07bb92c_1    conda-forge
numpy                     1.23.0          py310h0a343b5_0    conda-forge
openssl                   3.1.0                h53f4e23_3    conda-forge
pip                       23.1.2             pyhd8ed1ab_0    conda-forge
python                    3.10.11         h3ba56d0_0_cpython    conda-forge
python_abi                3.10                    3_cp310    conda-forge
readline                  8.2                  h92ec313_1    conda-forge
setuptools                67.7.2             pyhd8ed1ab_0    conda-forge
tk                        8.6.12               he1e0b03_0    conda-forge
tzdata                    2023c                h71feb2d_0    conda-forge
wheel                     0.40.0             pyhd8ed1ab_0    conda-forge
xz                        5.2.6                h57fd34a_0    conda-forge

```

What could be its possible work around?

# Answer
It worked, I installed a the newer version of deeptools, from GitHub:

https://github.com/deeptools/deepTools
> 2 votes

---
Tags: python, chip-seq, atac-seq, deeptools
---


title: Where do I get a large reference VCF?

text:
# Question
Title: Where do I get a large reference VCF?
I would like to download a large .vcf file containing many (hundreds or thousands) of samples. Ideally, I would download different population-specific .vcf files, but the ability to sort/filter by ancestry group is fine. Where do I get such a file? I prefer GRCh37 for consistency with other files I'm using.

---

Update: I tried the HGDP file--I'm trying the 100 genomes one now, in addition to trying http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data\_collections/1000\_genomes\_project/release/20190312\_biallelic\_SNV\_and\_INDEL/ now.

I downloaded the file, but when I opened it it seemed to have no samples (instead of having many samples). Here is a full row:

```
 1  10145   .   AAC A   151.43  PASS    AC=4;AF=0.0127389;AN=314;BaseQRankSum=0.72;ClippingRankSum=0;DP=24144;ExcessHet=160.275;FS=0;InbreedingCoeff=-0.1558;MLEAC=5;MLEAF=0.002604;MQ=1.63;MQRankSum=0.764;NEGATIVE_TRAIN_SITE;QD=3.52;ReadPosRankSum=0.085;SOR=0.458;VQSLOD=-0.9391;culprit=QD;VQSRMODE=INDEL;NS=157;ExcHet=0

```

*This question has also been asked on Biostars*

# Answer
Solution: Use https://bochet.gcc.biostat.washington.edu/beagle/1000\_Genomes\_phase3\_v5a/b37.vcf/, then to get the ancestry information (i.e., the ancestry group of each sample, e.g., "HG00566" or "NA19758") use https://www.internationalgenome.org/data-portal/sample.
> 1 votes

# Answer
For anyone who arrives here looking for a large reference VCF but *not* wanting GRCh37, you might consider the HGDP+1KG subset of the Genome Aggregation Database's dataset. You can find a VCF file per chromosome aligned, I believe, to GRCh38, here.
> 1 votes

---
Tags: vcf, genome, 1000genomes
---


# Question
Title: How to convert fasta file to tab delimited file
I have a fasta file like

```
>sample 1 gene 1
atgc
>sample 1 gene 2
atgc
>sample 2 gene 1
atgc

```

I want to get the following output, with one break between the header and the sequence.

```
>sample 1 gene 1   atgc
>sample 1 gene 2   atgc
>sample 2 gene 1   atgc

```

# Answer
> 13 votes
If you have multi-line fasta files, as is very common, you can use these scripts1 to convert between fasta and tbl (sequence\_name `<TAB>` sequence) format:

* FastaToTbl

```
#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
      if (NR>1)
                     printf "\n%s\t", substr($0,2,length($0)-1)
      else
       printf "%s\t", substr($0,2,length($0)-1)
            else
                    printf "%s", $0
}END{printf "\n"}

```
* TblToFasta

```
#! /usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld

     if (fld == NF-1)
      {
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls)
  {
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}

```

Save those in your `$PATH`, make them executable, and you can then do:

```
$ cat file.fa
>sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
>sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC
$ FastaToTbl file.fa
sequence1   ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGCTAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
sequence2   GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGCATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

```

And, to get the Fasta back:

```
$ FastaToTbl file.fa | TblToFasta
>sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
>sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

```

This can be a very useful trick when searching a fasta file for a string:

```
TblToFasta file.fa | grep 'foo' | FastaToTbl

```

---

If you really want to keep the leading `>` of the header (which doesn't seem very useful), you could do something like this:

```
$ perl -0pe 's/\n//g; s/.>/\n>/g; s/$/\n/;' file.fa
>sequence1 ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGCTAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATAC
>sequence2 GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGCATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

```

But that will read the entire file into memory. If that's an issue, add an empty line between each fasta record, and then use perl's paragraph mode to process each "paragraph" (sequence) at a time:

```
perl -pe  's/>/\n>/' file.fa | perl -00pe 's/\n//g; s/.>/\n>/g; s/$/\n/;'

```

---

1Credit to Josep Abril who wrote these scripts more than a decade ago.

# Answer
> 12 votes
There is a very simple BioPython solution, that is minimal, readable, and handles multi-line fasta:

```
from Bio import SeqIO

for record in SeqIO.parse('example.fa', 'fasta'):
    print('>{}\t{}'.format(record.description, record.seq))

```

# Answer
> 9 votes
assuming there is only one sequence line per record, use paste with two 'stdin'

```
cat your.fasta | paste - -

```

# Answer
> 8 votes
You can use these commands:

```
perl -pe 's/>(.*)/>\1\t/g; s/\n//g; s/>/\n>/g' sequences.fa | grep -v '^$'

```

Explanation:

1. Append a tab to every header line
2. Join all lines
3. Split the single obtained line by the '\>' character
4. Remove the empty line (the first line is empty due to the fact that '\>' is the first character of the FASTA file)

# Answer
> 7 votes
A very useful tool for this kind of data manipulation is bioawk:

```
$ bioawk -c fastx '{print ">"$name" "$comment"\t"$seq}' test.fa
>sample 1 gene 1    atgc
>sample 1 gene 2    atgc
>sample 2 gene 1    atgc

```

bioawk is based on awk, with added parsing capabilities. Here, we tell that the format is fasta or fastq with `-c fastx`, and this makes the `$name` (between "\>" and the first blank character), `$comment` (after the first blank character) and `$seq` (the sequence, in one line) variables available within awk instructions.

See for instance this answer for another use case.

# Answer
> 4 votes
Where possible, I recommend using a dedicated parsing library, rather than hacking a parser together: as you can see in the other answers, parsing even simple formats gets complex pretty quickly if you value correctness.

Here’s a small R script that does what we need, using ‘seqinr’:

```
#!/usr/bin/env Rscript
suppressPackageStartupMessages(library(seqinr))
parsed = read.fasta(file('stdin'), as.string = TRUE)
table = data.frame(unlist(parsed), row.names = sapply(parsed, attr, 'Annot'))
write.table(table, stdout(), sep = '\t', quote = FALSE, col.names = FALSE)

```

Save it as `fasta-to-tsv`, make it executable, and use it as follows:

```
fasta-to-tsv < input.fasta > output.tsv

```

Equivalent code of similar length can be written in Python or Perl.

# Answer
> 3 votes
This could be easily done by seqkit fx2tab

```
seqkit fx2tab seq.fa

```

However, seqkit will not print the "greater than" symbol ("\>"). If you do need the symbol:

```
seqkit fx2tab seq.fa | sed 's/^/>/g'

```

# Answer
> 2 votes
Remove empty records (description without sequence):

```
awk '$2{print RS}$2' FS='\n' RS=\> ORS= f1.fa > f2.fa

```

Remove blank lines:

```
sed '/^$/d' f2.fa > f3.fa

```

Convert multi-line fasta to single-line fasta:

```
awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);} END {printf("\n");}' f3.fa > f4.fa

```

Finally, @Pierre solution:

```
cat f4.fa | paste - - > f.txt

```

# Answer
> 2 votes
In cases where there is no sequence wrapping and each sequence occupies only a single line, the following shell command is probably going to be fastest, easiest, and most convenient.

```
paste - - < your.fasta > your.new.fasta

```

# Answer
> 2 votes
Just suggested improvements to the awk scripts in @terdon's answer, using any POSIX awk:

```
$ cat FastaToTbl
#!/usr/bin/env bash

awk -v OFS='\t' '
    {
        if ( /^>/ ) { out = (NR>1 ? ORS : "") substr($0,2) OFS }
        else        { out = $0 }
        printf "%s", out
    }
    END { print "" }
' "${@:--}"

```

```
$ cat TblToFasta
#!/usr/bin/env bash

awk -F'\t' '
    {
        gsub(/.{60}/,"&"ORS,$2)
        sub(ORS"$","",$2)
        print ">" $1 ORS $2
    }
' "${@:--}"

```

```
$ ./FastaToTbl file.fa
sequence1       ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGCTAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
sequence2       GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGCATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

```

```
$ ./FastaToTbl file.fa | ./TblToFasta
>sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
>sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

```

`FastaToTbl` will actually work in any awk but `TblToFasta` requires a POSIX awk for support of regexp intervals like `{60}`. If you have a very old, pre-POSIX awk that doesn't support regexp intervals then **get a new awk** but if that's impossible for some reason then change `TblToFasta` to the following and then it'll also work in any awk:

```
$ cat TblToFasta
#!/usr/bin/env bash

awk -F'\t' '
    BEGIN {
        dots = sprintf("%"60"s","")
        gsub(/ /,".",dots)
    }
    {
        gsub(dots,"&"ORS,$2)
        sub(ORS"$","",$2)
        print ">" $1 ORS $2
    }
' "${@:--}"

```

# Answer
> 1 votes
This is an old post I have noticed and there are many offered solutions. Since it’s a frequently asked question, I thought it’s worth for me to mention that there is an overlooked tool set which contains a stand-alone program called `faToTab` in addition to many other useful bioinformatics tools.

`faToTab inputFile.fasta outFileFasta_tab.txt`

It’s a gold-chest in my opinion. Here are the links to utilities folder and details: Description and Download instructions \- Binaries by machine type \- Link to the github page.

Anaconda installation is:

```
conda install -c bioconda ucsc-fatotab
conda install -c bioconda/label/cf201901 ucsc-fatotab

```

# Answer
> 1 votes
In python I would do:

```
#...Suppose you have header information and
#...sequences stored in lists

header_info1 = [elements]
header_info2 = [elements]
sequences = [sequences]
index = enumerate(sequences)
table = open(pathtoyourfile.tsv,'w+')
for h1,h2,s,i in zip(header_info1,header_info2,sequences,index):
    table.write(f">{h1}\t{h2}\t{s}\n")
    if i+1==len(sequences):
       table.write(f">{h1}\t{h2}\t{s}")
table.close()

```

So basically I use f strings and I iterate over these three vectors that are of the same length. At the end of the iteration, I remove new line(\n) since it won't be needed anymore since you have not to write anything further.

---
Tags: fasta, format-conversion
---


# Question
Title: Software for microbial profiling from 16S rRNA gene sequence
I have hundreds of GBs metagenomic 16S rRNA gene sequence data. I want to do **microbiome composition profiling** (with relative abundance) from the data. Also after that, I will do functional profiling (**gene and metabolome level**).

Can you please suggest me the **best software(s)** for these two tasks?

Thank You

# Answer
I tried dada2 and is not bad (if you know R).

QIIME2 is also an option. Many other are available, the choice might also depend on your sample and your exact question.

For functional profiling you may try picrust2
> 2 votes

# Answer
For comprehensive microbial profiling of your 16S rRNA gene sequence data from Illumina MiSeq, I recommend using the web app Microbioma16S (www.microbioma16s.it). This powerful tool offers an end-to-end pipeline that covers everything from preprocessing raw data to conducting alpha and beta diversity analysis. With Microbioma16S, you can explore interactive results of taxonomy associated with features, taxa relative abundances, and perform comparisons between sample groups. It's a user-friendly solution for your microbial profiling needs.
> 1 votes

# Answer
The best thing to do is a phylogenetic tree in my opinion.

The trendy approach and the one deemed acceptable in publication is deep learning for precisely this application, such as the pipeline published here. This approach is attractive in 16S NGS metagenomics classification because it uses advanced approaches in deep learning, notably convolution neural networks (CNN) and deep belief networks (DBN). I don't really understand DBN, but I suspect it is a variant of recurrent neural networks (RNN) such as long-short term memory (LSTM).

Having said all that the accuracy is low, viz.

> For instance, at the genus level, both CNN and DBN reached 91.3% of accuracy with AMP short-reads, whereas RDP classifier obtained 83.8% with the same data.

91% accuarcy for a deep learning neural network is minimal, but it will be quick and will do precisely what you want with little fuss. A phylogenetic tree on the otherhand can take several hundred hours of total processor time. The thing about a 16S tree is that it is 100% accurate, because you can manually identify "grey zones" to assign their classification.

You can retrain CNN models and in future this will likely resolve the accuracy issue. However, I don't know about about DBN and I don't think RSS can retrained, more specifically I don't know how to do it.
> 0 votes

---
Tags: software-recommendation, metagenome, multi-omics, systems-biology
---


# Question
Title: Selection analysis based on each leaves as foreground with FastCodeML
I am using positive selection analysis for each leaves in a phylogenetic tree.

In the routine calculation via `codeml`, I will be running model zero for all leaves and then model 2 for each leaves as foreground (marking leaves as #1 in phylogenetic tree), then calculate selection using a likelihood ratio test.

Is there a method for using `FastCodeML` for each leaves as foreground instead of the branch?

# Answer
First and foremost, the fastcodeml guys didn't replace codeml via refactorising all the functionality of codeml, they refactorised one particular calculation and hence the confusion. fastcodeml isn't a particularly name helpful in truth, because it is not a faster version of all the functionality codeml, its just one of the calculations.

In the authors defence they (notably) improved the runtime of the most intensive calculation in codeml and thats their point - its proof of concept.

Thus, basically no fastcodeml can not be used in this context. fastcodeml runs the branch-site model of codeml only.

What the OP describes isn't how the branch-sites model works, they are describing the branch model.

The branch-site null model sets dN/dS to 1, and the LRT compares dN/dS likelihood between the designated branch(es) between dN/dS fixed to neutral against allowing dN/dS to be any value. Thus it's measuring whether the branch (branches) deviate from neutrality. What the OP describes via, model 0 one rate for the tree versus model 2 meaning two or more rates for the tree (basically two rates) and that is a different analysis and it is called a branch model.

---

If fastcodeml can be shoehorned to perform sites and branch-models that would be good to know.

---

There is an issue about 'tip selection' but the key answer concerns the basic difference between codeml and fastcodeml.
> 0 votes

---
Tags: phylogenetics, phylogeny, positive-selection
---


# Question
Title: Create bed file where each column is the number of intersections with another bed file
For example I have 3 bed files:

window.bed:

```
chr1    0   1000
chr1    1000    2000

```

lp1.bed:

```
chr1    20  50
chr1    30  40
chr1    1500    1900

```

lp2.bed:

```
chr1    68  786
chr1    899 987

```

How can I create this output:

```
chr1    0   1000    2   2
chr1    1000    2000    1   0

```

In this output columns 1, 2, 3 are the standard bed format cols, of the chromosome, start, and end of the features taken from window.bed above. Column 4 is the number of intersections between window.bed and lp1.bed per feature in window.bed. For example, the first element in column 4 is '2' because there are two intersections between window.bed and lp1.bed. Column 5 is the same but for the number of intersections between window.bed and lp2.bed for each feature in window.bed.

In the case that a feature is a subset of another feature in either lp1.bed, or lp2.bed, they should be counted as two intersections. For example, the second feature of lp1.bed is a subset of the first feature. This is reflected in the first element of column 4 of the output.

I suspect that it is doable with bedtools intersect, but I cannot work out the syntax.

# Answer
The relevant `bedtools` option seems to be `-c`:

```
    -c  For each entry in A, report the number of overlaps with B.
        - Reports 0 for A entries that have no overlap with B.
        - Overlaps restricted by -f, -F, -r, and -s.

```

We can get the number of overlaps per region per file with:

```
$ bedtools intersect -c -a window.bed -b lp1.bed
chr1    0   1000    2
chr1    1000    2000    1

$ bedtools intersect -c -a window.bed -b lp2.bed
chr1    0   1000    2
chr1    1000    2000    0

```

So we're almost there, we just need to combine them, and that's a job for `paste`:

```

$ paste  <(bedtools intersect -c -a window.bed -b lp2.bed) \
         <(bedtools intersect -c -a window.bed -b lp1.bed)
chr1    0   1000    2   chr1    0   1000    2
chr1    1000    2000    0   chr1    1000    2000    1

```

But that gives us extra fields, so we can use `cut` to remove them:

```
$ paste  <(bedtools intersect -c -a window.bed -b lp1.bed) \
         <(bedtools intersect -c -a window.bed -b lp2.bed | cut -f 4)
chr1    0   1000    2   2
chr1    1000    2000    1   0

```

The `paste` command just joins lines from its two input files in the order it receives them, but the bedtools commands will always output the same order as the lines in `window.bed` so this should be fine.

The `<(command)` construct is called "process substitution" and is a way of passing the output of a command as a file to another command that expects files as input. Note that it isn't supported by all shells, but if you're using bash or zsh which are the most common these days, it should work.
> 2 votes

# Answer
Via `bedmap --count`, which gives the option to specify a custom value when regions do not map:

```
paste <(bedmap --echo --count --unmapped-val 0 --delim '\t' window.bed lp1.bed) <(bedmap --count --unmapped-val 0 window.bed lp2.bed) > answer.bed

```

Ref. https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#element-and-overlap-statistics
> 2 votes

---
Tags: bedtools, bed
---


# Question
Title: Refactoring pandas using an iterator via chunksize
*This question was also asked on Stack Overflow*

**Bioinformatics rationale** eggNOG files can be very big and sump all available RAM for regular to medium sized desktops.

---

I am looking for advice on using a pandas iterator.

I performed a parsing operation using Python pandas, the size of the input files is eggNOG a program used in metagenomic gene function identification but has broader application. The files are large and parsing is resulting in 'RAM bottleneck' phenomenon and will not processing/parse the file.

The obvious solution is to shift to an iterator, which for pandas is the chunksize option

```
import pandas as pd
import numpy as np

df = pd.read_csv('myinfile.csv', sep="\t", chunksize=100)

```

Whats changed with the original code is the `chunksize=100` bit, forcing an iterator.

The next step is just to perform a simple operation, dropping a few columns and moving all '-' characters to np.nan then writing the whole file to disk.

```
df.drop(['score', 'evalue', 'Description', 'EC', 'PFAMs'],axis=1).replace('-', np.nan)
df.to_csv('my.csv',sep='\t',index=False)

```

How is the step performed under an iterator?

# Answer
> 0 votes
The solution posted by @el\_oso and @Corralien here was as follows:

> if you're already having memory issues reading the file, the below will actually read and write to the new file in chunks.

```
cols_to_keep = {'col1', 'col2', 'col3'}
df = pd.read_csv('myinfile.csv', chunksize=100, usecols=cols_to_keep)
with open('my.csv', 'w') as f:
    f.write('\t'.join(cols_to_keep) + '\n')
    for chunk in df:
        #do your processing here, appending that chunk to the file
        chunk.to_csv(f, header=None, sep='\t')

```

> f is \_io.TextIOWrapper i.e. a stream you can consistently write to in small memory chunks

---

This is very elegant because there are two components which are both cool and highlight why `pandas` is so widely used.

1. Instead of loading a large file and then dumping stuff thats not wanted: don't load junk at point source.
2. Simply open the write file outside the chunking loop will then stream the chunks to straight to disk and once the final chunking loop completes the file snaps shut.

The process is

* Load a chunk of data, but don't load any 'junk'.
* Process the chunk of data and write it to disk.
* Load another chunk of data repeating the above two steps.
* Continue until the entire file is processed.

The "do processing" is an OO module (there's alot going on) and the chunk is actually `self.chunk`. None of the original code needs to change its simply restructured inside the `open` loop to stream chunks to disk.

The beauty of this solution is `chunksize` is a variable, here `100` is just the example. Python's iterator called `yield` is one line at a time iteration, but this the automates a defined 'chunk' at a time. Thus the runtime speed can be tuned against the "RAM bottleneck".

**Testing** The chunked outfile was then reduced via `groupby` to remove the "duplicates" and frequencies summed (code not shown). The `chunking` OO module output file is identical to "non-chunked" OO module and written back into the module.

Its worth upvoting the SO answers.

---
Tags: python, phylogenetics, phylogeny, pandas, kegg
---


# Question
Title: Blast output file only shows 500 lines -outfmt 6
I had created databases of different sets of metagenome datasets - one with 6 runs, other with 48 runs, another with 100 runs, etc using the accession list for each of these datasets and makeblastdb command. Since I wanted to check the abundance of one specific 16S rRNA gene sequence, I ran a blastn command using this 16S sequence (F7) as the query against each of the created databases with a txt file (F7\_droughtS.txt, F7\_Dy.txt, so on) as the output. Ideally I should have got differently sized files with different number of OTU matches. I did get different files since the SRR ids were different for each database, but the results did not go beyond 500 lines for any database search. Can anyone explain why this is so?

Following is the command I used for blastn of: `sequence.fasta` against database `firstdb` :

```
blastn -query sequence.fasta -db firstdb -outfmt 6 -out output1.txt

```

Also ran the following to know the default number of lines : `blastn -help -max_target_seqs`

Number of database sequences to show one-line descriptions for

> Not applicable for outfmt \> 4 Default = 500' * Incompatible with: max\_target\_seqs -num\_alignments \<Integer, \>=0\> Number of database sequences to show alignments for Default = 250' * Incompatible with: max\_target\_seqs

# Answer
> 1 votes
I've reproduced the issue. The solution is to simply set `-max_target_seqs` to whatever number of hits you want. E.g. `-max_target_seqs 1000` gives 1000 and so on ... The problem is specific to the `-outfmt` option and not default searches.

**blastn version** I'm using v2.12.0, Feb 8 2022

---

**Proof** My database contains around 260000 fasta sequences.

```
makeblastdb -in mystuff.aln.0.fa   -parse_seqids -blastdb_version 5  -out test.db -title "test" -dbtype nucl
blastn -query query.fas -db test.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

```

Then `wc blasthits.txt`

> \>5000 hits

If we now use `-outfmt 6`, but remove the stringency of `-evalue 1e-6` then ...

```
blastn -query query.fa -db test.db -outfmt 6 -num_threads 4 -out blasthits.txt

```

Now perform `wc blasthits.txt`

> 500

So the OP is correct, I didn't know `-outfmt 6` restricts the output to 500.

If we now introduce the `max_target_seqs 1000` parameter, suggested by @Laura then,

```
blastn -query query.fa -db test.db -outfmt 6 -max_target_seqs 1000 -num_threads 4 -out blasthits.txt

```

Then `wc blasthits.txt`

> 1000

If its 550 `wc ...` is 550 etc...

Thats it. Note in my blasts I usually parallelise across 4 processes i.e. `-num_threads 4`.

---
Tags: phylogenetics, phylogeny, metagenome, blastn, makeblastdb
---


# Question
Title: Assign score to PyRanges from another PyRanges object
I have two PyRanges objects. PR1 is a set of regions of interest across a chromosome, PR2 tiles the genome and has a score associated with each range.

I need a way to get for each range in PR1 its score based on its overlap with PR2. If it is contained within just one range in PR2, it inherits the score of that range. If it overlaps multiple ranges on PR2, the score should be a weighted average of those ranges based on the size of the overlap.

Applying this to the image above, it would be:

```
score(A) = score(X)
score(B) = (OverlapBX * score(X) + OverlapBY * score(Y)) / sizeB

```

I've been browsing the PyRanges docs and I see a lot of functions that do almost that, but most of them don't take into account a score.

I had done this before in R using GRanges, where I could iterate over the overlap object to piece together scores, but I can't find a direct port of R's `findOverlaps()` in PyRanges.

---

The inputs I'm using for this are a list of mutation positions and a map of recombination rates over the chromosome they're on. So if two mutations are in a section with a constant recombination rate, the section between them inherits that recombination rate, but if the map says there's different rates between them, I need to take that into account.

Here's a pastebin with the data I'm working with: https://pastebin.com/sBYsjRQB

---

The Unfun Cat's first suggestion:

```
dist_to_pos = distances.join(positions, report_overlap=True).drop(like="_b")
df = dist_to_pos.df
g = df.groupby("Id2")
avg_score = g.cMperMb.sum() / g.Overlap.first()
avg_score.name = "AvgScore"
j_with_avg_score = pr.PyRanges(df.set_index("Id2").join(avg_score))

```

The results this give are... curious. Here's what it gives when applied to my data:

```
j_with_avg_score

Start   cMperMb     cMatstart   End     Chromosome  Id1     Overlap     AvgScore
0   2916415     1.223681    3.272236    2917484     chr1    0   1   12.236813
1   2916415     1.223681    3.272236    2917484     chr1    1   24  12.236813
2   2916415     1.223681    3.272236    2917484     chr1    2   5   12.236813
3   2916415     1.223681    3.272236    2917484     chr1    3   74  12.236813
4   2916415     1.223681    3.272236    2917484     chr1    4   110     12.236813
...     ...     ...     ...     ...     ...     ...     ...     ...
947     2960324     0.269208    3.291457    2960462     chr1    913     138     0.001951
948     2960462     0.352644    3.291494    2960790     chr1    913     77  0.018319
949     2960462     0.352644    3.291494    2960790     chr1    914     6   0.018319
950     2960462     0.352644    3.291494    2960790     chr1    915     15  0.018319
951     2960462     0.352644    3.291494    2960790     chr1    916     71  0.018319

952 rows × 8 columns
## input data contains 917 sections --> mismatch!

res_dic = {}
for x,y in zip(j_with_avg_score.Start, j_with_avg_score.AvgScore):
    res_dic[x] = y
res_dic

{2916415: 12.236813262999997,
 2917484: 0.02265139573888889,
 2917538: 0.06125764894,
 2917608: 0.10651243958260867,
 2918281: 0.6792060563708334,
 2918836: 1.1527867390640016,
 2921187: 4.176119918333334,
 2921274: 1.2543028034,
 2921300: 0.9058800047555555,
 2921976: 3.26162101122,
 2922306: 0.4522015679581395,
 2923593: 0.7409518310923077,
 2924172: 8.386433179800006,
 2925614: 0.31748233525000014,
 2926340: 0.17712928483599996,
 2926675: 1.4567026153333333,
 2926730: 5.063891720272731,
 2927172: 2.42991683355,
 2927410: 0.9370814839,
 2927411: 0.1303501307571428,
 2928125: 0.20534225654799998,
 2928557: 0.19228154158888885,
 2928849: 1.5796778519999994,
 2929269: 1.9946231915199983,
 2931597: 0.15826380113170732,
 2933043: 0.010025524,
 2934387: 0.05559167711694925,
 2938888: 0.0023231045055555556,
 2940194: 0.01882597932,
 2941104: 0.1670264778,
 2941694: 0.0197562495886076,
 2942700: 0.3931709166833331,
 2946526: 1.6788888500999994,
 2947460: 0.19597836537931007,
 2951262: 5.532616244500001,
 2956899: 0.44889279189999975,
 2960324: 0.0019507828376811595,
 2960462: 0.018319161953246753}

```

Something's getting chewed up here, due to the number of unique sections in the output, the duplicates, and the start and end positions not matching the input data.

However, the first value calculated is exactly the same as what I had gotten in my R script (with a few magnitude adjustments): `2916415: 12.236813262999997` vs `2917178: 0.00012236813263`

So the calculations happening are most likely correct, I think it's mainly the start and end positions that are getting messed up.

# Answer
I have some inspiration below. It computes a weighted score based on overlap.

I need to know what exactly you want the overlap BX and BY to be. Could you post some example data and expected output?

```
# initialize test data
import numpy as np
np.random.seed(0)

import pyranges as pr
a1 = pr.data.aorta().drop()
a1.Id1 = list(range(len(a1)))
a2 = pr.data.aorta2().drop()
a2.Id2 = list(range(len(a2)))
a1.Score1 = np.random.random(len(a1)) * 100
a2.Score2 = np.random.random(len(a2)) * 100

# compute weighted average based on overlaps
j = a1.join(a2, report_overlap=True).drop(like="_b")

df = j.df

g = df.groupby("Id2")
avg_score = g.Score2.sum() / g.Overlap.first()
avg_score.name = "AvgScore"
j_with_avg_score = pr.PyRanges(df.set_index("Id2").join(avg_score))

```
> 2 votes

---
Tags: python, genomicranges, pyranges
---


title: editing dataframe in r

text:
# Question
Title: editing dataframe in r
One of the columns in my R data frame looks like this:

```
[1] "SVLEN=1007"   "SVLEN=."      "SVLEN=."      "SVLEN=."      "SVLEN=204"    "SVLEN=-379"
 [7] "SVLEN=-492"   "SVLEN=-825"   "SVLEN=-713"   "SVLEN=."      "SVLEN=3553"   "SVLEN=432391"
[13] "SVLEN=."      "SVLEN=-52"    "SVLEN=>1615"  "SVLEN=."      "SVLEN=6899"   "SVLEN=."
[19] "SVLEN=."      "SVLEN=4468"   "SVLEN=."      "SVLEN=."      "SVLEN=."      "SVLEN=627"
[25] "SVLEN=."      "SVLEN=."      "SVLEN=."      "SVLEN=1171"   "SVLEN=."      "SVLEN=625"
[31] "SVLEN=."      "SVLEN=."      "SVLEN=-237"   "SVLEN=."      "SVLEN=."      "SVLEN=-30"
[37] "SVLEN=39"     "SVLEN=."      "SVLEN=-30"    "SVLEN=-377"   "SVLEN=120497" "SVLEN=."

```

and I would like it to look like this:

```
[1] "1007"   "."      "."      "."      "204"    "-379"
 [7] "-492"   "-825"   "-713"   "."      "3553"   "432391"
[13] "."      "-52"    ">1615"  "."      "6899"   "."
[19] "."      "4468"   "."      "."      "."      "627"
[25] "."      "."      "."      "1171"   "."      "625"
[31] "."      "."      "-237"   "."      "."      "-30"
[37] "39"     "."      "-30"    "-377"   "120497" "."

```

I could go in to the excel file I got it from and edit it but is there an easy way to do this in R?

**ANSWER**:

I found a solution to this probelm going back to the R documentation (https://www.statology.org/gsub-r/) use `gsub` to section off characters in a string. My command for this was: `gsub(".*=","\\1",df$SVLEN_INFO)`

And that worked!

# Answer
Adding an answer for now:

You can use `sub` to remove the `SVLEN=` part. You don't need `gsub` and you don't need regular expressions, let alone capture groups.

```
df$SVLEN_INFO <- sub(x = df$SVLEN_INFO, pattern = "SVLEN=", replacement = "")

```
> 2 votes

# Answer
Following gringer's advice I am also adding this solution:

Used gsub to section off characters in a string.

My command for this was: `gsub(".*=","\\1",df$SVLEN_INFO)`

Thanks for all the help!
> -1 votes

---
Tags: r, dataframe
---


title: change the order in ggplot

text:
# Question
Title: change the order in ggplot
I'm trying to make a plot that looks like this:

This is a rudimentary plot I did in excel. It's a plot of reads by length. Very simple nothing crazy to it. The important thing to note is that it's plotting the data in increasing order of length.

When I try to use ggplot to make this same plot in R it does something like this:

Where it orients the read names in alphabetical order not in ascending order of the lengths. The dataframe/table I have read into R is already listed in ascending length order. For example:

```
1   c6dd66b4-dab5-484b-92b6-c4ecdcc155b1   1570
2   ecf15d4f-15b1-4636-a97d-e0928281e311   1739
3   461f05c7-0fc6-4564-95df-6188966fe616   1790
4   e989c32c-491e-4af1-a37d-4d4f5c87a3d9   1797
5   ce94c2b3-442b-4b8c-bd59-6d7c23d5fe3d   1824
6   ed15ea16-a28d-42ca-9a07-2818d90ddd4b   1832
7   481550cd-6f87-44c7-915f-f0f3d638ffe4   1853

```

As you can see the lengths increase as you go down and the read names are in no particular order. Can someone help me read in the correct order to make this plot actually work in ggplot?

# Answer
> -2 votes
You need to enforce the order you need on the read-name column. Let's say `df` is your data frame with `read.name` and `length` as the two columns.

This is a fairly easy operation.

```
df <- df[order(df$length),]
df$read.name <- factor(df$read.name, levels = df$read.name, ordered = TRUE)
ggplot(df, ...) #Your current ggplot code goes here

```

---
Tags: r, ggplot2
---


# Question
Title: How to get a GISAID account? I registered months ago, still no reply!
Inspired by amateur variant hunters, I would like to join the Pango lineage proposal community and help contribute to variant surveillance.

However, I cannot seem to get access to GISAID, the platform that exclusively hosts a large share of SARS-CoV-2 sequences from around the world.

I applied/registered for an account back in summer 2022 (almost a year ago) using this form: https://gisaid.org/register/

But since then, nothing has happened. I have send an email and filled out the contact form on the website but never heard back.

I am very confused how GISAID can say they are a "public" database, while at the same time not giving accounts to people like me, not even responding to emails.

How does everyone else get access? Is there a secret email you need to contact? Or do you need to pay? I've read here which describes the founder of GISAID in detail.

Any help would be greatly appreciated!

# Answer
> 3 votes
I've got a GISAID account and it's subject to an approval process. It was pretty simple and is awarded quickly. There are definitely no fees needed. I personally didn't see the advantages over NCBI's really very impressive SARS-CoV-2 database system and thats what I'd recommend.

Their resource is here: https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?SeqType\_s=Nucleotide&VirusLineage\_ss=SARS-CoV-2,%20taxid:2697049&ids=MZ039713%20MZ039682%20OV522182%20

Its perfect for interrogation, it has the full spectrum of epidemiological metadata, and further details are in this post here:

How to programatically download SARS-CoV-2 fasta from NCBI/Genbank via API?

The Science journal page seemed a very good write up of GISAID BTW.

---

If you are wanting GISAID for networking, thats a bit tricky because NCBI is not going to do that. You can reapply for GISAID access, but at the time I was much more impressed by NCBI. The speed at which NCBI produced an awesome database was truly jaw-dropping and they'd implemented full pangolineages classification at an early stage. Thus I just use NCBI.

---
Tags: phylogenetics, phylogeny, public-databases, sars-cov-2, variants
---


title: h5ad file format filter

text:
# Question
Title: h5ad file format filter
I'm trying to do simple filter for single data that is stored in `h5ad` file format using this

```
import scanpy as sc

sc.pp.filter_cells(adata, min_genes=500)
sc.pp.filter_genes(adata, min_cells=10)
print(adata.n_obs, adata.n_vars)

```

when I saved the filtered data and see the size of my original data compare it, the size of the filtered data is more, but the dimension etc are reduced due to filtering.

So is it something wrong am I doing?

**UPDATE**

```
adata = sc.read_h5ad('/Single_cell/cerebellar_development/508cf892-174c-45ab-b2dc-05f54f1ee7ed/aldinger20.processed.h5ad')
adata

AnnData object with n_obs × n_vars = 69174 × 49166
    obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'sample_id', 'percent.mito', 'S.Score', 'G2M.Score', 'Phase', 'CC.Difference', 'nCount_SCT', 'nFeature_SCT', 'age', 'figure_clusters', 'sex', 'type', 'experiment', 'fig_cell_type'
    var: 'name'
    obsm: 'X_pca', 'X_umap', 'X_tsne'

sc.pp.filter_cells(adata, min_genes=500)
sc.pp.filter_genes(adata, min_cells=10)
print(adata.n_obs, adata.n_vars)

results_file = 'adata.h5ad'
adata.write(results_file)

```

**FULL CODE**

# updated version: h5ad filter, added layers to make the filtered file usable for seurat to read and make a correct seurat object

```
#https://github.com/satijalab/seurat/issues/2485

import scanpy as sc
import numpy as np
import pandas as pd
import time
from ipywidgets import interact, interactive, widgets
from IPython.display import display, HTML

def filter_and_save_h5ad(h5ad_file, min_genes, min_cells):
    # h5ad file read
    adata = sc.read_h5ad(h5ad_file)

    # cells and genes filter parameter from user
    sc.pp.filter_cells(adata, min_genes=min_genes)
    sc.pp.filter_genes(adata, min_cells=min_cells)

    # Add 'counts' layer
    adata.layers['counts'] = adata.X.copy()

    # add 'winsorization' layer
    adata.layers['winsorized'] = adata.X.copy()

    # filtered dataset save
    filtered_file = h5ad_file.replace('.h5ad', '_filtered.h5ad')
    adata.write(filtered_file)

    print(f"Filtered dataset saved to {filtered_file}")

# input widgets
h5ad_file_widget = widgets.Text(description='File_add:')
min_genes_widget = widgets.IntText(description='Min_Genes:')
min_cells_widget = widgets.IntText(description='Min_Cells:')
button = widgets.Button(description='Filter_Save')

# timer widget
time_widget = widgets.HTML(value="<p>Time passed away: 0 seconds</p>")

def on_button_clicked(click_me):
    h5ad_file = h5ad_file_widget.value.strip()
    min_genes = min_genes_widget.value
    min_cells = min_cells_widget.value

    # timer
    start_time = time.time()

    filter_and_save_h5ad(h5ad_file, min_genes, min_cells)

    # Calculate and update the elapsed time
    elapsed_time = time.time() - start_time
    time_widget.value = f"<p>Time passed away: {elapsed_time} seconds</p>"

button.on_click(on_button_clicked)

# Display the widgets and button
display(h5ad_file_widget, min_genes_widget, min_cells_widget, button, time_widget)

```

# Answer
> 1 votes
You didn't show us how the original .h5ad object was generated, how you read it in, how you wrote it, and if you did anything else in-between besides filtering. `AnnData.write()` takes a `compression` arg that might explain the difference. As also mentioned in that link, ensuring you are writing a sparse matrix can make a substantial difference to file size (the matrix in `adata.X` is generally by far the largest part of the object)

# Answer
> 1 votes
Scanpy's `h5ad` is in fact hdf5 format. This format is used for rapid read/write (i/o) large amounts of data and thats its purpose. Its fairly easy to use from scratch.

Its purpose is not to compress data. If you are concerned with the contents of `h5ad` why not just read them directly via the `hdf5` viewer. You can download it here. https://www.hdfgroup.org/downloads/hdfview/

The file will then open up as a "Explorer" (Windows) or "Finder" (OSX) GUI file system and you can open the various files and see whats in them.

My suspicion is that the original data is been written into the file. 'cause it can i/o rapidly and access precisely what it needs then file size could be a secondary consideration. Anyway you can peer in and find out by looking at the original file and the hdf5 save you've done.

An example from the website is here

Cool right?

Personally, I think HDF5 is overrated and `feather` outperforms it for speed (i/o) - but thats just for dataframes. Its not much faster than pickling in case you are interested and pickling in addition does compression. parquet is the format that does rapid i/o but specialises in compression, again mostly for dataframes.

HDF5 however does permit very pretty GUI interfaces.

---

Also if you want to load up and look at various files, that just uses the Python `hdf5` module.

---
Tags: scrnaseq, single-cell
---