title: Genome QC + Assembly Pipeline semantics

text:
# Question
Title: Genome QC + Assembly Pipeline semantics
I’m trying to create a pipeline for genome assembly. How best can I “redirect/pipe” from existing fasta files (or files in general) to other steps of the pipeline?

I was thinking of going from the SRA Download step (prefetch and fasterq-dump) to fastqc, then to trimmomatic or bbduk to remove the adapters, one more fastqc check and then to SPAdes for the assembly.

In my (amateur) mind it would be cool to be like:

“SRRXXXXXX” | prefetch | fasterq-dump | fastqc | bbduk.sh | fastqc (second) | spades.py

Any and all help is greatly appreciate! I know this stuff is pretty basic but I am so new to everything it’s painfully obvious to even me.

# Answer
> 1 votes
Maybe take a look at workflow management systems, like

Admitted, each of those adds their additional learning curve. However, especially snakemake is not so much more complicated than a standard bash workflow. Those systems do help a lot organize and orchestrate your workflows

# Answer
> 0 votes
I understand the desire to pipe everything like that but I don't think it will work, even for such a small workflow. And as your pipe grows it'll make it harder to troubleshoot/debug.

These days everyone uses workflow languages to write pipelines. It's a bit more work to get set up but worth it in the long run. I use WDL. Snakemake is very popular. BASH will of course work for small pipes but not with | for every step.

# Answer
> 0 votes
Directly piping will definitely not work, as Liam McIntyre suggests. However, if you don't feel up to learning a workflow language just to do this, it is perfectly possible to just put the various commands in a bash driver script (or whatever you're comfortable with).

Each of those tools will take parameters to define output filenames, or have systematic ways of writing FASTA files and defining inputs. You will have to look at the command line interface of each of those tools to figure it out (for example, sratools here, bbduk.sh here and spades here). an example would be as follows:

```
#!/usr/bin/bash
SRA=$1

# get reads
prefetch $SRA  # writes to a known path locally
fasterq-dump $SRA --split-3  # looks for path_to_sra_file and outputs it. assumes you are splitting.

# trim reads
bbduk.sh in1=$SRA''_1.fastq in2=$SRA''_2.fastq out1=$SRA''_clean_1.fq out2=$SRA''_clean_2.fq

# re-QC the files- this step does not feed into the other steps!!
fastqc $SRA''*.fastq $SRA''*.fq  # runs fastqc on all fastq files

# run assembly
spades.py -o my_assembly --pe1 $SRA''_clean_1.fq --pe2 $SRA''_clean_2.fq

```

and so on.

I suppose that from a formal perspective it's better to use a workflow engine as Liam McIntyre says, but with such a straightforward workflow I'd say start simple and get more elaborate only if you need to; you'd spend 10x as long just figuring out the workflow language. This is ok if it's a goal to learn it or if you already know it, but if you just want some results putting it into a script is fine.

---
Tags: fasta, genome, assembly, sratoolkit, fastqc
---


# Question
Title: How can I use pandas agg here to avoid iteration?
can anyone help me improve this pandas code?

```
import pandas as pd

df = pd.DataFrame(
    [
        [
            'chr1', 222
        ],
        [
            'chr1', 233
        ],
        [
            'chr1', 2123
        ],
        [
            'chr2', 244
        ]
    ], columns = ['chrom', 'pos']
)

df2 = pd.DataFrame(
    [
        [
            'chr1', 221, 223
        ],
        [
            'chr1', 230, 240
        ],
    ], columns = ['chrom', 'start', 'end']
)

```

Gives me 2 dfs with genomic coordinates. The first one is an exact position:

```
    chrom   pos
0   chr1    222
1   chr1    233
2   chr1    2123
3   chr2    244

```

and the second is a range:

```
    chrom   start   end
0   chr1    221 223
1   chr1    230 240

```

I need to find the count of exact coordinates that are in one of the ranges (in the same chrom)

This works but is slow:

```
c=0
for chrom, data in df.groupby('chrom'):
    tmp = df2.query(f'chrom == "{chrom}"')
    for p in data.pos:
        for s, e in zip(tmp.start, tmp.end):
            if s < p < e:
                c+=1

```

Then c = 2

I think I can use agg to do this without iteration (and hopefully faster) but I can't get it working. Can anyone show me how?

PS I am also asking this on stackoverflow.

# Answer
> 3 votes
```
# pip install pyranges or conda install -c bioconda pyranges

import pyranges as pr

g1 = pr.from_string("""Chromosome Start End
chr1   222 223
chr1   223 224
chr1   233 234
chr1   235 236
chr1   2237 238
chr1  2123 2124
chr2   244 245""")

g2 = pr.from_string("""Chromosome Start End
chr1    221  223
chr1    230  240
chr2    0 1000""")

r = g2.count_overlaps(g1)
r
# +--------------+-----------+-----------+------------------+
# | Chromosome   |     Start |       End |   NumberOverlaps |
# | (category)   |   (int32) |   (int32) |          (int64) |
# |--------------+-----------+-----------+------------------|
# | chr1         |       221 |       223 |                1 |
# | chr1         |       230 |       240 |                2 |
# | chr2         |         0 |      1000 |                1 |
# +--------------+-----------+-----------+------------------+
# Unstranded PyRanges object has 3 rows and 4 columns from 2 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

r.df
# 0       chr1    221   223               1
# 1       chr1    230   240               2
# 2       chr2      0  1000               1

```

---
Tags: python, pandas, chromosomes
---


# Question
Title: How to filter out partially similar strings from two lists and make one list?
I have following two lists. List1

```
m64041_200717_231916/100664277/0_42364 id=2440
m64041_200717_231916/100008029/0_42112 id=244006
m64041_200717_231916/100665802/0_46854 id=2440067
m64041_200717_231916/100008044/0_60161 id=24400
m64041_200717_231916/100928073/0_11804 id=2440152
m64041_200717_231916/10094492/0_44157 id=2440199
m64041_200717_231916/100008373/4102_50601 id=2440234

```

List2

```
m64041_200717_231916/100664277/0_42364
m64041_200717_231916/100665802/0_46854
m64041_200717_231916/100928073/0_11804
m64041_200717_231916/10094492/0_44157

```

These are sequence read ids from a genome file. List1 contains all the ids from the List2 but with some id numbers (id=2440). So the sequence ids are only similar to some extent. I want to compare these two files and I need to filter out List2 sequence ids from List one (with id numbers) as below. Output.

```
m64041_200717_231916/100664277/0_42364 id=2440
m64041_200717_231916/100665802/0_46854 id=2440067
m64041_200717_231916/100928073/0_11804 id=2440152
m64041_200717_231916/10094492/0_44157 id=2440199

```

I will be grateful if someone can help me with this issue.

# Answer
You really need to learn some python :)

```
with open('/path/list2.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip())

with open('/path/list1.txt', 'r') as fin:
    with open('/path/list3.txt', 'w') as fout:
        for line in fin:
            if line.strip().split(' ')[0] in to_keep:
                fout.write(line)

```
> 3 votes

# Answer
There is no need for a highlevel language in such a simple case, an `awk` oneliner does just fine:

```
$ awk -F " " 'FNR==NR {a[$1]; next} $1 in a' List2 List1
m64041_200717_231916/100664277/0_42364 id=2440
m64041_200717_231916/100665802/0_46854 id=2440067
m64041_200717_231916/100928073/0_11804 id=2440152
m64041_200717_231916/10094492/0_44157 id=2440199

```
> 2 votes

# Answer
This can be done in python easily. If you have different text files then you can load them into a python list and follow from step 2

```
# Step 1
list1 = ["m64041_200717_231916/100664277/0_42364 id=2440",
"m64041_200717_231916/100008029/0_42112 id=244006",
"m64041_200717_231916/100665802/0_46854 id=2440067",
"m64041_200717_231916/100008044/0_60161 id=24400",
"m64041_200717_231916/100928073/0_11804 id=2440152",
"m64041_200717_231916/10094492/0_44157 id=2440199",
"m64041_200717_231916/100008373/4102_50601 id=2440234"]

list2 = ["m64041_200717_231916/100664277/0_42364",
"m64041_200717_231916/100665802/0_46854",
"m64041_200717_231916/100928073/0_11804",
"m64041_200717_231916/10094492/0_44157"]

# Step 2
split_list1 = [ele.split(" ")[0] for ele in list1 ]
mapdict = {split_list1[i]:list1[i] for i in range(len(list1))}
common = [value for value in split_list1 if value in list2]
required = [mapdict[x] for x in common]
print(required)

```
> 0 votes

---
Tags: bash, filtering
---


# Question
Title: Upgrading guppy basecaller on Oxford Nanopore Gridion via PPA
I am trying to update the `guppy_basecaller` (current version installed is 3.2.10+aabd4ec) on a ONT Gridion device. I am initially trying to update it via the nanopore PPA, as described here (requires login access to the nanopore community)

After I issued `sudo apt update`, I got a GPG error, saying that the public key was not found in the system.

```
W: An error occurred during the signature verification. The repository is not updated and the previous index files will be used. GPG error: https://europe.oxfordnanoportal.com/apt xenial-stable InRelease: The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: An error occurred during the signature verification. The repository is not updated and the previous index files will be used. GPG error: https://mirror.oxfordnanoportal.com/apt xenial-stable-gridion InRelease: The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: Failed to fetch https://mirror.oxfordnanoportal.com/apt/dists/xenial-stable-gridion/InRelease  The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: Failed to fetch http://mirror.oxfordnanoportal.com/apt/dists/xenial-stable/InRelease  The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: Some index files failed to download. They have been ignored, or old ones used instead

```

I came across a possible solution shown here, so I tried to add the key to the system with the following command:

```
sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys AD9D09AF2EBBA5A8

```

The terminal was busy for too long and, after some time I got an timeout error:

```
Executing: /tmp/tmp.EWnOg4cPJk/gpg.1.sh --keyserver
keyserver.ubuntu.com
--recv-keys
AD9D09AF2EBBA5A8
gpg: requesting key 2EBBA5A8 from hkp server keyserver.ubuntu.com
gpg: keyserver timed out
gpg: keyserver receive failed: keyserver error

```

I then came across this post where one potential solution was to change the port, so I issued:

```
sudo apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys AD9D09AF2EBBA5A8

```

But then I got:

```
Executing: /tmp/tmp.sfqFuJpVb2/gpg.1.sh --keyserver
hkp://keyserver.ubuntu.com:80
--recv-keys
AD9D09AF2EBBA5A8
gpg: requesting key 2EBBA5A8 from hkp server keyserver.ubuntu.com
gpgkeys: key AD9D09AF2EBBA5A8 not found on keyserver
gpg: no valid OpenPGP data found.
gpg: Total number processed: 0
gpg: keyserver communications error: key not found
gpg: keyserver communications error: bad public key
gpg: keyserver receive failed: bad public key

```

Has anyone managed to update the `guppy_basecaller` via PPA, as described in the nanopore community page?

An alternative would be to install via the `.tar.gz` file, though I am first trying to upgrade the current version from the PPA to avoid conflicts with the current installation.

# Answer
This is an ONT problem, not an ubuntu problem. Keys can expire over time, and need to be updated.

The keys for the Guppy basecaller are not present in the ubuntu repository; you need to get them from ONT. This is indicated by the installation instructions for Guppy (which you've already linked to):

```
sudo apt-get update
sudo apt-get install wget lsb-release
export PLATFORM=$(lsb_release -cs)
###
wget -O- https://mirror.oxfordnanoportal.com/apt/ont-repo.pub | sudo apt-key add -
###
echo "deb http://mirror.oxfordnanoportal.com/apt ${PLATFORM}-stable non-free" | sudo tee /etc/apt/sources.list.d/nanoporetech.sources.list
sudo apt-get update

```
> 1 votes

---
Tags: nanopore, long-reads, base-calling
---


title: Dataset showing evolution of gene

text:
# Question
Title: Dataset showing evolution of gene
As a way to exemplify how the hierarchical clustering of particular DNA sequences can be used to derive phylogenetic trees, I'm looking for a (even very small) dataset of the same genes (or other biological sequences) from different organisms (for example animals where the evolutionary tree is known). Do you have any pointers?

Thanks in advance

# Answer
If you're building a tree of different species, it's probably better to use protein sequences. A common resource is Uniprot.

If you still want DNA sequences you can turn to NCBI Nucleotide.

Some examples that give decent trees: p53, TOP1, HBA.

A common way to construct a phylogenetic tree is to start with a multiple sequence alignment (MSA).
> 0 votes

---
Tags: gene, phylogeny, phylogenetics, clustering
---


# Question
Title: How to search for high coverage SRA entries
### The Question

I want to find high coverage SRA entries, e.g., above 100x.

I guess the best way is to use https://www.ncbi.nlm.nih.gov/sra with an appropriate search term. I don't mind if the search results contain some "false-positives". (I.e., if some of the found entries have coverage above 80x, while I aimed for 100x, that's fine.) I want to avoid using Entrez.esummary on a very large number of SRA entries (e.g., through Biopython), as it would be significantly slower (if I understand correctly).

How can I do that?

---

### My Attempt

I tried using the `Mbases` field.
Assuming the reference genome is at most 100kbp (e.g., in case I am interested in viruses), I thought a good approximation would be to search for SRA entries with more than 100\*100k bases, i.e., at least 10M bases. According to https://www.ncbi.nlm.nih.gov/books/NBK3837/, I thought I could use the search term `10:9999999[Mbases]`, but disappointingly, this didn't work. Just to make sure my syntax is right, I verified that `100:101[ReadLength]` works. So it seems that `Mbases` doesn't support the range syntax.

# Answer
> 1 votes
If processing a huge table with lots of "false positives" is OK, then you could enter your search term, from top right choose "Send to" \> "File" Format: "RunInfo".

---
Tags: ngs, database, public-databases, phylogenetics
---


# Question
Title: How to identify genes from a genome assembly of C. Elegans?
I have two full genome assemblies for C. Elegans samples collected from two different geographical areas that I found on WormBase. These are in fasta format. I want to go gene-by-gene and compare the nucleotide sequences corresponding to the gene between the two samples. My goal with this is to count how many single-base differences there are in each gene. However, I'm not sure how to match the sequences in the assembly to known reference sequences (by the way is there a good resource to download all of these at once?) for C. Elegans genes.

My thought was to run BLAST between the reference sequence and each genome and find the match with the best score, if any. Is this a reasonable approach, or is there some better way to do this?

# Answer
You don't want BLAST here, or at least not regular BLASTn. A better method would be to use tBLASTn to map the known *C. Elegans* proteins to the translated genome. However, I would recommend using a slightly more sophisticated approach:

1. Collect fasta sequences of all known *C. Elegans* proteins (amino acid sequences, not nucleotide). This should be easy enough to do on WormBase. For example, the protein sequences for *C. elegans* strain VC2010, assembly PRJEB28388, can be found here: ftp://ftp.wormbase.org/pub/wormbase/releases/current-production-release/species/c\_elegans/PRJEB28388/c\_elegans.PRJEB28388.WS279.protein.fa.gz
2. Use a tool that can model introns and exons to map a protein into a genome. I haven't worked in this field in more than 10 years now, so it is very likely that there are other programs around today, but back in the day, I would do this using either `exonerate` or `genewise`. With `exonerate`, the command would be:

```
exonerate -m protein2genome -n 1 -t assembly.fasta -q allproteins.pep > out

```

That will give you a list of best-hits for each of your input proteins and should accurately map the protein to the genome. Once you have that, you can extract the matches and start comparing.
> 2 votes

---
Tags: sequence-alignment, gene, sequence-annotation, blast, public-databases
---


# Question
Title: How can I get bedtools to tell me which genes are being expressed?
I'm trying to align the Acinetobacter baumanii genome to a genome. I've already done the alignment, and I want to use bedtools to see which genes are being expressed exactly. When I try running the following command,

```
bedtools coverage -a outCSEQSAligned.out.bam -b GCF_000746645.1_ASM74664v1_genomic.bed

```

unfortunately, it doesn't tell me which genes are being expressed at each place but instead tells me that the sequence aligns to the single Acinetobacter chromosome which isn't too helpful:

```
NZ_CP009257.1   282134  282254  GWNJ-1012:340:GW200702637th:1:1105:15519:10019/1        255     +       282134  282254  0,0,0   1       120,    0,4       120     120     1.0000000

```

I'm assuming that 282134-282254 are the positions on the chromosome that correspond to that alignment. I suppose I could extract these positions and figure out which genes they correspond to, but I'm wondering if there's an easier way to use bedtools to figure out which genes are being expressed. For reference, outCSEQSAligned.out.bam is my bamfile which I aligned the Acinetobacter genome to, and GCF\_000746645.1\_ASM74664v1\_genomic.bed is a bedfile corresponding to the Acinetobacter genome.

Thanks for the help!

# Answer
`Bedtools` is not the best tool for this. You can use `featureCounts`, which gives you the number of read that overlap with features in a genome. In this case, your features are genes.

`bedtools coverage` is used to calculate the overlap between each read in BAM and each line in BED. It is not very straightforward to quantify gene expression. Plus bedtools does not account for exon in each gene and summarize them to gene level, like `featureCounts` does.
> 0 votes

---
Tags: sequence-alignment, gene-expression, bedtools
---


title: ATAC seq density calculation

text:
# Question
Title: ATAC seq density calculation
This paper Fig2b

> b Changes of chromatin accessibility in AMD (n = 14) relative to normal (n = 11) in all retina samples. Each dot represents one ATAC-Seq peak. Blue line in the left panel indicates average fold changes of peaks with the same ATAC-Seq intensity. **The percentage of reduced peaks is shown under the density curve in the right panel.**

How is the density calculated or what is being shown there? average ATAC seq signal or Fold change?

# Answer
> 1 votes
They are not referring to density of reads. It is the density of the fold change, the y variable in the plot to left of Fig2b. It tells you the distribution of the foldchange and 92% of the fold change is \<0.

See this page for more details on density and histograms

Suppose your data is something like this (purely making this up):

```
set.seed(555)
fc = rnorm(1000,-0.5,2)
averATAC = -0.5*(fc - 0.5)^2 + 0.2*fc + 3  + rnorm(100,25,10)

```

You fit a density to it and you get the plot:

```
par(mfrow=c(1,2))
plot(averATAC,fc)
plot(dens$y,dens$x,xlab = "density",ylab="fc")

```

# Answer
> 0 votes
Just count the number of reads overlapping an ATAC-seq peak and divide it by the length of the peak gives you the read density. Each dot is a fold change, but to calculate fold change, average density from all the samples are taken, that is what the "average" is referring to.

---
Tags: atac-seq
---


# Question
Title: Python equivalent to R GRanges and IRanges
**Problem:**

I am trying to convert some codes written in R to Python and part of that conversion process is find classes equivalent to the GRanges and IRanges from the GenomicRanges R package in Python. https://bioconductor.org/packages/release/bioc/vignettes/GenomicRanges/inst/doc/GenomicRangesIntroduction.html#granges-genomic-ranges

I couldn't find any equivalent library written in Python that does the same operations as the aforementioned R classes. There is however pyranges library in Python but it is not as flexible as GRanges and doesn't have IRanges implemented in it. https://github.com/biocore-ntnu/pyranges

For example in R I could add the following in the GRanges class:

```
gr = GRanges(
seq = Rle(DF$chr),
ranges = IRanges(DF$Start,DF$End),
ref = DF$ref,
alt = DF$alt,
sampleID=DF$sampleID,
seqlengths = chromSize)

```

However, pyranges don't allow you to add arguments other than (chr, start, end, strand). However, in my case, I need to have the option to add sample\_id, seqlengths and other arguments.

Thank you in advance. Best

# Answer
> 3 votes
I am not sure if your question has anything to do with IRanges. If I get it correct, the limitation you have with pyranges now is adding extra columns to the meta data, like adding values to `@values` slot in `granges`.

You can do it using `insert` or `setattr` :

```
import pyranges
import pandas as pd
import numpy as np
gr = pyranges.random(10)

setattr(gr,'alt',np.random.choice(['A','T','G','C'],10))

gr.insert(pd.DataFrame({'ref':np.random.choice(['A','T','G','C'],10)}))

+--------------+-----------+-----------+--------------+------------+------------+
| Chromosome   | Start     | End       | Strand       | alt        | ref        |
| (category)   | (int32)   | (int32)   | (category)   | (object)   | (object)   |
|--------------+-----------+-----------+--------------+------------+------------|
| chr1         | 27431195  | 27431295  | -            | C          | G          |
| chr2         | 197045893 | 197045993 | -            | A          | G          |
| chr4         | 86316012  | 86316112  | -            | C          | A          |
| chr8         | 4560598   | 4560698   | -            | G          | G          |
| ...          | ...       | ...       | ...          | ...        | ...        |
| chr11        | 130395526 | 130395626 | -            | A          | A          |
| chr12        | 28557156  | 28557256  | +            | T          | G          |
| chr13        | 55519337  | 55519437  | -            | T          | C          |
| chr14        | 1807028   | 1807128   | +            | G          | A          |
+--------------+-----------+-----------+--------------+------------+------------+
Stranded PyRanges object has 10 rows and 6 columns from 10 chromosomes.
For printing, the PyRanges was sorted on Chromosome and Strand.

```

Or directly adding them like this (thanks to @ChrisRands for pointing it out):

```
gr.alt = np.random.choice(['A','T','G','C'],10)

```

# Answer
> 0 votes
In pyranges you would not need a call like the one you show. You can just convert a dataframe with the columns you mention above to a pyranges like this:

```
import pyranges as pr
gr = pr.PyRanges(your_df)

```

This avoids the code duplication shown in your example. The column names in the df are preserved in the PyRanges.

---
Tags: r, python, bioconductor
---


# Question
Title: How to convert BAM file to bigWig in Python?
I am manipulating some BAM files using `pysam` package which seems to be very fast and handy. However, I ran into problem when I am trying to generate bigWig files from BAM files in Python. I used `bamCoverage` from `deeptools` to do this previously, but it is a standalone tool that cannot be used from within my Python scripts.

Therefore, I am wondering what Python package would allow me to read in BAM files and convert them to bigWig files.

Thanks

# Answer
> 4 votes
deepTools has a (somewhat poorly documented) API, since it's a python package too. The basic code framework is:

```
from deeptools import writeBedGraph
from deeptools.getScaleFactor import get_scale_factor

wr = writeBedGraph.WriteBedGraph(options and input files)
func_args = {'scaleFactor': get_scale_factor(...options...)}
wr.run(writeBedGraph.scaleCoverage, func_args, output_file_name, ...)

```

# Answer
> 3 votes
pyranges:

```
import pyranges as pr
gr = pr.read_bam("your_file.bam")
gr.to_bigwig("out.bw", chromosome_sizes=pr.data.chromsizes()) # for hg19

```

---
Tags: python, sam, pysam, deeptools, bigwig
---


title: genotyping or variant calling in R

text:
# Question
Title: genotyping or variant calling in R
Are there any good options for calling variants or even just a pileup from R? R is not ideal for this, but I'd like to integrate with other functions.

I found deepSNV::bam2R which roughly does what I am looking for. Are there any alternatives?

# Answer
Samtools/bcftools/tabix is available in R via Rsamtools.

SEQprocess also provides wrappers for different tools such as GATK and bowtie.

More utility packages are noted here.
> 1 votes

---
Tags: r, sam, variant-calling, genotyping, variants
---


# Question
Title: Create new range based from existing range (python)
I have a pyranges object that I am looking to iterate over. For each range, I would like to "resize" it to its centre (so `start+end/2`) and save that into a separate pyranges object. There will be additional processing associated based on what ranges already exist nearby, so I can't just use a `resize` function, I have to iterate over the object.

First, how would I go about manually creating a new pyranges object, then inserting one by one new ranges into it?

Second, when I try to iterate over my current pyranges object (using this line `for myspot in pr.itergrs(hotspots, strand=False) :`), it gives me this error: `AttributeError: 'tuple' object has no attribute 'stranded'`. This error occurs whether or not I have strands in my pyranges, and whether or not I specify `strand=False`. Why does this occur?

---

The data I am working with: a pyranges covering hotspots in the human genome (all chr1 for the moment)

```
+--------------+-----------+-----------+-----------+
| Chromosome   | Start     | End       | Score     |
| (category)   | (int32)   | (int32)   | (int64)   |
|--------------+-----------+-----------+-----------|
| chr1         | 154421820 | 154421821 | 26413     |
| chr1         | 30749469  | 30749470  | 14833     |
| chr1         | 55259503  | 55259504  | 9900      |
| chr1         | 14985097  | 14985098  | 7735      |
| ...          | ...       | ...       | ...       |
| chr1         | 247975751 | 247976497 | 1         |
| chr1         | 247977599 | 247981391 | 1         |
| chr1         | 248018012 | 248019454 | 1         |
| chr1         | 248177298 | 248178490 | 1         |
+--------------+-----------+-----------+-----------+

```

The hotspots are sorted by their score, and I'd like to go through and for each one, resize it to its centre, keeping the score, but if there is already a hotspot within a certain range (say +/-5kb), to discard it instead, so that I end up with clear chunks between hotspots. This is the "pseudocode" I was planning to do to implement this:

```
add centre of hotspot 1 to a new pyranges object
create ranges with score = -1 for 5kb flanking each side
for myhotspot in hotspots:
  resize to centre
  if not overlap with new pyranges object
    add centre
    create flanking regions
discard ranges with score<0

```

This may not be the most elegant solution, but I'll never have more than ~2000 ranges to process so it should work fine. Is this implementable with pyranges?

# Answer
> 1 votes
Quick and dirty (not to mention atrociously bad `O(n)`) solution:

```
import pyranges as pr
import numpy as np
np.random.seed(42 * 10)

# create large df to test on
gr = pr.random(int(1e5), length=10000, chromsizes={"chr1": 249250621})
gr.Score = np.random.randint(250, size=len(gr))

def remove_worst_scores_until_no_overlap(gr):
    df = gr.df
    old_length = -1
    new_length = len(df)
    while old_length != new_length:
        df = df.drop(["Cluster", "Count"], axis=1, errors="ignore")
        df = pr.PyRanges(df).cluster(count=True).df.sort_values("Score", ascending=False)
        df = df[df.duplicated('Cluster', keep='last') | (df.Count == 1)]
        old_length = new_length
        new_length = len(df)

    return pr.PyRanges(df.drop("Cluster Count".split(), axis=1)).sort()

hotspots_no_overlaps = remove_worst_scores_until_no_overlap(gr)

```

# Answer
> 4 votes
Answer to main question:

```
import pyranges as pr
assert int(pr.__version__.split(".")[2]) >= 93, "pip install pyranges==0.0.93"
import numpy as np
np.random.seed(42 * 10)

# create large df to test on
gr = pr.random(int(1e5), length=1, chromsizes={"chr1": 249250621})
gr.Score = np.random.randint(250, size=len(gr))

result = gr.max_disjoint(slack=5000)

# small example that is easy to validate
gr2 = pr.from_string("""Chromosome Start End Score
chr1 1 2 5
chr1 4999 5000 10
chr1 10001 10002 1
chr1 15000 15001 2
chr1 20100 20101 5
""")

gr2.max_disjoint(slack=5000)

# +--------------+-----------+-----------+-----------+
# | Chromosome   |     Start |       End |     Score |
# | (category)   |   (int32) |   (int32) |   (int64) |
# |--------------+-----------+-----------+-----------|
# | chr1         |         1 |         2 |         5 |
# | chr1         |     10001 |     10002 |         1 |
# | chr1         |     20100 |     20101 |         5 |
# +--------------+-----------+-----------+-----------+
# Unstranded PyRanges object has 3 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

```

---

Could you please post an example of what data you have and what you want to end up with?

1. Inserting row after row is going to be slow, whether it is pandas or pyranges. It is best to collect all the data and then create a new dataframe or PyRanges.
2. `pr.itergrs` is for iterating over multiple pyranges at the same time. It takes a list of pyranges while you seem to have a collection of tuples. The function iterates over every chromosome/strand pair in the list and for it extracts the data from every pyranges. See the code example at the end.

To iterate over a pyranges you can do `for k, df in gr: # ...`.

Anyways, to do what you want:

```
import pyranges as pr
import numpy as np

f1 = pr.data.f1()
f1.Start = ((f1.Start + f1.End)/2).astype(np.int32) # int is 64 bit
f1.End = f1.Start + 1

```
3. I understand I have not done the addititonal processing you speak of, but if you can post an example I'm sure I'd be able to help you :)

Code example `pr.itergrs`:

```
import pyranges as pr
f1 = pr.data.f1()
f2 = pr.data.f2()
for k, dfs in pr.itergrs([f1, f2], keys=True):
    print(k)
     for df in dfs:
         print("df:")
         print(df)

('chr1', '+')
df:
  Chromosome  Start  End       Name  Score Strand
0       chr1      3    6  interval1      0      +
2       chr1      8    9  interval3      0      +
df:
  Chromosome  Start  End Name  Score Strand
0       chr1      1    2    a      0      +
('chr1', '-')
df:
  Chromosome  Start  End       Name  Score Strand
1       chr1      5    7  interval2      0      -
df:
  Chromosome  Start  End Name  Score Strand
1       chr1      6    7    b      0      -

```

---
Tags: python, pyranges
---


# Question
Title: KeyError when getting features from a genbank file with biopython with some accessions but not others
I'm very new to python but i've been using it to extract the sequence of a gene from a genbank file. The issue is is that sometimes i'll get the output i want (prints the sequence to a file) and sometimes it will return a key error. This depends on which accession i'm using. Does anyone know why it sometimes might give a key error? I thought it might be something to do with the genbank records themselves, but they look pretty similar and the gene is there (in the gene feature qualifier). EG works with HG738867.1 but not AP019703.1. Here's my code -

```
from Bio import Entrez, SeqIO

gi_genome = 'accession'
name = 'acrA'
Entrez.email = 'email'
handle = Entrez.efetch(db="nucleotide", id=gi_genome, rettype="gbwithparts", retmode="text")
record = SeqIO.read(handle, "gb")
handle.close()
element = 0
for feature in record.features:
    if feature.type == 'CDS' and name in feature.qualifiers["gene"]:
        report = 'record.features[%s]' % str(element)
        gene_sequence = feature.extract(record.seq)
        with open('output.fasta', 'a') as f:
            print('>' + gi_genome + ' ' + name, file=f)
            print(gene_sequence, file=f)
        break
    else:
        element = element + 1

```

Here's the traceback -

```
Traceback (most recent call last):
  File "/home/ubuntu/Documents/Git_Branches/Project_planning/Learning/In_progress/utils/data.py", line 11, in <module>
    if feature.type == 'CDS' and name in feature.qualifiers["gene"]:
KeyError: 'gene'

Process finished with exit code 1

```

Thanks in advance!

# Answer
> 3 votes
Good question.

feature.qualifiers is a dict, if the dict doesn't have that key it will throw a KeyError. The way your code works it expects feature.qualifiers to have 'gene' as a key every time feature.type == 'CDS'. However, it is possible to have feature.type == 'CDS' without a gene entry and in the case of AP019703.1 this comes up before you reach the gene acrA.

As simeple fix is to simply use get instead:

```
feature.qualifiers.get("gene",[])

```

Now returns a empty list if there is not key "gene", which also satisfies the 'name in' part of your if statement.

If you get stuck on stuff like this use more print statements, break up your if statements and look at the whole file (example below) to understand the structure.

```
gi_genome = 'AP019703.1'
name = 'acrA'
Entrez.email = 'email'
handle = Entrez.efetch(db="nucleotide", id=gi_genome, rettype="gbwithparts", retmode="text")
record = SeqIO.read(handle, "gb")
handle.close()
with open('AP019703.gb' ,'w') as fout:
    SeqIO.write(record, fout, 'genbank')

```

---
Tags: python, genbank
---


# Question
Title: isSpike function in SingleCellExperiment package is deprecated?
I am trying to follow a tutorial from Sanger institute (from May 2019) on analysis of single cell RNA Seq data. They use isSpike function to filter out ERCC (control) and MT (mitochondrial RNA) reads, but I am getting an error message that isSpike is a deprecated function. I am trying to use altExp instead as recommended by some people. I can't seem to figure this out. Can anyone provide some code to use instead of isSpike? I am new to bioinformatics.

I tried to swap out the isSpike function with the altExp function, but I am getting stuck on how to translate the grepl function.

# Answer
`isSpike` was provided by `SingleCellExperiment`. It is now deprecated. The notes in `SingleCellExperiment` indicate that one should use `altExps`

I am struggling with a similar problem. Can't compile a package (`powsimR`) because a dependency uses `isSpike` (not used in the code directly), so I have to hunt down where it is used by one of a very long list of dependencies.
> 2 votes

# Answer
You are on the latest version of SingleCellExperiment most likely. You should not go back to using the old version. One thing is to check out the new tutorial. In any case, the purpose of altExp() is to provide another slot of ERCC reads and MT.. etc.

So in your example (super not recommended to watch video btw) :

```
molecules = read.delim("https://raw.githubusercontent.com/jdblischak/singleCellSeq/master/data/molecules.txt",stringsAsFactors=FALSE,row.names=1)

annotation = read.delim("https://raw.githubusercontent.com/jdblischak/singleCellSeq/master/data/annotation.txt",stringsAsFactors=FALSE)

```

We can pull out the spike in and mt genes:

```
library(SingleCellExperiment)
library(EnsDb.Hsapiens.v86)

Spikein_names = grep("^ERCC-",rownames(molecules),value=TRUE)
SpikeIn = molecules[Spikein_names,]

g = genes(EnsDb.Hsapiens.v86)
MTgene_names = g[seqnames(g)=="MT"]$gene_id
MTgenes = molecules[rownames(molecules) %in% MTgene_names,]

```

Have a vector of genes to keep:

```
keep = setdiff(rownames(molecules),c(SpikeIn,MTgenes))

```

Then make the single cell object:

```
sce <- SingleCellExperiment(
assays=list(counts=as.matrix(molecules[keep,])),colData=annotation)

```

Use altExp to slot in the stuff:

```
altExp(sce, "spike-in") <- SummarizedExperiment(SpikeIn)
altExp(sce, "MTgenes") <- SummarizedExperiment(MTgenes)

```
> 1 votes

# Answer
If you already have everything in a SingleCellExperiment called `sce`, and the old code called for something like this:

```
# old syntax
isSpike(sce, "ERCC") <- grepl("^ERCC", rownames(sce))

```

You can do the following instead:

```
# new syntax
is.spike <- grepl("^ERCC", rownames(sce))
sce <- splitAltExps(sce, ifelse(is.spike, "ERCC", "gene"))

```

This will remove your ERCC counts from the regular counts assay and put them in an 'alternative experiment', calling `sce` will print `altExpNames(1): ERCC`.
> 0 votes

---
Tags: r, scrnaseq, bioconductor
---


# Question
Title: How to subset genes and its nested features from a GFF file using a gene list
I would like to subset a GFF file (gene and nested features) from a gene list.

The GFF file looks like this

```
##gff-version 3
Scaffold_1      JGI     gene    22901   45904   .       +       .       ID=Genecv11000001m.g;Name=Genecv11000001m.g
Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24192   24254   .       +       0       ID=PAC4GC:50510902.CDS.3;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510902.CDS.4;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510902.CDS.5;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510902.CDS.6;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510902.CDS.7;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510902.CDS.8;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510902.CDS.9;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510902.CDS.10;Parent=PAC4GC:50510902
Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510902.three_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510903;Name=Genecv11000002m;longest=0;Parent=Genecv11000001m.g
Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510903.five_prime_UTR.1;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510903.CDS.1;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510903.CDS.2;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24198   24254   .       +       0       ID=PAC4GC:50510903.CDS.3;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510903.CDS.4;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510903.CDS.5;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510903.CDS.6;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510903.CDS.7;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510903.CDS.8;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510903.CDS.9;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510903.CDS.10;Parent=PAC4GC:50510903
Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510903.three_prime_UTR.1;Parent=PAC4GC:50510903

```

And a have the target genes in a list such as

```
Genecv11033552m
Genecv11003131m
Genecv11036683m
Genecv11012576m
Genecv11003654m
Genecv11012587m

```

I know that is possible to subset the gff using `grep -f gene_list.txt <gff_file>`. However this extracts only the gene and mRNA features, missing the CDS and UTR entries, while I would like to subset the gene together with all its children features (mRNA, five\_prime\_UTR, CDS, three\_prime\_UTR).

This happens because the ID in CDS and UTR features are the same of the mRNA and not as in the gene feature.

Any ideas?

# Answer
> 2 votes
I think this simple script that uses pyranges solves it. Half of the below is just setup to make a minimal reproducible example.

```
# pip install pyranges
# or
# conda install -c bioconda pyranges

import pandas as pd
from io import StringIO
import pyranges as pr

#########
# setup #
#########

gff_name = "gtfo.gtf"
contents = StringIO("""Scaffold_1   JGI gene    22901   45904   .   +   .   ID=Genecv11000001m.g;Name=Genecv11000001m.g
Scaffold_1  JGI mRNA    22901   45904   .   +   .   ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
Scaffold_1  JGI five_prime_UTR  22901   23284   .   +   .   ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1  JGI CDS 23285   23423   .   +   0   ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
Scaffold_1  JGI CDS 24031   24062   .   +   2   ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
Scaffold_1  JGI gene    22901   45904   .   +   .   ID=Gene2.g;Name=Gene2.g
Scaffold_1  JGI mRNA    22901   45904   .   +   .   ID=PAC4GC:WHATEVZ;Name=Gene2;longest=1;Parent=Gene2.g
Scaffold_1  JGI five_prime_UTR  22901   23284   .   +   .   ID=PAC4GC:WHATEVZ.five_prime_UTR.1;Parent=PAC4GC:WHATEVZ
Scaffold_1  JGI CDS 23285   23423   .   +   0   ID=PAC4GC:WHATEVZ.CDS.1;Parent=PAC4GC:WHATEVZ""")

pd.read_table(contents, sep="\s+", header=None).to_csv(gff_name, sep="\t", index=False, header=False)

genes_to_keep = set(["Genecv11000001m.g"])

############
# solution #
############

df = pr.read_gff3(gff_name, as_df=True)

def all_descendants(df, genes_to_keep):
    old_len, new_len = -1, len(genes_to_keep)
    while old_len != new_len:
        genes_to_keep.update(df[df.Parent.isin(genes_to_keep)].ID.drop_duplicates())
        old_len, new_len = new_len, len(genes_to_keep)
    return df[df.ID.isin(genes_to_keep)]

desired_gff = all_descendants(df, genes_to_keep)
print(desired_gff)
#    Chromosome Source         Feature  Start    End Score Strand Frame                                ID               Name longest             Parent
# 0  Scaffold_1    JGI            gene  22900  45904     .      +     .                 Genecv11000001m.g  Genecv11000001m.g     NaN                NaN
# 1  Scaffold_1    JGI            mRNA  22900  45904     .      +     .                   PAC4GC:50510902    Genecv11000001m       1  Genecv11000001m.g
# 2  Scaffold_1    JGI  five_prime_UTR  22900  23284     .      +     .  PAC4GC:50510902.five_prime_UTR.1                NaN     NaN    PAC4GC:50510902
# 3  Scaffold_1    JGI             CDS  23284  23423     .      +     0             PAC4GC:50510902.CDS.1                NaN     NaN    PAC4GC:50510902
# 4  Scaffold_1    JGI             CDS  24030  24062     .      +     2             PAC4GC:50510902.CDS.2                NaN     NaN    PAC4GC:50510902

```

# Answer
> 1 votes
Given the nested structure of a GFF file, it may be easier to make a python script using a GFF file parser, like gff3, gffutils, or BioPython. Dealing with nested data in `bash` or with some fancy `awk` script is probably going to trickier to write than it's worth.

You can start by ensuring your list of target genes is sorted in the same order as the GFF file, then iterate through the file. When you encounter a like with `mRNA` in the 3rd column, store that ID value and extract all the subsequent lines with that same ID.

If the line does not have that ID, you know you've collected all the genes, mRNAs, etc associated with that gene and can pop it from your queue.

Repeat until you reach the end of the GFF file or your list of target genes is empty.

---
Tags: python, bash
---


# Question
Title: Difference in Calculation through GEO2R and difference in Article
Here is microarray analysis where author uses GSE50632 series on an Agilent microarray GPL17660 platform with 4 eligible samples which I guess are GSM1224991, GSM1224992, GSM1224993 and GSM1224994 among which two are healthy and two are patient samples.

However I found the difference in the calculations-

Through GEO2R analysis table it is-

```
ID  adj.P.Val   P.Value t   B   logFC   miRNA_ID_LIST   SPOT_ID
148000  1.86e-02    5.82e-04    3.47    -0.402958   3.91    hsa-miR-3195

```

However in the paper the calculated values are far different specially logFC value.

What could be the reason for this?

# Answer
Here are two example of calculations of logFC validated in GEO2R analysis which depends if values are already log2 transformed. During GEOQuery analysis there is calculation to transform where expression set is checked for value distribution and log2 transformation is applied-

```
ex <- exprs(gset)
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
          (qx[6]-qx[1] > 50 && qx[2] > 0)
if (LogC) { ex[which(ex <= 0)] <- NaN
  exprs(gset) <- log2(ex) }

```

**Example 1:**

=\> Sample Values are already log2 Transformed

```
n(TumorSamples) = 45
n(ControlSamples) = 46

a = average(TumorSamples) = 9.199053111
b = average(ControlSamples) = 4.988139348

logFC = |a - b| = 4.210914

```

**Example 2:** =\> Sample Values are not log2 transformed

```
n(TumorSamples) = 2
n(ControlSamples) = 2

a = average(log2(TumorSample1) + log2(TumorSample2) ... ) = (log2(84) + log2(21.8)) / 2 = 5.419287

b = average(log2(ControlSample1) + log2(ControlSample2) ... ) = (log2(624) + log2(665)) / 2 = 9.331306

logFC = |a - b| = 3.91202

```
> 1 votes

---
Tags: geo
---


# Question
Title: Quantify low complexity regions in DNA sequences
I have fasta files with multiple sequences. They are reads that mapped to the genome outside of probe-targeted regions. From a quick perusal, they appear to be repetitive and have low complexity. Is there a way I can quantify their complexity and use that as an explanation as to why those regions were mapped to?

# Answer
> 1 votes
"Complexity" is also known as entropy. A quick and dirty way is to use a file compression utility - zip will get good compression for a repetitive sequence.

You can also look at the k-mer distribution for some suitable k: in a repetitive sequence it will be biased. Here is a tool that give a local plot: https://rdrr.io/github/vsbuffalo/qrqc/man/kmerEntropyPlot.html

# Answer
> 1 votes
Use Picard EstimateLibraryComplexity.

From the docs:

> Estimates the numbers of unique molecules in a sequencing library.

```
java -jar picard.jar EstimateLibraryComplexity \
     I=input.bam \
     O=est_lib_complex_metrics.txt

```

# Answer
> 0 votes
How long are your reads, and how many? Some colleagues of mine wrote https://github.com/eclarke/komplexity to quickly score complexity in large numbers of short reads. The score is derived from the number of unique k-mers normalized by the sequence length, so it's simple and fast, though the score inherently becomes less informative for longer and longer sequences (so, best for short reads). We've used it for filtering repetitive or other low-complexity sequences in metagenomics projects.

---
Tags: ngs, read-mapping
---


# Question
Title: Oxford nanopore promethion; does guppy3 separate fast5 into pass and fail? If not what does?
just a quick question about fast5 pass fail folders. I have been working under the assumption that the fast5 pass fail folders I was given as raw data from our ONT vendor came from them doing base calling with guppy 3.something. Is that right? (The person that worked for the vendor has since moved on so I can't ask them - so my question is really does guppy3 divide fast5s into pass/fail?)

I further assumed that given I am now doing my own basecalling with guppy4.4 that I should use all the fast5s (both pass and fail) as some of the fails might become pass with the improved NN. Again, am I right? If not could someone please explain?

I'm going to also put this question on the ONT community site, FYI

Thanks,

Liam

# Answer
> 1 votes
Only if they've enabled the `--qscore_filtering` flag when using Guppy; this flag isn't enabled by default when running Guppy from the command line.

It's a good idea to use all fast5 reads for recalling, because accuracy improves over time, so some that previously failed may slide into the pass folder.

FYI, I'd be much happier if people asked their nanopore sequencing questions here. It's a publicly-accessible website which is indexed by search engines, so makes it easier to find answers, and easier to see the answers when found (i.e. no login needed).

Note that when recalling, you set the input and output directories (I think via the `-i` and `-s` flags). These should be *different* directories, and ideally different output directories for each recalling attempt, to avoid duplicate file entries. The original fast5 files shouldn't be modified by recalling as long as the input and output directories are different.

---
Tags: nanopore, long-reads, qc, fast5
---


# Question
Title: How to import a phylogenetic tree (in nexus format) with node support values in R
I have a phylogenetic tree in nexus format that was generated with RAxML. In the tree file, there are node support values. I want to plot this tree in R with the support values showing. However, when I import the tree file in R with the read.nexus function in ape, the node values do not appear in $node.labels. Can someone recommend another option in R?

# Answer
> 1 votes
I recommend you this tutorial about ggtree. You will find the solution. \[https://guangchuangyu.github.io/ggtree-book/short-introduction-to-r.html\]\[ggtree\]

---
Tags: r, phylogeny, ggtree
---


# Question
Title: How is the t-statistic value calculated in GEO2R or Bioconductor?
I was trying to calculate t-statistics using Python's `scipy` and `numpy as np`-

```
scipy.stats.ttest_ind_from_stats(
    NP.mean(x[['GSM1224991', 'GSM1224992']]),
    NP.std(x[['GSM1224991', 'GSM1224992']]),
    len(['GSM1224991', 'GSM1224992']),
    NP.mean(x[['GSM1224993', 'GSM1224994']]),
    NP.std(x[['GSM1224993', 'GSM1224994']]),
    len(['GSM1224993', 'GSM1224994']),
)

```

For example in the dataset given below, `GSM1224991 & GSM1224992` are in first group and `GSM1224993 & GSM1224994` form second group.

```
ID_REF,GSM1224991,GSM1224992,GSM1224993,GSM1224994
148000,6.436150368369428,6.499787040655854,4.430816798843313,3.0819099697950434

```

Here is the calculation given in GEO2R analysis for the same reference but while I calculate t test it is `0.029630766122806237` instead of `-0.33768` something as shown below.

```
ID  adj.P.Val   P.Value t   B   logFC
148000  1   0.001426    -9.44   -0.33768    -2.513

```

How to calculate t-stat here? Or considering single row is wrong approach?

# Answer
> 1 votes
That is not trivial. GEO2R uses the Bayesian linear model-based framework `limma` https://bioconductor.org/packages/release/bioc/html/limma.html for the analysis. I doubt you can easily and meaningfully implement the relevant code in Python without a major effort and notable background knowledge. The user guide https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf has background information plus links to the relevant papers. Don't try to analyze things by hand, use expert software. If it is not significant, well, then it is not significant. A 2vs2 comparison is strongly underpowered, therefore lack of DEGs are not unexpected, especially for microarrays.

---
Tags: statistics, limma, geoquery
---


title: Find species from FASTA files

text:
# Question
Title: Find species from FASTA files
In a school project we got the task of finding the species from two FASTA files. We got some hints of what to look for, like it is a unicellular eukaryote and that we can use blast. We don't know if the two files are the same species.

So I wonder if there is a standard way of finding species from FASTA files using the terminal?

I have a background in informatics, so I'm not well known with bioinformatics tools.

# Answer
> 2 votes
Are you not allowed to use the web BLAST tool? That's what I use if I need to quickly find out the likely origin of a DNA sequence.

There are command-line ways to do the same thing, but I only use them if I need to search hundreds or thousands of sequences at once.

---
Tags: fasta, blast, taxonomy, homework
---


# Question
Title: What does genotyping data and beta values for SNPs mean for GWAS study?
I have a list of Genes with their alleles (RS IDs) associated with a specific condition. I have to find their genotyping data and the beta values from published studies. I looked through several studies and none of them had anything mentioned about the said terms?

Apparently its for calculating the polygenic score of the listed SNPs. Does that mean the genotyping data is the allelic possibilities and beta score as effect size? Effect size is the population size (n)?

# Answer
> 1 votes
The beta statistic is a common output field produced by programs that carry out genome-wide associations. There's a good explanation here.

I interpret a beta statistic as the degree to which a change in genotype from the reference allele alters the phenotype of interest. The association (and beta statistic) can be negative or positive. It should always be tied to a particular variant, usually called the "risk allele".

The beta statistic is related to the number of individuals found in each group with a particular genotype. In that sense it could be argued to be related to the population size (at least, assuming pre-grouping of populations by genotype), but that's not a typical interpretation.

The GWAS catalogue is a good place to start off looking for data associated with GWAS:

https://www.ebi.ac.uk/gwas/home

However, I think that's for summary statistics. If someone has asked you for raw genotype data, that is trickier to get, and tends to depend on the condition. Genotype data is personally-identifying information, and there are lots of issues associated with a blanket public release of such information, which is why it's so hard to find.

Without more information on the particular phenotype, I can't give much guidance, but useful places to look for most traits are dbGaP and UK Biobank. Both of those data sources need forms filled out and approved (to demonstrate an ethically-responsible research aim) in order to access genotype-level data.

There are a few public datasets around that include data uploaded by individuals (who have explicitly consented to their genotype data being public), but it tends to be trait-specific.

---
Tags: gwas
---


# Question
Title: lower mapping rates in salmon v0.13 compared to previous versions
Hi there :) Thanks for the tool!

I recently updated to the new salmon (from 0.8... its been a couple years) and I noticed that my mapping percentages change dramatically between the two versions.

For example, using the default settings in v0.8, I see a mapping rate of 96.25% however upon upgrading to v0.13, these mapping rates go down to 82-84% depending on which of the parameters I play with (I've permutated the validateMappings, incompatPrior, consensusSlack, minScoreFraction, ma, etc). Below I have a snippet of one of these tests. I'm currently assembling the unmapped reads to see if I can determine why these reads were not mapped.

What do you think is the most likely cause for this dramatic change in performance?

Perhaps it could be the transcript index? (I used a kmer of 31 and my reads are 100 bp).

**UPDATE:**

I have run the indexing and quasi-mapping using salmon v0.8, v0.10 (when validateMappings was released), v0.12, v0.13 and had mapping rates of 96, 90, 90 and 83 respectively. So it looks like a gradual decay of mapping rates. Or alternatively, a gradual reduction in spurious mapping?

Seems like what ever improvements from 0.12 to 0.13 changed the mapping rates significantly. I will go ahead with both mappings (v0.12 and v0.13) to see if there are differences in the DESEQ2 outputs.

Thanks a lot! (while I have presented values here for one set of fastq files, the same pattern is seen throughout all 72 of my samples from different experiments)

```
for F in $FILES ; do
  R1=${F}_L001_R1_001.fastq.gz
  R2=${F}_L001_R2_001.fastq.gz
  salmon quant \
      -p 35 \
      -i $ref_dir/$transcript_index \
      -l A \
      -1 $R1 \
      -2 $R2 \
      -o $(basename $F).VM.rFB4.IP9.quant \
      --incompatPrior 9.9999999999999995e-21 \  # default of v0.8
      --writeUnmappedNames \
      --validateMappings \
      --rangeFactorizationBins 4

done

```

---

*Update* `salmon` is now released with a `--allowDovetail` option in 0.13.1 and inspected via a `-z` dump.

# Answer
> 2 votes
*\[Most likely answer based on comments attached to the question\]*

Hi Courtney, We've looked into the data and identified the source of the different mapping rates. Specifically, the cause is discordant reads; in this case reads that are dovetailing with respect to their positions on the reference. If you run Bowtie2 with "RNA-seq" flags (i.e. disallowing discordant mappings), then you get essentially the same mapping rate as salmon here. Is allowing dovetailing reads something that should be very important in your analysis? If so, we can consider what would be necessary to add back support, but generally, they provide evidence of tenuous quality.

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


# Question
Title: P-value correction when evaluating correlation between gene and miRNA expression
First of all I apologize without the question is very basic, I am taking my first steps in bioinformatics.

## Data information

We are evaluating the correlation (using the Pearson, Kendall or Spearman method) between gene expression and miRNA expression using the corAndPvalue function of WCGNA.

The resulting structure would be a DataFrame containing all combinations between each gene with each miRNAs, containing the following columns:

```
Gene     miRNA      Correlation P-value
Gen_1    miRNA_1    0,959       0.00311
Gen_1    miRNA_2    -0,039      0.1041
Gen_1    miRNA_3    -0,344      0.0021
Gen_2    miRNA_1    0,1333      0.00451
Gen_2    miRNA_2    0,877       0.07311
...

```

## Question

Considering the huge number of correlation tests we are going to evaluate, we need to adjust the p-values to avoid correlations due to chance. Bonferroni does not seem to be the best solution, so we would use Benjamini-Hochberg method (BH). The question is:

The BH correction for the `Gen_1` x `miRNA_1` combination, should consider the p-values of all combinations that include Gen\_1 (Option 1), or should consider all the p-values of all the genes x miRNA combinations (Option 2)?

For example, let's assume an expression dataset of 20,000 genes and another of 15,000 miRNAs

**Option 1:**

To adjust `Gen_1` x `miRNA_1` we would use 15,000 p-values (`Gen_1` x `miRNA_1`, `Gen_1` x `miRNA_2`, ..., `Gen_1` x `miRNA_15000`).

**Option 2:**

To adjust `Gen_1` x `miRNA_1` we would use 300,000,000 p-values (`Gen_1` x `miRNA_1`, `Gen_1` x `miRNA_2`, ..., `Gen_1` x `miRNA_15000`, `Gen_2` x `miRNA_1`, `Gen_2` x `miRNA_2`, ..., `Gen_2` x `miRNA_15000` and so on).

## Suplementary question

Documentation of the method fdrcorrection from Python Statsmodels library suggests that for negative correlations (that could be frequent in a mRNA x miRNA correlation analysis) Benjamini-Yekutieli would work better; is that right? Or Benjamini-Hochberg method would be appropiated for this case?

Any kind of help would be much appreciated, thanks in advance!

# Answer
> 0 votes
I made the same question in CrossValidated forum and got an excellent answer!

The important part:

> You need to correct for all of the comparisons you are doing. So if that's 300,000,000 comparisons you need to correct for that many multiple comparisons.

For more information check the answer in the link above

---
Tags: statistics, gene-expression, correlation, p-values
---


# Question
Title: Extracting the CIGAR string from a BAM via Python?
Is there a standard method in Python to extract a CIGAR string from the BAM?

There are great libraries which parse the CIGAR, e.g. https://pypi.python.org/pypi/cigar/0.1

```
>>> c = Cigar('10M20S10M')
>>> c.mask_left(10).cigar
'10S20S10M'
>>> c.mask_left(9).cigar
'9S1M20S10M'
>>> Cigar('10S').mask_left(10).cigar
'10S'
>>> Cigar('10H').mask_left(10).cigar
'10H'
>>> Cigar('10H').mask_left(11).cigar
'10H'
>>> Cigar('10H').mask_left(9).cigar
'10H'

```

It looks like pysam already parses the cigar string, if I'm not mistaken:

```
import pysam

bam = 'myfile.bam'
bamfile = pysam.AlignmentFile(bam, 'rb')

for read in bamfile:
    if not read.is_unmapped:
        cigar = read.cigar
        print(cigar)

```

This outputs lists of tuples of the parse output:

```
[(5, 61), (0, 30), (5, 198)]
[(4, 11), (0, 30), (4, 248)]
[(4, 11), (0, 30), (4, 248)]
....

```

How could I simply output the CIGAR from `read`? I would prefer to create a column in a pandas DataFrame

# Answer
If you really do just want the cigar string then it's `read.cigarstring`. However, I'm not sure what you're trying to gain with the cigar package from Brent. Unless you want to get the string with the masking changed then the tuple you get from pysam is the same as what you get from cigar (with the exception of the numeric operations instead of character operations).
> 4 votes

# Answer
```
# pip install pyranges
# pip install bamread

# wget http://hgdownload-test.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeRikenCage/releaseLatest/wgEncodeRikenCageGm12878CellPapAlnRep1.bam # (half a gig)
# samtools index wgEncodeRikenCageGm12878CellPapAlnRep1.bam

f = "wgEncodeRikenCageGm12878CellPapAlnRep1.bam"

import pyranges as pr
df = pr.read_bam(f, sparse=False, as_df=True)

df[["Chromosome", "Start", "End", "Strand", "Cigar"]].drop_duplicates("Cigar")
#          Chromosome  Start    End Strand             Cigar
# 0              chr1  10579  10606      -               27M
# 36             chr1  29334  29360      +           8M1I18M
# 37             chr1  29334  29360      +           7M1I19M
# 38             chr1  29334  29360      +           5M1I21M
# 44             chr1  29339  29365      +           4M1I22M
# ...             ...    ...    ...    ...               ...
# 19655298       chrM  14550  14575      +       16M5I3M3D3M
# 19660213       chrM  14747  14775      -       15M5D5M4I3M
# 19663987       chrM  15374  15399      -   2M1I2M2D4M3I15M
# 19676864       chrM  16234  16271      +  15M1I4M12D5M1I1M
# 19677003       chrM  16264  16291      +   14M1D3M2D1M3I6M

```
> 0 votes

---
Tags: python, sam, pysam, cigar
---


title: Brackets in .Rmd code blocks

text:
# Question
Title: Brackets in .Rmd code blocks
I am making some R markdown documents for my stats class. I've noticed I get an error in some functions when I have a code block with {r} instead of just r. Like so...

``` r
this code will be fine
```

``` {r}
this code will have an error
```

When should I be using the brackets for codeblocks in .Rmd and when should I leave them out?

# Answer
Does it work without a space between the ````` and the `{r}`? It may be that the Rmd processor is getting confused about the spaces, but ignoring the bare `r`, and assuming the code is R by default.

The examples here don't have spaces:

https://rmarkdown.rstudio.com/lesson-3.html
> 2 votes

---
Tags: r
---


# Question
Title: Calling GEOquerry error while loading the library
*Note: this question has also been asked on Biostars*

I started getting a strange error when trying to load GEOquerry with `library(GEOquerry)`.

```
Error: package or namespace load failed for ‘GEOquery’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]): namespace ‘rlang’ 0.4.9 is already loaded, but >= 0.4.10 is required

```

When I run `> devtools::install_github('r-lib/rlang')` I get ...

```
Error in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]) :
  namespace ‘rlang’ 0.4.9 is already loaded, but >= 0.4.10 is required

```

Has anyone dealt with this before? Restarting my R session has no effect.

# Answer
> 1 votes
rlang 0.4.10 is the current version at CRAN: https://cran.r-project.org/web/packages/rlang/index.html

Try to install it via `install.packages`, not via GitHub directly.

---
Tags: r, bioconductor, geoquery
---


# Question
Title: Finding SNP/Indel Markers of Interspecific Variation
This is a revised version of my original question, edited to include the portions already answered by @Maximilian Press.

# Goal:

## Find SNP and indel alleles which are *fixed between* but *lack polymorphism within* two closely related species (I'm looking for the BLUE region of the Venn diagram)

# Desired Outputs:

1. All SNPs/Indels of Species A which are *exclusive to* and *lack variation within* Species A

---

2. All SNPs/Indels of Species B which are *exclusive to* and *lack variation within* Species B

---

# Scratch Materials:

1. Species A reference genome

---

2. Species B reference genome

---

3. Illumina short reads for Species A

---

4. Illumina short reads for Species B

---

My current plan is to...

1. Map Species B to Species A:

```
minimap2 -ax asm5 referenceA.fasta referenceB.fasta > BtoA.sam

```
2. Find all SNP/Indel alleles of Species B relative to Species A (Represents the *entire* circle in the above Venn diagram):

```
bcftools mpileup -f referenceA.fasta -s BtoA.sam | bcftools call -o B_compared_to_A.vcf

```
3. Find all SNPs/Indel of Species B reads relative to Species B reference (Represents the *red* region of the Venn diagram):

```
bcftools mpileup -f referenceB.fasta -b bam_listB | bcftools call -o B_compared_to_B.vcf

```
4. Perform a liftover on "B\_compared\_to\_B.vcf" to correct for coordinate differences between Species A and Species B references (I haven't completely figured out how to do this - I want to make sure I have the general pipeline nailed down before I get into the weeds)
5. Write a program which extracts all loci found in "B\_compared\_to\_A.vcf" but NOT in "B\_compared\_to\_B\_liftedOver.vcf" (This should yield the BLUE region of the Venn Diagram
6. Repeat reciprocally to get alleles exclusive to Species A

Does the logic for this new pipeline check out?

# Answer
Ok, with invaluable aid from @Maximillian Press, I've figured out how to find alleles which are variable between two species but fixed within one species.

Here is how you find SNP/Indel loci which are distinct between species A (hereafter sA) and species B (hereafter sB) but are fixed WITHIN sA:

```
# Maximillian Press's Contribution: Map sA reference to sB reference using minimap2, then convert to bam.

minimap2 -a sB_ref.fasta sA_ref.fasta > AtoB.sam
samtools view -bh AtoB.sam > AtoB.bam

# Sort bam by coordinates, then use bcftools to generate a pileup and call all alleles of sA relative to sB

java -jar $EBROOTPICARD/picard.jar SortSam I=AtoB.bam O=sorted_AtoB.bam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$TMPDIR
bcftools mpileup -O b -o AtoB.bcf -f sB_ref.fasta sorted_AtoB.bam
bcftools call -m -v -o AtoB.vcf AtoB.bcf

 # Call all alleles of sA relative to sA.
# This represents the SNP/Indel loci which are variable WITHIN sA and thus ought to be discarded from the final output.
# "bam_listA.txt" contains names of all the coordinate-sorted bam files for species A
bcftools mpileup -O b -o AtoA.bcf -f sA_ref.fasta -b bam_listA.txt
bcftools call -m -v -o AtoA.vcf AtoA.bcf

```

I then use the "sed" function to replace chromosome names with their analogs (happy to give more details if anyone else needs to do this) and also removed the VCF headers prior to making a C++ program which extracts all the rows found in AtoB.vcf but not in AtoA.vcf. My final output was thus all SNP/Indel loci which vary between sA and sB but NOT within sA. Happy to share the C++ program if anyone else needs to do this.
> 1 votes

# Answer
## Problem statement

You want to find a set of SNPs in populations of two species, for each of which you have a distinct reference genome. You then want to know which of those SNPs are *shared* between the two genomes (there is variation in each species), and which ones are *private* to only one of the species (in the other genome there is no variant at that location).

To accomplish this task, you need to first find SNP variation in each genome, and then figure out which SNPs are the same between the two genomes. In other words, line up two VCF files for different genomes and figure out which SNPs are the same between them. In comparing genome assemblies, this is sometimes called a *liftover* operation.

## Genomes can have different coordinate systems.

Let's first think about what exactly you are looking at with 2 different reference genomes. For example, here are outputs of `samtools faidx` for 2 different S. cerevisiae strains, S288C and Kyokai-7:

```
# S288C
ref|NC_001133|  230203  95  60  61
CHR2.19980913   813142  234173  60  61
CHR3.19980521   315339  1060907 60  61
CHR4.19990210   1531929 1381540 60  61
CHR5.19970727   576870  2939039 60  61
CHR6.19970727   270148  3525562 60  61
CHR7.19970703   1090936 3800252 60  61
CHR8.19970727   562638  4909411 60  61
ref|NC_001141|  439885  5481523 60  61
CHR10.19970727  745440  5928778 60  61
CHR11.19970727  666448  6686681 60  61
CHR12.19970730  1078172 7364277 60  61
ref|NC_001145|  924430  8460517 60  61
ref|NC_001146|  784328  9400452 60  61
CHR15.19970811  1091283 10197892    60  61
ref|NC_001148|  948061  11307461    60  61
ref|NC_001224|  85779   12271444    60  61

# Kyokai first 10 lines
gi|347729946|dbj|BABQ01000001.1|    5491    137 80  81
gi|347729945|dbj|BABQ01000002.1|    4704    5834    80  81
gi|347729944|dbj|BABQ01000003.1|    203046  10734   80  81
gi|347729943|dbj|BABQ01000004.1|    1036    216456  80  81
gi|347729942|dbj|BABQ01000005.1|    126714  217642  80  81
gi|347729941|dbj|BABQ01000006.1|    121900  346077  80  81
gi|347729940|dbj|BABQ01000007.1|    20385   469638  80  81
gi|347729939|dbj|BABQ01000008.1|    243458  490415  80  81
gi|347729938|dbj|BABQ01000009.1|    189876  737054  80  81
gi|347729937|dbj|BABQ01000010.1|    81973   929441  80  81
# ...
# goes on for many thousands of lines

```

What you will notice is that there are significant differences between the 2 assemblies. The first is almost all in chromosomes, the second is largely *mapped* to various chromosomes but exists in thousands of tiny contigs. There is no way to just map reads to them, call variants, and assume that you (or bcftools) can figure out what is going on.

This is possibly an extreme example, but even in the case of two genomes assembled to chromosome pseudomolecules, there is basically no chance that those two chromosomes are on the same coordinate system or have similar metadata. That means again that bcftools probably won't be able to figure it out, because a SNP at position 112929 in genome A will be homologous to position 110201 in genome B.

## Liftover: a way to solve the mapping problem between 2 genomes.

You can create a mapping between the two genomes. This is 100% a solved problem in bioinformatics, and it is more or less the same as doing liftovers between different assembly versions of the same genome. I think that you can more or less just use the UCSC LiftOver tool to convert SNP coordinates, unless I am missing something. That how-to guide that I linked is going to give you much better instructions than I can give you, here is a more minimal version that might be simpler.

The only extra suggestion that I have is that the guide uses BLAT to align the two assemblies, and minimap2 is a *much* faster tool. You may need to fight the input formatting a little bit, but if you run into runtime issues you might make that change. This is why I suggested it in the comment, though I wasn't aware that the liftover docs were so complete.

## Short pipeline:

1. Align reads (individuals of species A to reference A, individuals of species B to reference B).
2. call variants--\>get VCF files using preferred workflow for each of A and B independently.
3. compute liftover mapping (chain file, usually) from reference A to reference B.
4. use computed liftover mapping to actually lift VCF over from one genome to another.
5. align A and B references, call SNPs between them (possibly can be done from (3) output). Ensure that you are using as reference the liftover target genome so that the VCF is informative.
6. For each SNP private to each species (output of 5), query the BAMs to ensure that no variation exists.
> 2 votes

---
Tags: vcf, genomics, variant-calling, software-recommendation
---


# Question
Title: Print specific columns in a matrix on the basis of sample id's in the header
I have a matrix file (expression.txt) in which the first column is a gene\_id and from the second column the sample id's start. This matrix has 20,000 columns with sample id's corresponding to different tissue. I want to extract the samples of one tissue. Example: SM-5GZZ7 SM-5GIEN SM-5EGGH are samples from brain tissue. So I want to separate the list of samples \[(file=brain.txt;n=5000\] corresponding to brain tissue from this matrix (with bash or R)

```
Name    SM-5GZZ7    SM-5EGHI    SM-5EGHJ    SM-5N9CS    SM-5GIEN    SM-5EGHH    SM-5EGGH
ENSG00000223972.5   0   0   0   0   0   0   0
ENSG00000227232.5   187 109 143 251 113 139 199
ENSG00000278267.1   0   0   1   0   0   0   0
ENSG00000243485.5   1   0   0   1   0   0   0
ENSG00000237613.2   0   0   0   0   0   0   0
ENSG00000268020.3   0   1   0   0   0   1   1

```

The file brain.txt looks like (Number of lines in this files = 5000)

```
SM-5GZZ7
SM-5GIEN
SM-5EGGH
SM-5GZXC
SM-5GZXO
SM-5GZY1
SM-5EGHK
SM-5EGHL

```

# Answer
> 1 votes
Example in R. If you have a file of column names, `names.txt`, such as:

```
mpg
cyl
disp

```

and a data.frame like `mtcars`

```
> mtcars
                   mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360        14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D         24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230          22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
 [ reached 'max' / getOption("max.print") -- omitted 23 rows ]

```

Then you can select the columns by:

```
library(dplyr)
list_names = read.table("names.txt", header=F) %>% dplyr::select(V1)
> mtcars %>% dplyr::select(list_names$V1)
                     mpg cyl  disp
Mazda RX4           21.0   6 160.0
Mazda RX4 Wag       21.0   6 160.0
Datsun 710          22.8   4 108.0
Hornet 4 Drive      21.4   6 258.0
Hornet Sportabout   18.7   8 360.0
Valiant             18.1   6 225.0
Duster 360          14.3   8 360.0
Merc 240D           24.4   4 146.7

```

# Answer
> 1 votes
There's a million different ways to do simple taks like this. Pick one and learn it:) Python is my go to as it can do pretty much everything you'd ever want. I saved your example as csv with no spaces called expression.txt.

```
import pandas as pd
df = pd.read_csv('expression.txt')
df[['SM-5GZZ7','SM-5GIEN','SM-5EGGH']]

SM-5GZZ7    SM-5GIEN    SM-5EGGH
0   0   0.0 0.0
1   187 109 143 251 113 139 199 NaN NaN
2   0   0.0 0.0
3   1   0.0 0.0
4   0   0.0 0.0
5   0   0.0 1.0

```

---
Tags: gene-expression, bash, data-retrieval, data-preprocessing
---


# Question
Title: How to extract all peptide sequences for a UniProt protein entry with an API?
I know how to generally pull down information for a UniProt entry using the REST API, for example:

https://www.uniprot.org/uniprot/?query=id:Q9UJL9&sort=score&columns=id,go(molecular%20function)&format=tab

I'm looking at this entry as an example. I'm able to pull down most of the attributes I want with a variation of the API query above, however, I'm stuck on one part.

For all human proteins, I want to pull down the peptide sequences that come from that protein.

In the example above, it would be the sequences of the peptides associated with this section:

Can this be done with the API or do I have to come with a scraping method for all human proteins where I (1) get the protein sequence of each entry, (2) get this table, (3) calculate where the peptides start and end from the sequence in step 1? I feel like this alternative will be much slower, so I want to check I'm not missing anything obvious.

# Answer
If downloading the complete database is not an option, you may able to use the instructions from this help page (based on the gff format for any given query):

https://www.uniprot.org/help/download\_subsequences
> 1 votes

---
Tags: uniprot, api
---


# Question
Title: How to programmatically classify a protein according to its genbank feature
Say I found an interesting protein in a genbank file, e.g.:

```
 gene            273792..275552
                 /locus_tag="ECs_0242"
 CDS             273792..275552
                 /locus_tag="ECs_0242"
                 /codon_start=1
                 /transl_table=11
                 /product="type IV secretion protein Rhs"
                 /protein_id="NP_308269.1"
                 /translation="[...]"

```

(from https://www.ncbi.nlm.nih.gov/nuccore/NC\_002695.2)

I would like to programmatically get a classification of this protein.
The (not new) motivation is that my pipeline finds many proteins, and I want to identify enriched or depleted protein classes (e.g., enrichment for DNA-binding proteins).

By googling for "type IV secretion protein Rhs", I found a UniProt entry which has a link to a QuickGO page, which includes three categories: "transposase activity", "transposition, DNA-mediated" and "DNA binding". From this page I can also get an ancestor chart for each of these categories, e.g.:

Optimally, the code would work quickly (preferably locally), and for a description from a CDS feature (e.g.: "type IV secretion protein Rhs") return the GO categories it belongs to. It would be even better to also have the ancestry of each category in a computer-readable format.

What is the best way to do that?
Is there a specific UniProt subset that fits my need (finding GO categories given a protein description)?
Is there a specific GO subset that fits my need (extracting ancestor GO categories as explained here)?

(I figured that GO categories fit my need of classifying proteins, but if a better option exists, please let me know)

# Answer
> 2 votes
Using the UniProt website API (https://www.uniprot.org/help/api) this could be done with a query like

https://www.uniprot.org/uniprot/?query=name%3A%22type%20IV%20secretion%20protein%20Rhs%22&columns=id%2Centry%20name%2Creviewed%2Cprotein%20names%2Cgenes%2Cgo(biological%20process)%2Cgo(molecular%20function)%2Cgo(cellular%20component)&sort=score

or, in tab-separated format:

https://www.uniprot.org/uniprot/?query=name:%22type%20IV%20secretion%20protein%20Rhs%22&format=tab&limit=10&columns=id,entry%20name,reviewed,protein%20names,genes,go(biological%20process),go(molecular%20function),go(cellular%20component)&sort=score

(Remove the "&limit=10" part from the URL to obtain the full set instead of the first 10 entries.)

You can add colums as described here: https://www.uniprot.org/help/customize, https://www.youtube.com/watch?v=ado1r8IDm3U (video, at about 1min25) or https://insideuniprot.blogspot.com/2015/03

Note that while UniProt tries to improve protein name annotation, the unreviewed section of UniProtKB (i.e. UniProtKB/TrEMBL) is often far from perfect and its accuracy relies on the information provided by the submitter of the nucleotide entry, unless an automatic annotation rule is available to assign a more suitable name (see https://www.uniprot.org/help/protein\_names).

Please don't hesitate to contact the UniProt helpdesk if you have any additional questions.

---
Tags: sequence-annotation, uniprot, genbank, gene-ontology
---


# Question
Title: Bulk RNA-Seq Read Length Normalization across different samples
I have 20 samples out of which 14 are 100 bp in length and 6 are 150 bp. Is there a way to normalize the read length for cross-sample differential expression comparison? What would be the best way to approach this?

Note: The samples are NOT from reference/model species.

# Answer
> 1 votes
The simplest thing to do is to trim the 150 bp fasts so that they are 100 long. I don't think there is an easy way to correct for the fact that the 150 bp long reads will have a higher unambiguous rate of alignment and gene assignment than the 100 bp long reads.

If you have a mix of all the experimental conditions across all the lengths, you can include length as an element in the design, so that DESeq2 may be able to model the difference between the two batches.

# Answer
> 0 votes
I've done this by doing a RPkM-like transformation which I call VSTPk. Transformations that adjust for transcript length will correct a little bit for the noise associated with mapping short reads to transcripts of varying lengths:

https://bioinformatics.stackexchange.com/a/15218/73

Another, probably better, alternative is to use Salmon to map reads to transcripts, which applies its own mapping-based correction to read counts. Salmon produces a statistic that relates to the mappable length of a transcript, and normalises counts based on how many molecules the model predicts came from each transcript isoform.

Unfortunately all these normalisation techniques require good gene models, which are difficult to obtain for non-model organisms.

---
Tags: rna-seq, deseq2, normalization, edger, salmon
---


# Question
Title: How to extract column from a matrix matching the another file (sample file)
I have two files master file.txt and sample.txt

master.txt

```
Name    GTEX.1117F.3226.SM.5N9CT    GTEX.111FC.3126.SM.5GZZ2    GTEX.111FC.3326.SM.5GZYV    GTEX.1128S.2726.SM.5H12C    GTEX.1128S.2826.SM.5N9DI    GTEX.117XS.3026.SM.5N9CA
ENSG00000223972 1   0   0   1   0   0
ENSG00000227232 298 168 197 106 221 184
ENSG00000278267 0   1   1   0   0   0
ENSG00000243485 0   2   0   1   0   1
ENSG00000237613 0   0   0   1   0   1
ENSG00000268020 1   0   0   1   0   1
ENSG00000240361 1   3   0   4   2   1
ENSG00000186092 3   1   1   1   1   2

```

Sample file.txt

```
GTEX.1117F
GTEX.111FC
GTEX.111XS

```

Desired output

```
Name    GTEX.1117F.3226.SM.5N9CT    GTEX.111FC.3126.SM.5GZZ2    GTEX.111FC.3326.SM.5GZYV    GTEX.117XS.3026.SM.5N9CA
ENSG00000223972 1   0   0   0
ENSG00000227232 298 168 197 184
ENSG00000278267 0   1   1   0
ENSG00000243485 0   2   0   1
ENSG00000237613 0   0   0   1
ENSG00000268020 1   0   0   1
ENSG00000240361 1   3   0   1
ENSG00000186092 3   1   1   2
ENSG00000238009 2   2   1   1
ENSG00000233750 0   3   3   0
ENSG00000268903 103 44  76  24

```

I tried using the R-code form this link Print specific columns in a matrix on the basis of sample id's in the header

However, I am not able to figure out how to extract columns matching only half part of ID (like for GTEX.111FC I want both GTEX.111FC.3126.SM.5GZZ2 and GTEX.111FC.3326.SM.5GZYV)

Thank you

# Answer
Assuming your "master" files looks like this:

```
df <- data.frame(Name = LETTERS[1:5],
                 Col1 = sample(1:10,5),
                 Col2 = sample(1:10,5),
                 Col11 = sample(1:10,5),
                 Col3 = sample(1:10,5))

  Name Col1 Col2 Col11 Col3
1    A    5    4     3   10
2    B    1    6     5    6
3    C    3    1     9    4
4    D    9    9     8    9
5    E    8   10     1    1

```

For a sequence of pattern you want to extract, you can use base r function such as `grep`:

```
match_pat = c("Col1","Col2")

m <- unlist(sapply(match_pat,function(x) grep(x,colnames(df))))

Col11 Col12  Col2
    2     4     3

```

And then, extract your columns of interest by doing:

```
> df[,m]
  Col1 Col11 Col2
1    5     3    4
2    1     5    6
3    3     9    1
4    9     8    9
5    8     1   10

```

In your example, it will be something like:

match_pat <- c("GTEX.1117F","GTEX.111FC","GTEX.111XS")
m <- unlist(sapply(match_pat,function(x) grep(x,colnames(df))))

master[,m]
```

> 1 votes

---
Tags: r, data-preprocessing, gtex
---


title: What is a PacBio "movie file"?

text:
# Question
Title: What is a PacBio "movie file"?
I came across references to "movie files" from PacBio sequencing in this paper:

https://www.jimmunol.org/content/204/12/3434

Specifically:

> Movie files used to generate results presented in this article have been submitted to the National Center for Biotechnology Information BioProject (https://www.ncbi.nlm.nih.gov/bioproject/) under accession number PRJNA389440.

Looking at the SRA entries, it looks like it's just the raw reads as I expected. **But what is meant by "movie files"?** (This is probably just my ignorance about PacBio as opposed to Illumina or others but I've never seen that term before and can't find much online.)

# Answer
> 1 votes
\[disclaimer: not a PacBio person, I'm just vaguely aware of the technology\]

Movie files are the PacBio raw data format, representing observed signal intensities in each sequencing well over the course of the sequencing run. The movie files can be used to re-call sequences with updated statistical models, which can improve accuracy, or test out hypotheses about how epigenetics alters the rate of nucleotide incorporation (i.e. methylation).

Illumina has a similar static version of this which is an image (or multiple images, depending on chemistry) of the flow cell after each synthesis cycle. Illumina's raw data format consumes tens of terabytes per run, and is of minimal value after base calling, so is almost always discarded after base calling.

Anything uploaded to NCBI SRA must contain FASTQ information. While the original uploaded files can be kept (which may include time/signal information), the files are converted into a FASTQ-like format for NCBI's internal representation. As such, there's a chance that movie files could be uploaded to SRA, but those files would *also* need to include the called FASTQ information (e.g. in a HDF5 container). In contrast to this, ENA *does* allow uploading raw data files without called information.

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


# Question
Title: Why are TPMs per 10k or 100k in many scRNA-seq studies?
I noticed that many scRNA-seq papers normalize TPMs to 10k or 100k as opposed to 1M (as the abbreviation defines them). It doesn't really matter since you are just moving the decimal point, so why mess with an established convention?

Additionally, the values might be log-transformed. In that case, log(x+1) would affect the data less if the total is 1M.

# Answer
When studying demographics, you often need to account for the population size of a city or region (GDP per capita, is a good example). Normalizing by $10^6$ people is a nice round number, but lots of cities and regions don't have that many people. So if we want to keep numbers in the range \[1, 100\], 1M is too big a denominator. That's why a lot of studies use $10^5$ people as a scaling factor.

The same goes for RNA-seq and other sequencing measurements. Your readout will typically contain $10^7 - 10^8$ reads for a bulk sequencing run. If your sample has thousands of cells, this means very few cells will have \> $10^6$ reads assigned to it.

But the TPM unit was designed for bulk measurements and that's what most people are used to. So we can keep the same idea of scaling by number of reads, but just use a different convenient scaling factor. This way we get both scaling by sequencing depth and nice numbers to work with.

As for your $\log(1+x)$ comment, yes, if $x$ is in TP10K, the low read counts will be biased more than if $x$ was in units of TPM. But consistently lowly-expressed genes often aren't of interest. So this won't be too much of a problem. The bigger issue is the zero count over-inflation that is common in scRNA-seq measurements.
> 2 votes

# Answer
Collecting all the question comments as a community answer (feel free to edit this to make it more readable):

There tend to be far fewer counts in single cell experiments. By reducing the scaling factor you are avoiding having massive expression values for your samples.

It isn't uncommon to scale single cell experiments by the median cell library size. That way small values typically won't be too small and large values won't be too large. People also choose 10,000 because that's roughly around the library size for a single cell.

I haven't seen an explanation yet. Here is an example paper from last week (reputable journal and lab): doi.org/10.1016/j.celrep.2018.11.056 . "... We also excluded outlier cells from the MPPs group which were likely caused by cell contamination. To correct for batch issues, the two experimental batches were centered to have the same mean log(TP100K+1) (hereafter referred to as log(TPM), where TPM stands for ‘transcripts per million’) per gene. PCA was done in R with scaled..."
> 1 votes

---
Tags: rna-seq, scrnaseq, normalization
---


# Question
Title: BAM to gene expression matrix (UMI counts per gene per cell),10X
I am trying to reproduce some results of a scRNASeq experiment. However I am new to the server-side aspect of such analyses and am very confused at the moment.

The data provided by the authors of the paper is in .BAM format and from there I wish to derive a gene expression matrix (UMI counts per gene per cell). The authors stated that they did this using the 10X genomics pipeline, however it seems that .BAM outputs are sort of the "final version" of the pipeline and, the pipeline only really deals with truly raw sequencing files (bcl/fastq). I am trying to learn how to do all this processing myself as I will need to eventually. Unfortunately, most of the documents I have read end up making me more confused. I know BAM files are technically already processed and aligned but, how would I get the counts from the file?

I have considered converting it back to `fastq` format and following the pipeline from the beginning, but it just seems like I would be moving backwards by doing such.

I have also found some R packages that can integrate BAM files in the environment like `Rsamtools`, but I would much rather do this on the server and then upload a csv into R for the downstream analysis as the files are +60GB each.

Essentially I am just asking for some advice, or a vignette/article I can read that would clear up how to proceed with the BAM outputs to derive the gene expression data. I have already read the `cellranger` pages on the 10X genomics website but they were of no help in answering this quesion.

Thank you for your time,

# Answer
> 5 votes
### How to re-analyze 10X BAM files?

This is a great question and honestly, I don't think there was an easy way to do this at the time the question was asked. The reason is that if you want to re-do the authors' analysis to get the gene expression (gene-cell / feature-barcode) matrix and they used `cellranger` (the official 10X pipeline software) to do that, you will need FASTQ files as input - because "10x pipelines require sequencer FASTQs (with embedded barcodes) as input."

**Therefore, to get FASTQ files from 10X BAM files, there's a relatively new tool provided by 10X Genomics: `bamtofastq`. Then you can simply use those FASTQ files as input to `cellranger`.**

More info here: https://support.10xgenomics.com/docs/bamtofastq

### How to correctly download 10X BAM files from NCBI?

The only caveat is that you need a BAM file generated directly by 10X's `cellranger` (or the respective 10X pipeline, if not dealing with gene expression) - that means that a BAM file obtained by downloading an SRA from NCBI and converting to BAM won't work; you need to get the original BAM file directly (often found among the originally submitted files, under "Original format").

To give you some additional help here (because figuring this out is also not trivial): an original 10X BAM file (if originally submitted in this format) can be downloaded from NCBI using `prefetch` from `sra-tools`, specifying the `--type` option - for example:

```
prefetch --type TenX --max-size 100000000 SRR5167880

```

Because the file was larger than the default limit of 20GB, I also had to increase it by specifying the `--max-size` option.

Alternatively, you can follow the official 10X Instructions to Download and Process BAM files of 1.3 Million Brain Cells.

*I know this is an old thread but I'm hoping my answer will be helpful to anyone who has the same question (like me).*

# Answer
> 4 votes
You should be able to parse out what you need using the tags in the .bam. 10xGenomics' website says what tags they add.

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam

Going backwards would also involve parsing the tags, because you have to make two fastq files, and a simple bam -\> fastq pipeline won't do that correctly.

---
Tags: sam, 10x-genomics, cellranger, bash, umi
---