title: Long-form gene summary

text:
# Question
Title: Long-form gene summary
The following (thx to these answers) retrieves the info, in particular a long summary (such as "The membrane-associated protein encoded by this gene is a member of the superfamily of ATP-binding cassette... etc."), for a few mouse genes.

```
library(dplyr)

org.Mm.eg.db::org.Mm.egSYMBOL2EG %>%
  as.data.frame() %>%
  pull(gene_id) %>%
  head() %>%  # !!!
  rentrez::entrez_summary(id = ., db = "gene") %>%
  { do.call("rbind", .) } %>%
  as.data.frame() %>%
  select(name, description, summary)

```

Note the `head`. The protocol/server does not allow large queries. How can I get the same information (incl. the long summary) for all the mouse genes at once?

# Answer
> 1 votes
You could wrap the call too `entrez_summary` in a loop.
For example to cut the data into chunks of five genes:

```
gene.vector # this would be the genes you are interested in
gene.chunks <- ceiling(length(gene.vector) / 5)
# initializing an empty list
summary.list <- vector(mode = "list", length = gene.chunks)
for(i in gene.chunks){
   # select the ith set of five genes from your gene vector
   goi <- gene.vector[(i - 1) * 5 + 1:5]
   # if there are less than five elements there will be NAs that need to be removed
   goi <- goi[!is.na(goi)]
   # adding your code
   summary.list[[i]] <- org.Mm.eg.db::org.Mm.egSYMBOL2EG %>%
     as.data.frame() %>%
     filter(gene_id %in% goi) %>% # here the current five genes are selected
     pull(gene_id) %>%
     rentrez::entrez_summary(id = ., db = "gene") %>%
     { do.call("rbind", .) } %>%
     as.data.frame() %>%
     select(name, description, summary)

   Sys.sleep(10) # wait for 10 seconds
}

```

There will also be a way to do this using `dplyr` `group`. In this case the loop is a little easier to debug.

I added `Sys.sleep` at the end to reduce the frequency with which you are making requests from entrez. If you are posting too many requests in a short time this will be like a DDOS attack. You (and potentially your whole IP range including other people at your institution) could get banned from the service for a couple of hours. Maybe try it with a subset first.

If you really want to get the information for all mouse genes this probably will be rather slow. Alternatively, you can go to the entrez webpage and try to download the information as a list. Unfortunately, I don't know exactly where you will find this information.

# Answer
> 1 votes
The following bash script downloads the summaries into an XML file.

```
#!/usr/bin/env bash

# This requires the ncbi tool `gene2xml`
# sudo apt install ncbi-tools-bin

DIR="download_cache"
mkdir -p "$DIR"

URL="ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/ASN_BINARY/Mammalia/Mus_musculus.ags.gz"
wget -S -O - "$URL" | gzip -d | gene2xml -b T -l T | gzip -9 > "$DIR/mm.xml.gz"

URL="ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/ASN_BINARY/Mammalia/Homo_sapiens.ags.gz"
wget -S -O - "$URL" | gzip -d | gene2xml -b T -l T | gzip -9 > "$DIR/hs.xml.gz"

```

These XML files can be converted to tables along these lines, viz.:

```

import xml.etree.ElementTree as ET
import gzip

from pathlib import Path

from tcga.utils import mkdir, unlist1, relpath

def parse(src):
    def find_text(e, path):
        e = e.findall(path)
        if e:
            return unlist1(e).text
        else:
            return ""

    yield ('gene_id', 'description', 'summary')

    with gzip.open(src, mode='r') as fd:
        e: ET.Element
        for (x, e) in ET.iterparse(fd, events=['end']):
            if e.tag == "Entrezgene":
                gene_id = find_text(e, "**/Gene-track_geneid")
                descptn = find_text(e, "**/Gene-ref_desc")
                summary = find_text(e, "Entrezgene_summary")

                yield (gene_id, descptn, summary)

                e.clear()

def main():
    for src in Path(__file__).parent.glob("download_cache/*.xml.gz"):
        trg = mkdir(Path(__file__).with_suffix('')) / Path(src.stem).with_suffix('.tsv.gz')
        if trg.is_file():
            print("Skipping existing file:", relpath(trg))
        else:
            print("Writing:", relpath(trg))
            with gzip.open(trg, mode='wt') as fd:
                for t in parse(src):
                    print(*t, sep='\t', file=fd)

if __name__ == '__main__':
    main()

```

Btw., a table of vertebrate gene homology is available here (it's just a tab-separated table), and can be used to impute summaries by homology.

---
Tags: r, phylogenetics, entrez
---


# Question
Title: MUMmer plot error. Line 884. What to do?
My third year comp sci friend got me this far, I hope you guys can help me with the last inch! For the record, I have practically no coding experience. I am using the 3.9.4 alpha version of MUMmer.

Here is the manual I have been following: http://mummer.sourceforge.net/examples/

I am at part 2.2.2 Running the Mummer plot, and I type in the final line to generate the plot.

I type:

`mummerplot -x "[0,8666124]" -y "[0,9025688]" -postscript -p mummer mummer.mums`

And the error is:

```
Can't use 'defined(%hash)' (Maybe you should omit the defined ()?) at /usr/bin/mummerplot line 884

```

Line at 884:

```
if ( defined (%$rref) ) {                                                                                #-- skip reference sequence or set atts from hash
        if ( !exists ($rref->{$idR}) ) { next; }
        else { ($refoff, $reflen, $refdir) = @{$rref->{$idR}}; }

```

Its perl

# Answer
> 2 votes
I had the same error and I solved by removing `defined` and writing only `if (%$rref)` in the mummerplot code. You have to do this for all the lines that you have it. There should be six lines to correct in this way .

You can find an explanation of why here, that tells that the function `defined`:

> Returns a Boolean value telling whether EXPR has a value other than the undefined value undef. If EXPR is not present, $\_ is checked.

And then:

> Use of `defined` on aggregates (hashes and arrays) is no longer supported. It used to report whether memory for that aggregate had ever been allocated. You should instead use a simple test for size: `if (%a_hash) { print "has hash members\n" }`

---
Tags: errors
---


# Question
Title: Combining read counts from three separate GEO studies
I want to do differential expression analysis with DESEQ2. I have three read counts files downloaded from GEO (small RNAseq based) where the number of miRNAs and id is nearly the same. These studies have used the same tool for alignment but different tools for counting aligned reads.

My question is

1. Can I combine these read count files
2. Is there any recommended meta-analysis based tool for this purpose
3. Is there any other way to do this other then redoing RNA-seq analysis.

Thank you

# Answer
1. Yes, but there will be a batch effect.
2. Not that I know if, you'll need to use sva or something along those lines to handle the batch effect.
3. Crossing your fingers and hoping the batch effect is small is pretty much your only other option if the data isn't setup in a way that you can use sva or similar. So yeah, redoing the experiment may be the only real option.
> 1 votes

---
Tags: rna-seq, reads, meta-analysis
---


# Question
Title: Automatize classical debugging in a bash script? (normal debugging by hand)
And thanks in advance. I am debugging a text file which is the one the two inputs on the program RGmatch https://bitbucket.org/pfurio/rgmatch/src/master/.

The code works. Is that some records give division by zero.

I am doing it by hand.

But is there a way by making a bash script, to see which records show the problem?

I do not want to touch the python code.

```
python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b ../../test.bed -o myassoSort.txt

```

I have tried with classical debugging, but my actual question is how to automatize that?

Also there are no empty records. The bed file was created based on refseq.txt from the latest human assembly. And also the GTF file. On the refseq.txt I processed to obtain a bed file.

I have checked if this works on a single record, it does, but not in some specific ones. The error I get is:

```
  File "rgmatch.py", line 1311, in <module>
    main()
  File "rgmatch.py", line 287, in main
    run(gtf, dhs, outputfile)
  File "rgmatch.py", line 1153, in run
    pctg_area      = (float(region_overlap)/intron_length)*100
ZeroDivisionError: float division by zero

```

The record I found so far is this one.

`chr19 14087850 14090751`

I did grep on both GTF files and Refseq Files and in both files there is information about the record. Because the program search for areas between region of interest, If I put inside just the record only in GTF file and bed Files no error occurs.

`grep "14090751" refseq.bed`

`692 NM_138352.3 chr19 - 14087850 14090751 14088416 14090420 6 14087850,14088670,14088879,14089067,14089784,14090084, 14088575,14088782,14088976,14089362,14090084,14090751, 0 SAMD1 cmpl cmpl 0,2,1,0,0,0,`

`grep "14090751" hg38.ncbiRefSeq.gtf`

```
chr19   ncbiRefSeq.2021-02-10   exon    14090085    14090751    .   -   .   gene_id "SAMD1"; transcript_id "NM_138352.3"; exon_number "1"; exon_id "NM_138352.3.1"; gene_name "SAMD1";
chr19   ncbiRefSeq.2021-02-10   5UTR    14090421    14090751    .   -   .   gene_id "SAMD1"; transcript_id "NM_138352.3"; exon_number "1"; exon_id "NM_138352.3.1"; gene_name "SAMD1";

```

# Answer
> 1 votes
If you want to just start incrementally adding lines from the bed file until you hit the one that reproduces the problem, you can do something like this:

```
for i in $(seq 20 30 $(wc -l < ../../test.bed)); do
    head -n $i ../../test.bed > slice.bed
    if python rgmatch.py -g ../../hg38.ncbiRefSeq.gtf -b slice.bed -o slice.out
    then
        true
    else
        echo "Failed when running with the $i first lines."
        break
    fi
done

```

**Explanation**

`seq 20 30 $(wc -l < ../../test.bed)` will iterate over all numbers between 20 and the number of lines in your bed file with increments of 30. Next, the `head -n $i ../../test.bed > slice.bed` will take the first `$i` lines of the bed file. This way, the loop will run the `rgmatch` command on the first 20 lines of the file, then the first 50, then the first 80 and so on as long as the command runs successfully. When the `rgmatch` command fails, it will report the value of `$i` to help you find the line that caused the problem.

---
Tags: python, bash
---


# Question
Title: How to specify the list of contigs in AGP gap file used for ENA submissions
I am battling the submission of scaffolded assembly to EBI, somehow I can not get AGP gap file right.

*background*

The only option seems to be using the Webin-CLI submission interface. For the scaffolded assembly submission, one needs the manifest file, the sequence (fasta or flat file), and AGP file, which denotes how contigs are glued together into scaffolds. I got already that the sequence file should contain only contigs, and EBI will glue them together using the AGP file.

*the problem*

So far so good, I checked the specification and the example provided and generated similarly looking files for my assembly, but then when I tried to validate the files using

```
ena-webin-cli -context genome -userName XXX -password XXX -manifest manifest_file -outputDir asm01_submission -inputDir ena_submission_input -ascp -validate

```

I got tons of errors with type

```
ERROR: Scaffold "NODE_667656" has only "1" component, minimum two components expected.
ERROR: Scaffold "NODE_353389" has only "1" component, minimum two components expected.
...

```

supposedly for each contig that actually was not scaffolded. I am quite confused about this, as the example file contains tons of scaffold with one component only (for example the first one `EG1_scaffold1`).

What am I missing?

Sample of my AGP file (incl 1 scaffold and 3 contigs):

```
Path_2003       1       1672    1       W       NODE_32694_length_1672_cov_36.024451    1       1672    +
Path_2003       1673    1772    2       U       100     scaffold        yes     align_trnscpt
Path_2003       1773    2521    3       W       NODE_61726_length_749_cov_25.233631     1       749     +
Path_2003       2522    2621    4       U       100     scaffold        yes     align_trnscpt
Path_2003       2622    5705    5       W       NODE_20424_length_3084_cov_26.290988    1       3084    +
NODE_667656     1       86      1       W       NODE_667656_length_86_cov_136.111111    1       86      +
NODE_353389     1       155     1       W       NODE_353389_length_155_cov_25.230769    1       155     +

```

# Answer
Oh you silly sausage! Don't use AGP gap files.

This is not a direct answer but a workaround. However, if your scaffolded genome is annotated, you might as well submit a flat file (.embl) instead of fasta and the gap file. The thing is, making that file is actually a lot more straightforward thanks to amazing EMBLmyGFF3. You will also need the locus tag (a string used as prefix for loci in your genome) you must pre-register in the web interface (it's among "Project" options). Note that although the script does handle some EBI specifics (like automatic filtering of gff3 tags they don't like, or filtering sequences \< 100 bases), some there are things that need to be taken care of. For example, I had to remove duplicates (using `agat_sp_fix_features_locations_duplicated.pl`) and overlapping genes (using `agat_sp_fix_overlaping_genes.pl`), both available in AGAT collection. This might be different for genomes annotated with a different tool (mine was `braker2`).

Anyway, once I was advised to use this tool, it took like 2 hours for 400Mbp genome including the uploading time, so I certainly recommend it over wrestling with AGP files.
> 1 votes

---
Tags: database, public-databases
---


# Question
Title: Label and labelbox mismatched in Pymol after ray
I would like to ask how to overcome the problem of mismatched boxes and labels (pic). Is there an option focused on this? thanks

Radim

# Answer
It is a glitch in PyMOL.

**Caveat**. I would say that, whereas `cmd.label` is great for adding labels en-mass for internal figures, it is not great for figures for dissemination (which require few *strategically* placed labels possibly with a faint white outer glow —cf. your D330): most figures in papers I would say are Powerpoint or Photoshop labelled.

You called it 'labelbox', this is not a normal term, so you may be using a custom script. The box is normally called a `bounding box` and `connector line`. To add it you need to set a few values. You may have done the following:

```
label name CA and resi 44, 'Foo'
set label_outline_color, 0xff8c69  # font colour to salmon
set label_connector_color, 0xff8c69 # connect or bounding box
set label_bg_outline, true  # adds a border (lowercase bool)
set label_connector, true   # adds a line
set label_connector_mode, 1 # how the line looks label_position
set label_position, [5, 5, 5]  # offset the label
set ray_label_connector_flat, true # cylinders as line? No thanks.

```

The issue is the `set label_position, [5, 5, 5]`. If you give the camera-based reference vector `[0,0,5]` you get no line as Z is depth. This ray traces fine. Otherwise the glitch you found happens. My guess is that the box (but not the label) gets calculated based on the viewport window size, not the ray trace window size. So my guess is that there is a resolution and ray trace window size combination just right that they do align... but on my machine at least this is low res.

To make matters worse, the shadow is horrendous and in a random location and `label_shadow_mode` does not stop it. Let's not even discuss `ray_label_connector_flat`. So the best way to avoid this is to simply label in Photoshop, which is much more convenient anyway.
> 1 votes

---
Tags: visualization, pymol, 3d-structure
---


# Question
Title: Is there an efficient way to extract CIGAR strings for read pairs from bam files with python?
I am working with bam files and I have to check if reads of a specific position or their mates are soft clipped. So, I am looking for a fast way to extract the read pairs from a bam file in python. So far, I use pysam and fetch reads of a given position. From those reads, I have access to the CIGAR string and can check them for soft clipping. My problem is that I cannot access the CIGAR of the read mates. Is there a way to access information of read mates with pysam or another tool, especially CIGAR strings, without iterating over the whole .bam file?

Unlike in this question Extracting the CIGAR string from a BAM via Python? I am looking for the CIGAR of both reads of a pair.

# Answer
> 5 votes
I think that the easiest way is to work with query-name sorted groups of reads. In that case, mates will be adjacent in the sort and you can use that to extract the paired CIGARs.

If you are depending on position sorting to extract the alignments of interest efficiently, you can maybe do the following:

1. extract alignments in the region of interest plus a lot of padding to ensure that you are also getting the mates (this will obviously only work with standard short-insert shotgun libraries).
2. Query-name sort the list iterator using built-in `sort()` or `sorted()` with a key function that returns the query name.
3. Iterate over the pairs using adjacency in the iterator to extract paired CIGARs. Discard unpaired reads (first reads whose query name does not match supposed mates).

Code would look something like this (I have not run or tested this in any way!):

```
import pysam

def get_query_key(read):
  """read is pysam.AlignedSegment"""
  return read.query_name

def get_paired_cigars(bam, chr, start, end, padding=1000):
  """bam is pysam.AlignmentFile"""
  reads = list(bam.fetch(bam, chr, start-padding, end+padding))  # may fail at chr ends
  reads.sort(key=get_query_key)
  first_cigar = None
  mate_cigar = None
  cigars = list()
  qn = ""
  for read in reads:
    if first_cigar is None:
      first_cigar = read.cigarstring
      qn = read.query_name
    elif mate_cigar is None:
      if read.query_name != qn:
        first_cigar = read.cigarstring
        qn = read.query_name
      else:
        mate_cigar = read.cigarstring
    if first_cigar is not None and mate_cigar is not None:
      cigars.append((first_cigar, mate_cigar))
      first_cigar = None
      mate_cigar = None
  return cigars

```

I agree that I would like this to be simpler! Unfortunately BAM files are not terribly well-designed for using pairing information at that level of detail IMO, I have always had to use something like this. I would be interested to see a straightforward alternative to this approach from someone with more experience.

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


# Question
Title: Finding smaller sequences from within larger sequences
I am currently working with fastq files which have hundreds of thousands of lines of text. However, not all of them are sequences I am interested in. My sequences are in one line and have a fixed length of 272 characters.

I possess a list of smaller sequences, around 2 million text strings in the following formats (which are fragments from the various fastq files I have):

Format 1

```
ATCGATCG
ATCGATCGATCG
ACGATCGTGTTACG
ATTTACGTACGTA
AATCG
ACGATACGATACG

```

Or I also have them in Format 2

```
ATCGATCGATACGAT|ATACGTGTTACGAT|ATCGATACG|ATACGATGA

```

(I just wrote random sequences just to give an idea of how the pattern.txt files look like, I know Format 1 and 2 are not the same in this post, however, the files are in the same order - only the syntax is changed.).

What I've been trying to do is use grep in a Linux environment and use the pattern.txt files to search for the bigger sequences in my fastq files where the pattern is present, and then output it into a .txt file to further work on.

I've used this command so far (with the Format 2 of patterns):

```
grep -E -f patterns.txt target.fastq > output.txt

```

**Now, this works.**

**However,** I have to limit the number of patterns in my key file and the number of fastq files I can search at once, otherwise I get "grep: memory exhausted" error.

I also want to print out all the sequences up to and including my pattern, but not past it (I don't care what's after my pattern sequence).

And for that I have **successfully** used the following:

```
grep -oE "\S+singlepattern" target.fastq > output.txt

```

However, I can only manage to get this to work if I manually write in the command the pattern sequence I want to look for (thus I can only use a couple). I can't seem to get it to work with the file list.

So I thought of using a while loop. This is what I've written so far in the terminal (it's one continuous line, but I formatted it here for simplicity on the eyes):

```
while read line; do
grep $line target.fastq > output.txt;
done < patterns.txt

```

**However,** this does absolutely nothing and just spits out a blank output.txt

I have tried changing the syntax, with quotations on arguments, no quotations, I tried the various -E -f -o -F, egrep, etc but it's all the same. No error message, just blank file.

I've also tried implementing the \S to the loop as follows:

```
while read line; do
grep -oE "\S+$line" target.fastq > output.txt;
done < patterns.txt

```

Yet, the result is the same.

Any suggestions on how I can fix this? Or even other methods to get the cut-down sequences (up to and including my pattern) from the fastq files into an output.txt file?

Thanks in advance!

# Answer
> 1 votes
Thank you all for your help. I managed to fix it in the end. The error was with how my pattern.txt file was formatted.

I initially exported the .txt from R Studio by specifying the EOL as \r\n. However, I think this was the issue as then from my laptop to the computer cluster the files seemed to have a windows EOL.

Thus on the linux server I did

```
dos2unix patterns.txt

```

And since then all of the suggested loops, both while read and the for chunk work beautifully and I get a nice `\S+$line` exported into the desired output.txt.

Thanks again all for your advice and help!

# Answer
> 0 votes
The reason your loop is giving a blank file is because you are using `>` which will overwrite whatever was previously in the file. This means that your `output.txt` will only have the result of the last `grep` and if that last `grep` has no results, you'll get an empty file.

Now, you're on the right track here. You could use a loop, although it will be horribly, horribly slow (by the way, you can use this exact format in the terminal as well, it doesn't need to be on the same line):

```
## create output.txt if it doesn't exist and empty it if it does exist
> output.txt

while read line; do
    ## Use >> instead of > so you append to the file
    grep -oE "\S+$line" target.fastq >> output.txt;
done < patterns.txt

```

Your attempt with the `-f` option is the most reasonable. You might get it to work if you tell `grep` to treat the input as fixed strings and not regular expressions (`-F`). Also, you don't need the `-E` (and that might be causing issues), so try this in case it works:

```
grep -Ff patterns.txt target.fastq > output.txt

```

If that doesn't work, just split your `patterns.txt` file into smaller files and then use grep in a loop with each of them:

```
split -l300000 patterns.txt pattern_chunk

```

That will split `patterns.txt` into multiple files of 3000000 lines each, named `pattern_chunkaa`, `pattern_chunkab` etc. For example, this is the result of running that command on a file with 5 million lines:

```
$ wc -l patterns.txt
5000000 patterns.txt
$ split -l300000 patterns.txt pattern_chunk
$ ls
pattern_chunkaa  pattern_chunkad  pattern_chunkag  pattern_chunkaj  pattern_chunkam  pattern_chunkap
pattern_chunkab  pattern_chunkae  pattern_chunkah  pattern_chunkak  pattern_chunkan  pattern_chunkaq
pattern_chunkac  pattern_chunkaf  pattern_chunkai  pattern_chunkal  pattern_chunkao  patterns.txt

```

So, you can now use each of those as the patterns:

```
## create output.txt if it doesn't exist and empty it if it does exist
> output.txt
for chunk in pattern_chunk*; do
    grep -f "$chunk" target.fastq >> output.txt;
done

```

This doesn't help you get the sequence only up until the matched pattern. For that, you can add the `\S+` to the beginning of each pattern in `patterns.txt` (this assumes you are using the file which has one pattern per line) and then split the file and loop over the chunks again:

```
## add the \S+
sed -i 's/^/\\S+/' patterns.txt
## create output.txt if it doesn't exist and empty it if it does exist
> output.txt

for chunk in pattern_chunk*; do
    grep -Ef "$chunk" target.fastq >> output.txt;
done

```

Note how I keep using `> output.txt`. This is not actually necessary, you can use `command > file` even if `file` doesn't exist. However, it is useful here to make sure that the `output.txt` file doesn't have the result of a previous run since we are now appending (`>>`) instead of overwriting (`>`) its contents.

# Answer
> 0 votes
I would recommend not do this using grep. There is a program within the KAT suite called KAT filter for filtering reads by a set of kmers (that might require all your sequences to be of constant length, but I am not sure). Another software you might want to check out is Cookiecutter. These tools specialised in subsetting sequencing files are usually a lot more efficient than generic regular expression tools.

---
Tags: fastq, linux, grep
---


# Question
Title: How to compute a p-value for the similarity measurement of two DNA sequences/contigs?
I'm attempting to come up with a statistical model to determine the p-value for the similarity value $ 0 \le s \le 1 $ of two DNA contigs. The idea is to get a numerical measure to estimate if two contigs are homologous.

I want to account for base composition (and/or other structural features) if possible, but not necessarily because I might want to do the analysis only from the information about the length of the sequences/contigs and largest overlaps. Maybe something inspired from other statistical significance in bioinformatics.

Assume we have two DNA sequences: $a$ and $b$ with the basepair lengths $a\_{\text{len}}$, $b\_{\text{len}}$. Their overlap $o\_{a-b}$ is the number of contiguous bases that appears in both sequences. If we define the similarity $s = \frac{o\_{a-b}}{\min{(a\_{\text{len}}, b\_{\text{len}})}}$ how can we define a p-value for $s$?

# Answer
> 1 votes
Well, here my thoughts:

The equation $s$ somehow tries to represent the ratio $\frac{equality}{equality+inequality}$. But $s$ returns only a value. And to assing the p-value to a single value... is complicated.

How could you approach this?

* Using some kind of Hamming distance for each base pair at the position $n$ ($1 \leqslant n \leqslant len$) between the `sequence i` ($s\_{i}$) and `sequence j` ($s\_{j}$) in order to get a set of binary data:

$$s\_{ij} = \begin{Bmatrix} 1\ if\ seq\_{i\_{n}}=seq\_{j\_{n}} \\ 0\ if\ seq\_{i\_{n}}\neq seq\_{j\_{n}} \end{Bmatrix}$$

With this, you get a set of binomial data (eg $\\{0,1,0,1,0,1,1,1,0,1,1,1,0,...\\}$). Now, you should be able to perform some (Binary?) test and get the p-value.

PD: Consider aligning the two sequences if there is interest in discussing base pair homology.

---
Tags: sequence-analysis, phylogenetics, p-values
---


# Question
Title: How do you convert Raw Alignment Score to Bit Score?
I'm coding a pipeline where I make a lot of pairwise alignments, and I end up with raw alignment scores. But, I really need to look at my results in terms of bit scores.

I know that the formula is:

$′ = (ƛ\* − ()) / (2)$

But, I don't know what `ƛ` or `K` would be. I'm currently using **BLOSUM62**.

I've looked around online for how to look up (and even empirically deduce) lambda and k, and I all I've come across were convoluted explanations.

I also know that this question from almost 3 years ago attempts to address the same topic, but the best response only leads to an NCBI link that requires login credentials.

Any help would be much appreciated.

# Answer
> 2 votes
*Answer from @terdon, converted from comment:*

The link just points to the NCBI ftp site where you can download the source package and look at the code. There is no need for login credentials, you can log in as anonymous. It's just that the version has changed so that answer's URL no longer works. Try this instead:

```
wget https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.11.0+-src.tar.gz

```

Then, extract the files and look at the file `c++/src/algo/blast/core/blast_stat.c` which, again according to the answer you linked to, should have the lamda and k tables at line 589.

---
Tags: python, sequence-alignment, bitscore
---


title: About the log2 fold change

text:
# Question
Title: About the log2 fold change
It seems that we have two calculations of log fold change:

`Actual log2(FC) = log2(mean(Group1)/mean(Group2))`

`Limma's "Log(FC)" = mean(log2(Group1)) - mean(log2(Group2))`

For volcano plots, can we use the second one?

https://www.biostars.org/p/100460/

# Answer
> 1 votes
Yes, you can use the second one for volcano plots, but it might help to understand what it's implying.

The difference between these formulas is in the `mean` calculation. The following equations are identical:

```
log(a/b)

```

and

```
log(a) - log(b)

```

However, in your example, the first equation calculates the mean prior to doing the log, whereas the second example does the log first (i.e. `log(mean(Group1))` vs `mean(log(Group1))`).

As a consequence, the equations are *similar*, but not identical, and that difference depends on the underlying assumptions of the expression model.

In the first equation (what you call 'Actual'), the genes are assumed to follow a normal distribution in linear space, whereas in the second equation they're assumed to follow a normal distribution in *log* space.

Which one of these is correct is a subject of debate. It's pretty much established that gene expression for *all genes together* is log-normal (or, at least, much closer to log-normal than to linear), but whether this applies at an individual gene level is less clear.

In most cases, the results will be similar when looking at individual genes, and the equations can be used interchangeably. Where they differ substantially, there will be large expression differences between genes, so it's unlikely to alter the interpretation of results.

# Answer
> 0 votes
You can use what limma reports for everything, including volcano plots. In the post you link to I describe in a comment why this is done in every package that I know.

---
Tags: r
---


# Question
Title: How to extract subset of protein structure (PDB format) file based on a subsequence of that protein
I looking at a particular protein structure called **2LY4** accessible from RSCB PDB website. The corresponding fasta sequence for that structure is this:

```
>2LY4_1|Chain A|High mobility group protein B1|Homo sapiens (9606)
GKGDPKKPRGKMSSYAFFVQTCREEHKKKHPDASVNFSEFSKKCSERWKTMSAKEKGKFEDMAKADKARYEREMKTYIPPKGE
>2LY4_2|Chain B|Cellular tumor antigen p53|Homo sapiens (9606)
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPL

```

And the PDB format file can be downloaded here (1.7MB). The entire is to large to be pasted here.

What I want to do is to extract the subset of PDB format based on the subsequence in fasta above. Namely `Chain A` starting from 1st residue to 30th residue.

```
GKGDPKKPRGKMSSYAFFVQTCREEHKKKH

```

How can I do that in R or Python?

# Answer
In general, if you simply want to extract that part of the PDB file, you could loop over it (it's plain text) and check the fields you're interested in:

```
with open('2ly4.pdb') as pdb:
    for line in pdb:
        if line[:4] == 'ATOM':
            chain = line[21]
            res_idx = int(line[23:26])
            if chain == 'A' and 1 <= res_idx <= 30:
                print(line[:-1])
        else:
            print(line[:-1]) # alternatively, skip with continue

```

However, this doesn't quite work with the file you've linked because it actually contains ten different conformations (models). This is better handled with a proper parser like prody:

```
import prody

model = prody.parsePDB('2ly4.pdb', model=1)
extract = model.select('chain A resnum 1to30')
prody.writePDB('extract.pdb', extract)

```

This will also allow you to do more complex atom selection.
> 4 votes

# Answer
Here is an example how to do it using another Python library with another selection syntax.
Requires `pip install gemmi`.

```
import gemmi

st = gemmi.read_structure('2ly4.pdb')
subset = gemmi.parse_cid('/*/A/1-30').copy_structure_selection(st)
subset.write_minimal_pdb('output.pdb')

```
> 3 votes

# Answer
If you need to do it only once, for this one PDB file, it's easier to do it manually. When you open the file you can see it has 10 models from NMR. Do you want to extract a subset of one model or of all 10? In this file the residues are numbered according to the sequence. So you find the line where residue 31 begins in the first model:

```
ATOM    500  N   PRO A  31      18.844  -9.699  -4.620  1.00  0.00           N

```

and delete all lines until

```
ENDMDL

```

Then you do the same for other models, or you just remove other models, depending on what you need. You may also want to remove headers – from the beginning of the file until the `CRYST1` record.
> 2 votes

---
Tags: r, python, fasta, protein-structure, pdb
---


title: Total pairs of amino acid substitution

text:
# Question
Title: Total pairs of amino acid substitution
I am reading a blog and it says:

> The numbers for identities and replacements used for calculating the overall alignment score in the expression above are usually presented in the form of a 20 x 20 matrix (20 is the number of the most common amino acids). In total there are 210 possible replacement pairs (residues replacing each other) of amino acids, which includes 190 pairs of different amino acid substitutions + 20 pairs of identical substitutions

My question is how are there 210 replacement pairs. What is the formula being used here? Insights will be appreciated. The link to blog is :

https://proteinstructures.com/sequence/amino-acid-substitutions/

# Answer
> 1 votes
There are 190 pairs of different amino acids: $$ {20\choose 2} = 190$$

and 20 pairs of identical amino acids (AA, TT, WW, etc). Combining these, you get all 210 pairs referenced in the article. Another way of seeing this is dividing the 20x20 matrix into 3 parts: the diagonal, all values above the diagonal, and all values below. We have 20 values on the diagonal with 380 values to split evenly across the upper and lower parts, leaving 190 for each. Therefore, the lower (or upper) part plus the diagonal comes out to 210.

---
Tags: sequence-alignment
---


# Question
Title: How do I increase the sensitivity of Salmon Alevin?
I'm in a little bit of a bind with targeted single-cell sequencing. I'm trying to match up our reads to the targeted amplicon panel (418 targets), and all but one have matched successfully with Salmon Alevin. I have a transcript that produces many 75bp reads that only match the first 28bp of a target, and Alevin is ignoring those reads as a match (the 28bp sequence is genome-unique).

Here's an example read, with matched sequence:

```
qual FFFFFFFFFFFFFFFFFFFFFFFFFFFF,,,::,,:,:,,,,,,,,,,,,,,,,,:,,,,,,,,:,,,,,:,:
read GCTGCATAATCTCTTCATTCCGAGGAGCACCCCACGCGCGCCGCGCCCCCCCGCCCCCGGGGGGGGGGGGGGG
SNPs ...................................G.G.G..G.GC.....CG...CCGGGGGGGGG..GG..
 ref GCTGCATAATCTCTTCATTCCGAGGAGCACCCCACCCCCACCCCCACCCCCACCCCATTCTTAAATTGTTTGG

```

As you may be able to see from the above sequence, the reads for this particular transcript get up to 28bp (sometimes with one or two errors), then the wheels start coming off the sequencer, showing an abundance of Gs. This was a NovaSeq run, and a polyG read is the same as no signal using its two-dye process, so I suspect that the sequencing stops after ~35bp, rather than actually synthesising Gs from thin air.

I have lots of these reads; there are about 400k reads that match the reference transcript (Cxcr4, just in case it matters), and about 250k of those reads have a G homopolymer that's at least 8bp in length.

Any ideas on what to do?

FWIW, the above match was taken from a Bowtie2 alignment I did on the reads. I could use Alevin for cell barcode correction only and do the UMI deduplication manually, but it seems a shame to do that when it's only one transcript (from 418) that's spoiling a fully-automated process.

I did try adding in an additional transcript that included the full sequence of a short-match read, but there seem to be too many errors in the non-matching sequence for any reads to count (at least... I *assume* that's why the counts are not coming through; I'm not quite sure how I can work out what's going on behind the scenes).

Is there any way to tweak Salmon Alevin to get it to accept these shorter matches? If I can change the kmer size, how do I do that?

End trimming is unlikely to help, because that's going to shorten the reads and make them less likely to map to the target sequence. I want a way to *include* these reads as valid hits to the target sequence.

# Answer
So, I think there are a few potential options here. Alevin is using selective-alignment internally to determine the mappings for the reads. So, even if you have a k-mer supporting the mapping, if you have gibberish for the rest of the read, the mapping score is going to be very poor and you're not going to recover that mapping locus. There are some parameters you can mess with to try to increase the sensitivity of mapping, though you should always be aware that upping the sensitivity a lot will mean there is less aggressive filtering early in the mapping process which can slow things down. The main parameter you should tweak to allow lower alignment scores is `--minScoreFraction`; the smaller this value the lower (worse) will be the permitted scores for valid mappings. Additionally, you can combine changes to the `--minScoreFraction` with the `--softclip` option to simply score the non-matching end of read as 0 (rather than a long deletion or large set of mismatches); this will lower the effect of bad read ends.

Another alternative you might consider is to try out alevin-fry for processing this data, and mapping the reads with the `--sketch` flag. This will only check certain structural constraints of the mapping and won't force a particular minimum alignment score. The caveat here is that dropping the requirement for alignment scoring can increase the rate of spurious mappings (reads deriving from elsewhere that spuriously map to some indexed target simply by virtue of sharing some sequence). In transcriptome wide analyses, this can be greatly mitigated by building an appropriate index and using a UMI resolution strategy that complements this expanded index. We've not examined this in detail for targeted single-cell sequencing. However, if you want to give it a try and run into any issues, feel free to drop an issue over on GitHub and we'd be happy to help work through it.
> 3 votes

# Answer
@ATPoint commented that *the default kmer size in Salmon indexing is 31*, which led me to a solution that involved altering the index generation step:

Instead of:

```
salmon index -t merged_salmon_${targetName}.fasta -i salmonIndex_${targetName}

```

I now have:

```
salmon index -t merged_salmon_${targetName}.fasta -i salmonIndex_k11_${targetName} -k 11

```

In other words, I changed the kmer size from its default of 31 to 11, which matches the BD-recommended seed substring length for the standard Bowtie2 approach ("-L 10", but I increased that to 11 because Salmon needs odd kmer sizes). I'm pretty sure that my 28bp-mapping reads won't be picked up with a kmer size of 31. Because this is targeted sequencing, I don't expect that a reduced kmer size will be an issue, but I'll keep a watch out for that.

I'm currently re-running the analysis. If it still doesn't work, then I'll see about adding in the tweaks suggested by @nomad.

Update: still not working, so I'll shift onto score tweaks.
> 0 votes

---
Tags: scrnaseq, transcriptome, salmon
---


# Question
Title: What is a good rule of thumb for the threshold of noise versus signal for RPK in RNA seq?
I have RPK values (RNA seq) and I'm wondering what is a good rule of thumb for what is considered to be noise versus what is considered to be signal? I.e what should I choose as a threshold value for what is considered noise versus not noise? (just as a rough guideline)

# Answer
Signal and noise are best determined at the count level, rather than after normalisation. I don't think there's an easy way to get a signal level from normalised data.

When I wanted to establish the noise level of counts from Illumina reads from lots of different samples, I looked at the total count distribution for genes that had counts in only one sample, and used the maximum value of that (or possibly an elbow of that distribution) for my signal threshold.
> 0 votes

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


# Question
Title: Programmatically retrieve Metadata from SRA Run Selector
I previously asked a question about how to retrieve the Accession List associated with a SRA project. The answer was:

`esearch -db sra -query 'PRJNA491191[bioproject]' | efetch -format runinfo`

where `PRJNA491191` is the bioproject that I'm interested in in this example. How can I do the same for the `metadata` file (the arrow in the picture below is pointing to accession list. It should point to metadata, which is what I'm interested in now)?

# Answer
Depending on exactly which information you seek, you can try:

```
esearch -db sra -query 'PRJNA491191[bioproject]' | efetch -format native -mode xml

```

This will return an XML with additional information that is not present in the `runinfo` table. You can use the Entrez Direct tool `xtract` to parse the output XML and extract the fields of interest into a tabular format.

Alternatively, you can try:

```
esearch -db sra -query 'PRJNA491191[bioproject]' | esummary

```

to retrieve a summary output that includes fewer fields than the `efetch` option shown above. Perhaps all of the fields that are of interest to you are present in the summary output... it will save you time and effort in parsing the XML and retrieving the information you want.

For example, you can use the following to retrieve some fields using `xtract` from the `esummary` output:

```
$ esearch -db sra -query 'PRJNA491191[bioproject]' \
  | esummary \
  | xtract -pattern DocumentSummary -element Summary/Title,Statistics@total_size,Study@acc
GSM3391430: hsf-1(OE) rep3; Caenorhabditis elegans; RNA-Seq  731402855  SRP161795
GSM3391429: hsf-1(OE) rep2; Caenorhabditis elegans; RNA-Seq  703822418  SRP161795
GSM3391428: hsf-1(OE) rep1; Caenorhabditis elegans; RNA-Seq  715482490  SRP161795
...

```
> 3 votes

---
Tags: phylogenetics, sra
---


# Question
Title: How to extract docked ligands from .maegz docking results file sorted by glide gscores in Schrodinger Maestro?
My VLS run completed yielding millions of ligands and I would like to extract top-scored 20000 ligands in a separate file. Incorporating the .maegz output in Maestro is not an option as it fills up all the memory and keeps crashing. Are there any command-line options for this task?

PS: I tried using rdkit in python to parse the output file after converting it to .sdf, however, rdkit skipped some ligands that it could not read.

# Answer
> 1 votes
It sounds like you are saying your file is already sorted? Is this correct?

RDKit Mae supplier is an iterator as is the gzip reader, so does not load all in memory, so that is the best option given that loading all 1e6 into memory crashes your system. What is the fraction of unreadable compounds in RDKit? And does it matter? These are likely failing sanitisation —the `sanitize=False` flag. If for some reason none can be ignored, circumventing this is hard. The `removeHs=True` flag is rarely/never the cause even though the description of several aromatics require protons, e.g. an indole ring is `c1ccc2c(c1)cc[nH]2`, the proton removal happens after loading. If you have an identifier and a list of SMILES you can correct dodgy bond orders etc. But making an unsanitisable molecule sanitary is tricky and rarely worth it as the output molecule is likely rubbish anyway. Plus loading in RDKit will allow you do some additional calculations, such as ligand efficiency (best for sorting anyway), partition coefficient etc etc. As opposed to simply slicing your sorted file.

Alternatively, I believe Open Babel does not load all into memory and has a `splitinto` argument to make the file smaller and can convert gzipped Mae file to SDF.

---
Tags: format-conversion, docking
---


# Question
Title: How to identify the DNA sequence of a gene in the complete genome sequence FASTA format file?
I need to identify the sequence of a gene in the complete genome sequence . I thought it was simple, instead it is not a straightforward task !

My method was the following: I downloaded the FASTA format file of the complete genome sequence of the archaea organism that I am studying, that is CP003685, and then I downloaded the CDS file from NCBI by clicking on "\>Send to\>Coding Sequencing\>FASTA Nucleotide\>Create File". In this way I have a txt file (always FASTA format file) in which I find a sequence for each gene (if I am not wrong).

So for example, I should find, in the complete genome sequence, the sequences of the CDS file (the last that I downloaded), right?

For example, if I consider:

```
>lcl|CP003685.1_cds_AFN03983.1_1007 [locus_tag=PFC_05180] [protein=DNA protection protein DPS] [protein_id=AFN03983.1] [location=complement(940023..940580)] [gbkey=CDS]
ATGCCAGAGCATAATAGGAGATTAGTTGAAAGGACAGGAATAGATGTAGAGAAGCTGTTGGAACTTCTCA
TTAAAGCGGCTGCAGCAGAGTTCACAACGTATTATTACTACACTATACTTAGAAACCATGCAACAGGTCT
GGAAGGAGAAGCAATTAAGGAGATTATCGAAGATGCTAGACTTGAGGACAGAAACCATTTTGAAGCCCTC
GTGCCAAGAATATACGAGCTTGGAGGAGAATTACCAAGGGATATTCGGGAATTTGCAGATCTAGCTTCTT
GTAGAGATGCCTACCTACCAGAAGAGCCCACTATAGAGAACATTCTAAAAGTCCTACTCGAAGCTGAGAG
ATGTGCAGTTGGAGTCTATACAGAAATATGCAACTACACTTTTGGAAAGGATCCAAGGACTTATGATTTG
GCCCTTGCCATATTACATGAAGAGATAGAGCATGAAGCATGGTTTGAAGAGCTCTTGACTGGAAAGCCGA
GTGGTCACTTTAGAAGAGGAAAGCCTGGGGAAAGTCCATATGTCTCAAAGTTTTTGAAAACTAGATAG

```

I should find in the complete genome sequence, this sequence, right ? But this is not. Why?

I suspect that I miss very important passages.

Thank you in advance.

# Answer
> 6 votes
I would just blast it. When blasting locally, you need to first make a database from your genome, so assuming you got the command-line version of blast installed you can do something like

```
makeblastdb -in my_study_genome.fa -dbtype nucl
blastn -max_target_seqs 10 -db my_study_genome.fa -evalue 1e-10 -outfmt 6 -query my_downloaded_gene_of_interest.fasta -out gene_in_my_study_genome.blast

```

blast search for anything that is "similar enough" and also considers both strands, so for example, if there is an SNP in the middle of the gene. Hence even if the searched gene and the reference genome are from different sequencing runs / individuals, blast will still find it. Blast will also handle linebreaks for you (in fasta files the nucleotide sequence is quite often wrapped by 60 or 100 characters a line). Long story short, there is not much that can go wrong with blast (while with full text search you can easily shoot yourself in foot).

-- edit 2 --

As @terdon suggested, potentially if you would locating a transcript sequence (excluding introns), you would do better with software like exonerate or genewise, because they model exon/intron boundaries. However, even in those cases blast will give you an approximate position of the gene.

# Answer
> 5 votes
In your sequence header you will see:

```
[location=complement(940023..940580)]

```

It's a reverse strand gene. Reverse complement the sequence, then search the genome sequence.

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


title: PCA plot in R coloured by sample type

text:
# Question
Title: PCA plot in R coloured by sample type
I'm a biologist, not a programmer so please be gentle.

So I have a dataset that looks like

```
Genes  Patient1   Patient2   Patient3
A          324      433         343
B          431       342        124
Z          232       234        267

```

then I have the sample sheet where it contains sample info like:

---

```
Patient1 - Healthy
Patient2 - Disease
Patient3 - Healthy
(110 samples)

```

I am using:

```
library(ggfortify)
df <- dataset
pca_res <- prcomp(df, scale. = TRUE)

autoplot(pca_res)

```

Then I want to do

```
autoplot(pca_res, data = ?, colour = '?')

```

I wish to use the info from the sample sheet to colour my PCA based on the state (healthy/disease) using the autoplot function. Is there a way to do this?

# Answer
I've found this the way to work with autoplot (ggfortify's prcomp autoplot) - Use a df that has the field you need for color as the first column (to make things easy), then do:

```
autoplot(prcomp(df[,-1], data=df, color='column1_name')

```

If the color-by field is not the first column, do:

```
autoplot(prcomp(df[,!(colnames(df) == "color_by_column_name")], data = df, color = "color_by_column_name")

```

EDIT:

That being said, you're going to need to process your data so patients are along rows instead of column. Transpose your dataset so rows are columns and columns become rows; change patient ID to rownames and Gene name to colnames, then `merge` the patient category dataset with the expression dataset. Once you merge, your dataset should look like:

```
         Category    A      B      Z
Patient1 Healthy     324    431    232
Patient2 Disease     433    342    124
Patient3 Healthy     343    124    267

```

And then,

```
autoplot(prcomp(df[,!(colnames(df) == "Category")]), data = df, color = "Category")

```
> 1 votes

# Answer
Can't say I've ever used autoplot, but this is fairly easy to achieve using "base" ggplot.

One thing autoplot does is put the axis on the same scale reflecting the percentage variation captured by each PC - I think this is a good thing, but most PCA plots you'll see don't do this.

As you haven't given us any info on your metadata, I'll just make a few assumptions - lets make it a data.frame called "colData" with a "Status" column.

prcomp stores the PC scores in "pca\_res$x"

```
# generate some metadata
colData <- data.frame(row.names=rownames(pca_res$x),status=as.factor(sample(c("diseased","healthy"),nrow(pca_res$x),replace = T)))

# get % variation for each PC
pca_res$percentVar <- pca_res$sdev^2/sum(pca_res$sdev^2)

# change scale to reflect % variation
d <- t(data.frame(t(pca_res$x)*pca_res$percentVar))

# merge metadata and modified PC scores by row.names (by=0)
d <- merge(d,colData,by=0)

# make the base plot (PC1 vs PC2)
g <- ggplot(d,aes(x=PC1,y=PC2,colour=status))
# add the points
g <- g + geom_point(size=2)
# add a "wow" factor to the points
geom_point(size=4,alpha=0.5)
# display the plot
g

```

There may be a way to do this with autoplot - e.g.

```
autoplot(pca_res,colour=as.numeric(colData$status))

```

gets you part of the way there, but I think using the standard ggplot syntax is easy enough anyway.
> 2 votes

---
Tags: r, ggplot2, pca
---


title: Extracting WDL map keys as a task

text:
# Question
Title: Extracting WDL map keys as a task
I want to get the array of keys from a map data structure, so I can perform array-based functions against the map in subsequent tasks. I am limited to `draft-2` Workflow Description Language (WDL), so I don't have access to the newer spec's extraction functions. I thought I could achieve this by writing a task:

```
task extract_keys {
  Map[String, String] map

  command {
    cut -f1 ${write_map(map)}
  }

  output {
    Array[String] keys = read_lines(stdout())
  }
}

```

This task works just fine and, as far as I can tell, returns the array of map keys. However, when I try to feed this array into subsequent tasks, I get an error:

```
workflow test {
  Map[String, String] my_map = {
    "a": "foo",
    "b": "bar",
    "c": "quux"
  }

  call extract_keys { input: map = my_map }

  scatter (k in extract_keys.keys) {
    call task_foo { input: k = k, s = my_map[k] }
  }
}

```

> Call input and runtime attributes evaluation failed for task\_foo

I don't have access to the full error message or the logs, unfortunately. The frontend I'm using hides these away.

(Incidentally, I'm aware that you can scatter over maps, although it's undocumented, per https://github.com/openwdl/wdl/issues/106#issuecomment-356047538. This code is just illustrative; I wish to do other things, such as `length` and `cross`, that *only* accept arrays.)

Is it not possible to scatter over the output of a previous task?

# Answer
> 4 votes
I seem to have solved this myself by virtue of the linked undocumented feature (scattering over a map). It turns out you don't need a separate task at all. This works:

```
workflow test {
  Map[String, String] my_map = {
    "a": "foo",
    "b": "bar",
    "c": "quux"
  }

  scatter (pairs in my_map) {
    String keys = pairs.left
    # String values = pairs.right
  }

  scatter (k in keys) {
    call task_foo { input: k = k, s = my_map[k] }
  }
}

```

---
Tags: wdl
---


# Question
Title: Positive marker and negative marker in scRNA-seq datasets
I am learning about the marker genes and clustering by Seurat in scRNA-seq datasets. However, I am confused about the term positive marker and negative marker. What I understand is the marker genes of a cluster is the genes which are highly expressed (high counts) in the cells of that cluster. Let's say cluster A has 200 genes. My questions are:

Can I say that those 200 genes are marker genes for cluster A? How to determine which marker genes are positive marker genes and which are negative marker genes?

# Answer
> 1 votes
In scRNA-seq clusters are groups of cells with similar transcriptomes typically defined by graph based methods (e.g. leiden, louvain) or more rarely k-means.

Marker genes could be genes that are known to be highly expressed (or better specifically expressed) in a particular cell type (e.g. CD3E is a T cell marker). You can visualize the expression of such marker genes on a UMAP/t-SNE without doing any clustering.

When talking about positive and negative markers, I assume the reference is too differentially expressed genes between different clusters. Typically simple differential expression tests are conducted (e.g. Wilcoxon rank sum) between each individual cluster and the other collective clusters. The upregulated differentially expressed genes could be considered positive markers and downregulated genes negative markers.

# Answer
> 0 votes
A positive marker is a gene that is overexpressed in a cell population compared to other populations.

A negative marker is a gene that is underexpressed in a cell population compared to other populations.

I think the Seurat FindMarkers function has an option to output Log2FoldChange as well as the marker names. In this case the positive markers will have a positive Log2FoldChange, and the negative markers will have a negative Log2FoldChange.

Both positive and negative markers can be used to describe a cell population. In the research that I've done in the past, there's a particular type of dendritic cell that is only distinguished from another dendritic cell type based on cell surface markers that it *doesn't* have; our own name for these DCs was "Triple negative".

---
Tags: scrnaseq, gene, seurat, clustering
---


# Question
Title: How do I remove allele annotations from SNP Ids in .bim file?
Hello every one i need help in editting my .bim file. So my .bim file looked like following. I want to convert the ids with chromosomal position from `chr1_867635_C_T` to `chr1_867635` while others remained as it is. I want to remove all the alleles from my SNP IDs, not just specific ones.

```
1       rs2880024             866893  A       G
1       chr1_867635_C_T       867635  A       G
1       chr1_869303_C_T       869303  A       G
1       chr1_873558_G_T       873558  C       A
1       exm2275405           879439  A       G

```

# Answer
> 4 votes
The simplest approach is to remove the first occurrence of a `_` followed by a capital letter, another `_` and another capital letter from each line:

```
$ sed 's/_[A-Z]_[A-Z]//' file
1       rs2880024             866893  A       G
1       chr1_867635       867635  A       G
1       chr1_869303       869303  A       G
1       chr1_873558       873558  C       A
1       exm2275405           879439  A       G

```

Alternatively, find the first occurrence of `chr` followed by numbers, then an underscore and more numbers and remove all non-whitespace characters after the last numbers (note that this assumes GNU `sed`, the default on Linux systems):

```
$ sed -E 's/(\schr[0-9]+_[0-9]+)_\S+/\1/' file
1       rs2880024             866893  A       G
1       chr1_867635       867635  A       G
1       chr1_869303       869303  A       G
1       chr1_873558       873558  C       A
1       exm2275405           879439  A       G

```

If your file is more varied than you show, you can use this to ensure you only make the change on the 2nd field:

```
$ awk -F'[ ]' '{ sub(/_[A-Z]_[A-Z]/,"",$2)}1' file
1       rs2880024             866893  A       G
1       chr1_867635_C_T       867635  A       G
1       chr1_869303_C_T       869303  A       G
1       chr1_873558_G_T       873558  C       A
1       exm2275405           879439  A       G

```

Or, the same but only matching A, C, T or G:

```
$ awk -F'[ ]' '{ sub(/_[ACTG]_[ACTG]/,"",$2)}1' file
1       rs2880024             866893  A       G
1       chr1_867635_C_T       867635  A       G
1       chr1_869303_C_T       869303  A       G
1       chr1_873558_G_T       873558  C       A
1       exm2275405           879439  A       G

```

---
Tags: snp, plink, linux
---


# Question
Title: Error: "fewer tokens than expected" while using --update-name flag in plink
I am trying to update my rsIDs in .bim file by using --update-name flag in plink but getting this error message. what could be the reason of it?

```
plink --bfile updated_ff  --update-name updated_ref.csv --make-bed --out test_updatedrs

PLINK v1.90b6.21 64-bit (19 Oct 2020)
Options in effect:
  --bfile updated_ff
  --make-bed
  --out test_updatedrs
  --update-name updated_ref.csv

Hostname: login1
Working directory: /gpfs/home/bsc08/bsc08856/phs001273_dbgap
Start time: Thu Jun 10 14:23:44 2021

Random number seed: 1623327824
96310 MB RAM detected; reserving 48155 MB for main workspace.
Allocated 15236 MB successfully, after larger attempt(s) failed.
Error: Line 1 of updated_ref.csv has fewer tokens than expected.

```

My updated\_ref.csv file looks like this

```
Name,chrpos
rs3793692,chr10_100008436
rs1983865,chr10_100016339
newrs150271620,chr10_100017402
rs10736126,chr10_100025924
chr10_100044928_A_G,chr10_100044928
rs7905091,chr10_100046078
rs6584182,chr10_100055644
chr10_100058710_C_T,chr10_100058710
rs4360594,chr10_10007547
rs927739,chr10_100077732
rs4917817,chr10_100085888

```

while my .bim looks like this

```
1 rs3094315 0 752566 G A
1 rs3131972 0 752721 A G
1 chr1770377 0 770377 T A
1 rs11240777 0 798959 A G
1 chr1834753 0 834753 C G
1 rs72890788 0 836924 A G
1 rs4970383 0 838555 A C
1 chr1840753 0 840753 G A

1 rs4475691 0 846808 A G

```

# Answer
The reason is (likely) that plink is expecting a space separated file for the `--update-name` argument:

```
--update-name <filename> [new ID col. number] [old ID col.] [skip]

```

You've given it a comma-separated file, so it thinks there is only one column, hence the error. You also don't want column names.

Change the sep by:

```
sed 's/,/ /g' updated_ref.csv > updated_ref.txt

```

Or, edit the file in place, and keep a backup of the original named `updated_ref.csv.bak` :

```
sed -i.bak 's/,/ /g' updated_ref.csv > updated_ref.txt

```
> 2 votes

---
Tags: plink
---


# Question
Title: Find selective markers for one cluster automatically
I detected all markers for each cluster and checked 20 top markers for my clusters one by one. Is there any way that I figure out the selective markers automatically from all the clusters?

For instance, in seurat tutorial, I found comparing one cluster two or three other clusters, but I want to see for instance what are the 10 top markers or genes that are expressed only in cluster of interest.

# Answer
> 2 votes
Directly copy-pasting from one of the Seurat vignettes:

```
# find markers for every cluster compared to all remaining cells, report only the positive ones
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>%
    group_by(cluster) %>%
    top_n(n = 2, wt = avg_log2FC)

```

Change `n` to 10 in the `top_n()` call and filter for your cluster of interest and you will get what you want.

---
Tags: seurat, clustering
---


# Question
Title: Filter rows of VCF file for Match=EXACT?
How do I remove rows in a VCF file on 1kGenome column, where Match=EXACT using bash ?

# Answer
> 2 votes
You can filter records in Python for example using `pyvcf.VcfFrame.filter_flagany()` method I wrote:

```
>>> from fuc import pyvcf
>>> data = {
...     'CHROM': ['chr1', 'chr1', 'chr1', 'chr1'],
...     'POS': [100, 101, 102, 103],
...     'ID': ['.', '.', '.', '.'],
...     'REF': ['G', 'T', 'A', 'C'],
...     'ALT': ['A', 'C', 'T', 'A'],
...     'QUAL': ['.', '.', '.', '.'],
...     'FILTER': ['.', '.', '.', '.'],
...     'INFO': ['DB', 'DB;H2;Match=EXACT', 'DB;H2;Match=EXACT', '.'],
...     'FORMAT': ['GT', 'GT', 'GT', 'GT'],
...     'Steven': ['0/0', '0/1', '0/1', '0/0'],
... }
>>> vf = pyvcf.VcfFrame.from_dict([], data)
>>> # vf = pyvcf.VcfFrame.from_file('your_vcf.vcf')
>>> vf.df
  CHROM  POS ID REF ALT QUAL FILTER               INFO FORMAT Steven
0  chr1  100  .   G   A    .      .                 DB     GT    0/0
1  chr1  101  .   T   C    .      .  DB;H2;Match=EXACT     GT    0/1
2  chr1  102  .   A   T    .      .  DB;H2;Match=EXACT     GT    0/1
3  chr1  103  .   C   A    .      .                  .     GT    0/0
>>> vf.filter_flagany(['Match=EXACT'], opposite=True).df
  CHROM  POS ID REF ALT QUAL FILTER INFO FORMAT Steven
0  chr1  100  .   G   A    .      .   DB     GT    0/0
1  chr1  103  .   C   A    .      .    .     GT    0/0
>>> # vf.to_file('filtered.vcf')

```

# Answer
> 1 votes
Hard to say without seeing VCF rows, but generally speaking with text files you can use `grep`, if I understand the question:

```
grep -vw "Match=EXACT" file.vcf > filtered.vcf

```

---
Tags: vcf, ngs, bash, filtering, data-preprocessing
---


# Question
Title: Theoretical limit of human genome compression
How small can a compressed file containing the human genome be? I'm aware that this question cannot actually be answered, since it is asking for the Kolmogorov complexity of the human genome, which is not computable. So reformulating: What is the largest known lower bound that we can put on the size of a file containing the human genome?

An example answer could be: on average we know that the human genome contain at least this amount X of randomness, therefore it is not possible to do better than Y bytes.

I'm agnostic on which human genome, but I would say without polymorphisms to make it simpler.

The compressed sequence would need to be stand alone, without using a reference (I asked a question on that topic here: Compressing the human genome to few megabytes).

There should be a program able to compress and decompress the genome in a lossless way. Apart from that, any kind of modelling and coding is allowed.

# Answer
> 9 votes
For the most modern CHM13 assembly (v2.0; Complete T2T reconstruction of a human genome; https://github.com/marbl/CHM13), the best known compression without any use of a reference is 567 MB (1.45 bits/bp; 1.80 bits/distinct canonical 31-mer).

Upper estimates of genome compressibility can be made using standard general or specialized compressors (see, e.g., this comparison). In my experience, `xz` can achieve great compression ratios on genomic sequences, especially when used with the maximum compression level (`-9`) and a single thread (`-T1`). The compression performance can be further boosted by using a large dictionary (`-e --lzma2=preset=9,dict=512MiB`). Importantly, the original FASTA file needs to be cleaned – linebreaks hamper compression as they act as noise; the reformatting can be done by `seqtk seq -U` (which also converts all nucleotides to the upper case). It the case of the CHM13v2 assembly, `xz` compresses the human genome to 634 MB.

However, even better performance can be achieved with specialized compressors, especially if they have particular support for reverse complements and mutational models. Among the best ones is GeCo3, which has been shown to compress CHM13v2, with a particular combination of parameters, to 567 MB .

Lower estimates are impossible in general (we can't reject the existence of a program that would create the whole sequence from small amount of starting conditions), but can be done for particular compression techniques, e.g., particular variants of dictionary compressors. However, this leads to intractable optimization problems (see e.g. the problem of minimal string attractors).

A useful auxiliary measure for estimating the compressibility of single genomes is the number of distinct k-mers for a "reasonably" selected k (e.g., k=31). The lower bound on the compressibility of k-mer sets can then be roughly approximated by 2 bits per distinct k-mer, however, this is limited by the embedded structural assumptions, e.g., that single genomes are in principle random strings with occasional (possibly distant) repeats. For instance, the CHM13v2 human genome assembly contains 2,512,390,070 distinct canonical 31-mers, which would then correspond to 628 MB assuming 2bits/kmers, while the best available specialized compressors can compress it down to 567 MB, i.e., by approx 10% better.

**Updates**

* May 27, 2022 – updated the upper bound based on the latest results with CHM13v2 (xz & GeCo3), changed the lower bound wording to an auxiliary estimate as the assumptions have shown not to hold

# Answer
> 5 votes
There's a benchmark that compared dozens of compressors on DNA data, including human genome: http://kirr.dyndns.org/sequence-compression-benchmark/

In particular, this table shows file size of compressed human genome, produced by the strongest setting of each compressor: http://kirr.dyndns.org/sequence-compression-benchmark/?d=Homo+sapiens+GCA\_000001405.28+%283.31+GB%29&cs=1&cg=1&com=yes&src=all&nt=4&only-best=1&bn=1&bm=size&sm=same&tm0=name&tm1=size&button=Show+table

The best result is currently 657 MB on that particular genome build (GCA\_000001405.28), when compressed with fastqz.

By the way, xz (mentioned in another answer) is currently number 20 in that table. (Though the genome is un-modified straight from the database, not the "cleaned" form.)

Disclosure: I worked on that benchmark.

Note: What is being compressed here is a FASTA file containing a particular assembly of a human genome. Compared to an actual genome it has this extra content: Sequence names and brakes between contigs, end of line characters, N characters representing unknown bases, and soft-masked repeats (represented as lowercase). On the other hand, it also misses some information from the actual genome: 1. It's an unphased haploid sequence, while the real genome is diploid (two copies of each chromosome). 2. It includes no epigenetic information (DNA methylation, histone modifications), nucleosome positioning, 3D structure of packed DNA, etc. So, it's just a digital model of a genome with a certain degree of accuracy, that is useful, simple, and relatively affordable to obtain. It's a format commonly used for representing genomes in databases.

---
Tags: phylogenetics, human-genome
---


# Question
Title: How to get KEGG pathway names from a list of 'ko' IDs?
I have a list of ko IDs from which I want to get the pathway names. How can I get that? I have tried this link, but it does not give me any result. What to do?

```
  ko00010   ko00030 ko00061 ko00120 ko00121 ko00190 ko00195 ko00230 ko00250 ko00253 ko00270 ko00300 ko00330 ko00400 ko00401 ko00410 ko00430 ko00450 ko00460 ko00471 ko00472 ko00500 ko00511 ko00520 ko00521 ko00523 ko00531 ko00550 ko00562 ko00604 ko00620 ko00625 ko00640

```

Thanks

# Answer
You can use the KEGG API:

```
ko_list="ko00010 ko00030 ko00061 ko00120 ko00121
         ko00190 ko00195 ko00230 ko00250 ko00253
         ko00270 ko00300 ko00330 ko00400 ko00401
         ko00410 ko00430 ko00450 ko00460 ko00471
         ko00472 ko00500 ko00511 ko00520 ko00521
         ko00523 ko00531 ko00550 ko00562 ko00604
         ko00620 ko00625 ko00640"

for i in ${ko_list}
do
    hit=$(curl -# http://rest.kegg.jp/list/pathway | grep -w "${i/ko/map}")
    echo -e "${i}\t${hit}" | cut -f1,3 >> ko_pathways.txt
done

```

Or if your list is too long you can change the input and use a file with a list of ids.

```
while read line
do
    hit=$(curl -# http://rest.kegg.jp/list/pathway | grep -w "${line/ko/map}")
    echo -e "${line}\t${hit}" | cut -f1,3 >> ko_pathways.txt
done < ko_list.txt

```

The result will be something like:

```
 ko00010    Glycolysis / Gluconeogenesis
 ko00030    Pentose phosphate pathway
 ko00061    Fatty acid biosynthesis
 ko00120    Primary bile acid biosynthesis
 ko00121    Secondary bile acid biosynthesis

```
> 2 votes

---
Tags: kegg
---


# Question
Title: querying metadata (GDC Data Portal) in python3
I am trying to query metadata from the Genomic Data Commons (GDC) Data Portal and used the sample code from its website:

```
import requests
import json

cases_endpt = 'https://api.gdc.cancer.gov/cases'

# The 'fields' parameter is passed as a comma-separated string of single names.
fields = ["submitter_id",
          "case_id",
          "primary_site",
          "disease_type",
          "diagnoses.vital_status"
         ]

fields = ','.join(fields)

params = {"fields": fields,
          "format": "TSV",
          "size": "100"
         }

response = requests.get(cases_endpt, params = params)

print(response.content)

```

However, I am not getting a tsv file in return. I am ultimately getting an FPKM instead. Any help is appreciated!

# Answer
> 1 votes
You can use `str()` to convert bytes to strings and `print()` with the proper separator.

```
import requests
from io import StringIO
import pandas as pd

cases_endpt = 'https://api.gdc.cancer.gov/cases'

fields = ["submitter_id",
          "case_id",
          "primary_site",
          "disease_type",
          "diagnoses.vital_status"
          ]

fields = ','.join(fields)

params = {"fields": fields,
          "format": "TSV",
          "size": "100"
          }

response = requests.get(cases_endpt, params=params)

content = str(response.content, 'utf-8')

print(content, sep='\t', end='\r')

```

If you want to convert to a data frame you can use `StringIO()` to load a data frame from an input string directly to a pandas data frame.

```
data = StringIO(content)
df = pd.read_csv(data, sep='\t')

```

The output:

|  | case\_id | disease\_type | id | primary\_site | submitter\_id |
| --- | --- | --- | --- | --- | --- |
| 0 | 375436b3-66ac-4d5e-b495-18a96d812a69 | Adenomas and Adenocarcinomas | 375436b3-66ac-4d5e-b495-18a96d812a69 | Rectum | TCGA-F5-6810 |
| 1 | 38f93ca8-4971-4932-b592-0a17eaece2ad | Ductal and Lobular Neoplasms | 38f93ca8-4971-4932-b592-0a17eaece2ad | Pancreas | C3L-02115 |
| 2 | 5ea1361a-050d-4fa0-be59-b530c456b6e3 | Ductal and Lobular Neoplasms | 5ea1361a-050d-4fa0-be59-b530c456b6e3 | Pancreas | C3L-00928 |
| 3 | 690378f3-e38a-4b25-8489-05a3b4bd40b6 | Adenomas and Adenocarcinomas | 690378f3-e38a-4b25-8489-05a3b4bd40b6 | Kidney | C3N-01175 |
| 4 | 6b4f086e-4ef9-4f30-874e-8a155b2e3bab | Squamous Cell Neoplasms | 6b4f086e-4ef9-4f30-874e-8a155b2e3bab | Other and ill-defined sites | C3N-04276 |

---
Tags: python
---


# Question
Title: Convert SNP data from GEO into PED and MAP files?
I can easily grab SNP-Chip data from NCBI's GEO database for a given Illumina SNP-Chip array, e.g. from here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52147

GEO provides data as SOFT formatted family file(s), MINiML formatted family file(s) or Series Matrix File(s).

I'm looking for a way to convert these into PED and MAP files for use in PLINK...

Or does PLINK read these formats natively?

# Answer
Sadly the 'contact us' email at GEO didn't help at all.

I ended up probably re-inventing the wheel with this script: https://gist.github.com/CholoTook/60968e3ab6d90cb8fd19be55a25592f1

YMMV
> 2 votes

---
Tags: plink, format-conversion, geo, snp-chip
---


title: Filtering VEP annotation file

text:
# Question
Title: Filtering VEP annotation file
I am filtering a VEP annotated vcf, trying to maintain just those variants classified as deleterious by `SIFT`and as damaging (probably or possibly included) by `PolyPhen`.

I am using:

```
filter_vep -i "$input" -o "$output" -filter "SIFT = deleterious and PolyPhen match damaging"

```

and I also tried:

```
filter_vep -i "$input" -o "$output" -filter "SIFT match deleterious and PolyPhen match damaging"

```

The output contains just variants with some value in these two fields, but some of the variants are classified as tolerated by `SIFT` and as benign by `PolyPhen`, so do you know why the filter is not working?

As an example: this two variants shouldn't exist in my filtered VCF, because the first one is classified as benign by PolyPhen and the second is classified as benign by both fields:

```
    CHROM  POS  REF  ALT  FILTER  INFO
    1       1934441 .       G       T       .       PASS    deleterious(0.01)|benign(0.138)
    1       207085100       .       G       A       .       PASS   tolerated(0.16)|benign(0.037)

```

As an example of a variant that has pass properly those filters:

```
    CHROM  POS  REF  ALT  FILTER  INFO
    2       3938449 .       T       A       .       PASS    deleterious(0.01)|probably damaging(0.67)

```

Those are 3 variants among those that I have in my filtered VCF, but I should be getting just those variants which are classified as deleterious and damaging by both SIFT and PolyPhen, as the third variant I've shown.

# Answer
> 4 votes
I think what you're trying to do could be achieved by adding parenthesis to the conditions and by specifying the format of the input file by using the flag `--format vcf`, as specified in the documentation. The final command should look like this:

```
filter_vep -i "$input" --format vcf --filter "(SIFT = deleterious) and (PolyPhen match damaging)" -o "$output"

```

# Answer
> 2 votes
Here is another way to achieve your goal using the `fuc vcf_vep` command I wrote:

```
$ fuc vcf_vep in.vcf 'SIFT.str.contains("deleterious") and PolyPhen.str.contains("damaging")' > out.vcf

```

For getting help:

```
$ fuc vcf_vep -h
usage: fuc vcf_vep [-h] [--opposite] [--as_zero] vcf expr

This command will filter a VCF file annotated by Ensemble VEP. It essentially wraps the `pandas.DataFrame.query` method. For details on query expression, please visit the method's documentation page (https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.query.html#pandas-dataframe-query).

examples:
  $ fuc vcf_vep in.vcf 'SYMBOL == "TP53"' > out.vcf
  $ fuc vcf_vep in.vcf 'SYMBOL != "TP53"' > out.vcf
  $ fuc vcf_vep in.vcf 'SYMBOL == "TP53"' --opposite > out.vcf
  $ fuc vcf_vep in.vcf 'Consequence in ["splice_donor_variant", "stop_gained"]' > out.vcf
  $ fuc vcf_vep in.vcf '(SYMBOL == "TP53") and (Consequence.str.contains("stop_gained"))' > out.vcf
  $ fuc vcf_vep in.vcf 'gnomAD_AF < 0.001' > out.vcf
  $ fuc vcf_vep in.vcf 'gnomAD_AF < 0.001' --as_zero > out.vcf

positional arguments:
  vcf         Ensemble VEP-annotated VCF file
  expr        query expression to evaluate

optional arguments:
  -h, --help  show this help message and exit
  --opposite  use this flag to return records that don’t meet the said criteria
  --as_zero   use this flag to treat missing values as zero instead of NaN

```

You can also do above with Python:

```
from fuc import pyvcf, pyvep
vf = pyvcf.VcfFrame.from_file('in.vcf')
expr = 'SIFT.str.contains("deleterious") and PolyPhen.str.contains("damaging")'
filtered_vf = pyvep.filter_query(vf, expr)
filtered_vf.to_file('out.vcf')

```

# Answer
> -1 votes
If you are working under a Linux environment, then you may easily use a piped grep search. For example:

```
grep deleterious file.vcf |grep damaging >filtered_file.vcf

```

---
Tags: sequence-annotation, variant-calling
---


# Question
Title: Understanding "Measurement of mRNA abundance using RNA-seq data" by Wagner et al 2012
I'm trying to understand why RPKM is not an appropriate way to normalize for RNA seq (I understand the general idea, but I'd like to gain a deeper understanding). So I'm reading the original paper by Wagner et al 2012 that suggests TPM as an alternative. However, I'm struggling to understand the following paragraph, and was wondering if anyone can shed some light on what the authors are trying to say:

> The reason for the inconsistency of RPKM across samples arises from the normalization by the total number of reads. While rmc as well as qPCR results are ratios of transcript concentrations, the RPKM normalizes a proxy for transcript number by $r\_g \times 10^3 / fl\_g$ the number of sequencing reads in millions, $R / 10^6$. The latter, however, is not a measure of total transcript number. The relationship between $R$ and the total number of transcripts sampled depends on the size distribution of RNA transcripts, which can differ between samples. In a sample with, on average, longer transcripts the same number of reads represents fewer transcripts.

One part that I'm particularly confused by is: "The RPKM normalizes a proxy for transcript number by $r\_g \times 10^3/fl\_g$ the number of sequencing reads in millions, $R/10^6$." What is the "proxy" here?

# Answer
That is why: https://www.biostars.org/p/9465851/#9465854.

All these naive per-million scaling methods may fail to correct for library composition differences.
> 1 votes

---
Tags: rna-seq
---


# Question
Title: How to tell if our ligand-protein docking is good from AutoDock Vina's result
I have perform a ligand-protein docking using **Autodock Vina**. The result of the docking looks like this:

```
WARNING: The search space volume > 27000 Angstrom^3 (See FAQ)
Detected 8 CPUs
Reading input ... done.
Setting up the scoring function ... done.
Analyzing the binding site ... done.
Using random seed: -1553787135
Performing search ... done.
Refining results ... done.

mode |   affinity | dist from best mode
     | (kcal/mol) | rmsd l.b.| rmsd u.b.
-----+------------+----------+----------
   1         -5.9      0.000      0.000.  Pose 1
   2         -5.7     22.945     25.492.  Pose 2
   3         -5.5      1.426      2.046.  Pose 3
   4         -5.5     23.669     25.616
   5         -5.4     25.783     29.152.  .....
   6         -5.3     21.146     23.357
   7         -5.2     20.323     22.545
   8         -5.2     23.864     26.064
   9         -5.1     23.422     26.585.  Pose 9

```

As far as I understand from these statistics Mode 1(Pose 1) is the best. However when I actually visualize them in Pymol, Pose 1 has no hydrogen bonding at all but Pose 2 has.

My question is how can we judge if which of those two Pose is the best to use?

Note in figure below Pose 2 has dashed line (Hydrogen bond).

# Answer
It is best to contextualise the numbers. -1 kcal/mol is about the potential energy gained from a hydrogen bond —technically described in the r^6 part of the Lenard–Jones term, it is also the average collision energy of a water molecule at 37°C as that is RT ($\frac{k\_b\cdot T}{N\_A}$, wiki)under a Maxwell–Boltzmann distribution. A salt bridge –2 kcal/mol (Columbic force term). So your scores are not very low, hence why you are counting two hydrogen bonds. Although you can also see a lovely sulfur-pi interaction which is –1 to –2 kcal/mol, so those bonds alone are probably making a –4 kcal/mol contribution, so I am guessing that some terms may be horrendous, such as repulsion forces etc. A nice metric is doing a conversion to ligand efficiency, which weeds out affinities driven by size of the molecule... in this case most of the molecule is doing nothing. Also, most programs have an accuracy of 1 kcal/mol or higher.

So one cannot say what ∆∆G is the best with the data at hand due to noise, but one can say that the ∆∆G is not low enough... Sorry.
> 2 votes

---
Tags: protein-structure, docking, structural-biology
---


title: Seurat - subsetting by genes expressed

text:
# Question
Title: Seurat - subsetting by genes expressed
Originally- I was looking if in at least one of these genes were expressed:

```
test1 <- subset(myData,
                  (KLRC1 > 0) | (KIR2DL1 > 0) |
                (KIR2DL3 > 0) | (KIR2DL4 > 0) |
                (KIR3DL1 > 0) | (KIR3DL2 > 0) | (KIR3DL3 > 0), slot = "data")

```

Now- I would like to subset if expressed in at least two of these genes- is this possible with the subset function?

# Answer
With a little bit of workaround:

i) Add a new column to the data slot (only because your original subset() call does so but it can be raw counts or any other data matrix in your Seurat object). `rowSums()` determines how many non-zero counts you have.

```
df$non_zero <- rowSums(cbind(df$gene1 > 0, df$gene2 > 0, df$gene3 > 0, ...)) > 1

```

Change df with the your data slot, should be sth like `YourSeuratObject@...$data`.

ii) Subset the newly modified data slot, basically treating the new column as another gene. Be careful not to use this column in preprocessing steps like normalization.

```
YourSeuratObject <- subset(YourSeuratObject, non_zero > 1), slot = "data")

```
> 1 votes

---
Tags: seurat, single-cell
---


title: Querying metadata (GDC) using a filter

text:
# Question
Title: Querying metadata (GDC) using a filter
I have been trying to access the GDC metadata with python through the filter below. However, this code I have been using always yields a blank tsv file. Any help on fixing the code is appreciated!

```
import requests
import json

cases_endpt = 'https://api.gdc.cancer.gov/cases'

# The 'fields' parameter is passed as a comma-separated string of single names.
fields = [
    "sample_id",
    "sample_uuid",
    "sample_type",
    "sample_type_id",
    "tissue_type",
    "tumor_code",
    "tumor_code_id",
    "oct_embedded",
    "shortest_dimension",
    "intermediate dimension",
    "longest dimension",
    "is_ffpe",
    "pathology_report_uuid",
    "tumor_descriptor",
    "current weight",
    "initial weight",
    "composition",
    "time_between_clamping_and_freezing",
    "time_between_excision_and_freezing",
    "days_to_sample_procurement",
    "freezing_method",
    "preservation_method",
    "days_to_collection",
    "portions"
    ]

    fields = ','.join(fields)

filters = {
    "op": "in",
    "content": {
        "field": "project_id",
        "value": ["TCGA-BRCA"]
    }
 }

params = {
    "filters": json.dumps(filters),
    "fields": fields,
    "format": "TSV",
    "size": "10000"
    }

response = requests.get(cases_endpt, params=params)

file = open("query10.tsv", "w")
file.write(response.text)
file.close()

```

# Answer
> 1 votes
I think you are using an outdated version of the API, you need to change the `"field"` with `cases.project.project_id`

```
import requests
import json

cases_endpt = "https://api.gdc.cancer.gov/cases"

fields = [
    "sample_id",
    "sample_uuid",
    "sample_type",
    "sample_type_id",
    "tissue_type",
    "tumor_code",
    "tumor_code_id",
    "oct_embedded",
    "shortest_dimension",
    "intermediate dimension",
    "longest dimension",
    "is_ffpe",
    "pathology_report_uuid",
    "tumor_descriptor",
    "current weight",
    "initial weight",
    "composition",
    "time_between_clamping_and_freezing",
    "time_between_excision_and_freezing",
    "days_to_sample_procurement",
    "freezing_method",
    "preservation_method",
    "days_to_collection",
    "portions",
]

fields = ",".join(fields)

filters = {
    "op": "in",
    "content": {"field": "cases.project.project_id",
                "value": ["TCGA-BRCA"]},
}

params = {"filters": json.dumps(filters),
          "format": "TSV",
          "size": "10000"}

response = requests.get(cases_endpt, params=params)
content = str(response.content, "utf-8")

with open("query10.tsv", "w") as file:
    file.write(content)

```

And if you use *with open*, you don't need to use *close*.

---
Tags: python
---


title: Masking sites in a vcf file

text:
# Question
Title: Masking sites in a vcf file
I need to mask all sites in a vcf file flagged by the 1000 Genomes Project as being unfit for population genetic analyses. The sites for all chromosomes are available at:

1000Genomes masked sites

From above mentioned link you can get both bed files and Fasta files.

I wanted to mask all sites using VCFtools but VCFtools is not longer being updated and therefore, can't me used on the my vcf file because you can only VCFtools upto VCF file format v.2. And file format for my file is v.3. Can someone suggest another tool that I can use. Insights will be appreciated.

# Answer
> 2 votes
I've tried many ways to do this with bcftools but don't think it's possible to exclude a bed file. What you want is

`bedtools intersect -v -a sample.vcf -b mask.bed -wa -header`

The `-v` option gives you only the parts of file a that are *not* in file b. `-wa` only keeps entries from file a, and `-header` preserves the header from file a.

# Answer
> 1 votes
If you are familiar with Python, check out the `pyvcf.VcfFrame.filter_bed()` method I wrote:

```
from fuc import pyvcf
vf = pyvcf.VcfFrame.from_file('your_vcf.vcf')
filtered_vf = vf.filter_bed('your_bed.bed', opposite=True)
filtered_vf.to_file('filtered.vcf')

```

# Answer
> 0 votes
Have a look at BCFtools, which I think does most of the stuff that VCFtools does. I notice there's a `--mask` option in `BCFtools consensus`, which may do what you want:

```
-m, --mask FILE
 BED file or TAB file with regions to be replaced with N
 (the default) or as specified by the next --mask-with
 option. See discussion of --regions-file in Common
 Options for file format details

```

---
Tags: vcf, bedtools, vcftools, bcftools, 1000genomes
---


# Question
Title: How to deal with wildcards (N) in sequence logos
This is an open question that came to my mind recently, though it might be that it is not a common use case. The purpose is to find common patterns in a collection of sequences by representing the patterns as PSSMs. However, some of the bases have sequencing errors and are thus missing. Is there a way of integrating these missing bases in the PSSM calculation?

For a collection of sequences, deriving a sequence logo is relatively straightforward. For the case:

```
ACTAGCGT
ACGTACTC
ACGATCTC
ACGTCAGT

```

Taking the frequencies of each letter at each position would yield the following (pseudocount-corrected) Position Frequency Matrix:

```
      A       G     T       C
0  0.85  0.0375  0.05  0.0625
1  0.05  0.0375  0.05  0.8625
2  0.05  0.6375  0.25  0.0625
3  0.45  0.0375  0.45  0.0625
4  0.25  0.2375  0.25  0.2625
5  0.25  0.0375  0.05  0.6625
6  0.05  0.4375  0.45  0.0625
7  0.05  0.0375  0.45  0.4625

```

And by including background frequencies, we can derive the total Infomation Content for each position, and sum them up to compute the total Infomation Content of the matrix.

But what would happen if the analyzed sequences had experimented errors in the sequencing and some of the letters were missing?

```
ACNNGCGT
ACGTACTC
ACGATNTC
ACNTCAGT

```
* Would you consider `N` a new letter of the alphabet and follow the same procedure? (counting ocurrences per position, calculating background frequencies, include it in the sequence logo...)
* Would you ignore them and normalize each column's frequencies to the number of non-N bases?
* Would you penalize positions/sequences bearing `N`'s when calculating the IC?
* Other?

I wanted to start with the case of `N` since it's the most generic wildcard, and therefore could be easier to address; but if you have ideas on incorporating other more specific wildcards, feel free to add it :)

# Answer
> 3 votes
Since the Position Frequency Matrix is an estimate for the probability to find a particular letter at each position, I think the best way to deal with sequencing errors that result in `N` is to ignore that base for that position and normalize based on the number of non-`N` bases.

The Information Content at each position is just a function of the (estimated) probability distribution. So, once you have the 'best' estimate for the Position Frequency Matrix, I think you should use it directly with no penalty when calculating the IC.

If there is a position with a lot of `N`s, though, you might get a biased estimate of the IC for that position - since using observed frequencies to estimate entropy/information is strongly biased when the number of observations is small.

---
Tags: sequence-analysis, motifs, pssm
---


title: 10X single cell sequencing issue

text:
# Question
Title: 10X single cell sequencing issue
Can you please help me in a technical issue? We sent some samples to single cell sequencing The company says

> Many of the reads were not assigned to cell-associated barcodes. This could be caused by high levels of ambient RNA or by a significant population of cells with a low RNA content, which the algorithm did not call as cells. The latter case can be addressed by inspecting the data to determine the appropriate cell count and using --force-cells.
>
> This would happen if the cells aren’t transcriptionally active – but as these are PBMCs this is unusual .

Have you ever heard about such a problem? Any solution in this case? sending more sample of each patient for sequencing or comparing single cell with matched bulk RNA-seq which one help?

10X say they would do deep sequencing (covering the cost) to get something for us but I don't know why my boss says we should compare single cell with matched bulk RNA-seq.

# Answer
> 1 votes
From what was described, mostly likely they ran cellranger and usually you can get a summary report like this

One of the metric would be `Fraction Reads in Cells` and I think this is what they are referring to this.

My suggestion would be to ask for the report from cellranger and check the metrics and the different plots. If what you have are indeed pbmcs, the reason is most likely not a lot of cells made it into droplets.

Submitting another sample will not help unless you can ensure the cells are of good enough quality, otherwise you will not get much cells back

Sequencing it deeper does not help as well because you are basically sequencing RNA of so called ambient DNA, and it does not give you useful information.

# Answer
> 0 votes
*Answer from @devon-ryan, converted from comment:*

Comparing with matched bulk RNA-seq isn't reasonably related to this, unless you want to ensure the pseudo-bulk is similarish to the actual bulk. Of course they'll be different, but I suppose if you had enough samples and enough background knowledge about what to expect maybe you could get an idea whether the pseudo-bulk way reasonable, which would hit at whether the scRNA-seq was OK. But realistically you'd need a lot of experience with these exact cells in bulk and scRNA-seq to do that. Otherwise, it's a waste of time.

---
Tags: 10x-genomics
---


title: Unifying these files together

text:
# Question
Title: Unifying these files together
I have two .csv files

One for raw read counts of each gene and the other 3114 single cells for four patients which should be the row names of the first file

As single cells per patients seems duplicated I am getting error and don't know how to read them together

Any way this is the files

https://www.dropbox.com/s/slkhspz7lgdgph6/ascites.csv?dl=0

https://www.dropbox.com/s/iq8lt7qrh7hdtq8/ascitesRownames.csv?dl=0

I tried

```
a=read.csv("ascites.csv")
Warning message:
In doTryCatch(return(expr), name, parentenv, handler) :
  restarting interrupted promise evaluation
b=read.csv("ascitesRownames.csv")
View(b)
View(a)
rownames(a)=b

```

> Error in `.rowNamesDF<-`(x, value = value) : duplicate 'row.names' are not allowed In addition: Warning message: non-unique value when setting 'row.names': ‘logical(0)’

Any help in having a matrix in row names are 3114 single cells and genes in columns?

# Answer
> 1 votes
There are a couple of problems:

* As Ram RS said you get a data frame when you load that second CSV. Since the information is stored in a row rather than a column you need a little more work to extract it too.
* R will silently mangle your column names from the first data frame and will also interpret the second file as having column names (and then mangle those too) unless you give extra arguments to `read.csv` to prevent it.
* The fact that you only have four identifiers (not 3114) in the second file means you can't use those as row names. Your rows in `a` are cells, not patients, right? You could make up something unique based on both patient and cell info but I'd suggest avoiding that.

First you could use something like this (going straight from Dropbox here) to read the files without losing information:

```
a <- read.csv("http://dl.dropboxusercontent.com/s/slkhspz7lgdgph6/ascites.csv", check.names = FALSE)
b <- read.csv("http://dl.dropboxusercontent.com/s/iq8lt7qrh7hdtq8/ascitesRownames.csv", header = FALSE)

```

Then, since you want to associate more information than can easily fit in a "wide-format" (matrix-like) data frame, and row and column names are so easily mangled by R, I'd suggest using a "long-format" data frame instead. The package reshape2 is very helpful for this. Here you can put the patient IDs as an additional column in your first data frame, explicitly label the cells, and then make a new data frame where each row is a cell and gene combination for a particular patient:

```
a <- cbind(patient = as.character(b[1, ]), cell_id = 1:nrow(a), a)
gene_counts <- reshape2::melt(
  a, id.vars = 1:2, variable.name = "gene", value.name = "count")

```

So you get something like this:

```
> head(gene_counts)
  patient cell_id gene count
1   7873M       1 A1BG     1
2   7873M       2 A1BG     0
3   7873M       3 A1BG     2
4   7873M       4 A1BG     0
5   7873M       5 A1BG     1
6   7873M       6 A1BG     0

```

(Note the gene column is a factor there.)

I'd be very careful with row and column names since R is so prone to altering them without warning, and in general with assuming things will match up in the right order like the rows of the one table and the columns of the other. Hadley Wickham has a good summary of the problems with row names in particular:

https://adv-r.hadley.nz/vectors-chap.html#rownames

# Answer
> 0 votes
`read.csv` produces a data.frame. You'll need `rownames(a) <- b$col_name`, not `rownames(a) <- b`.

---
Tags: r, single-cell
---


title: Hemoglobin subunits genes in scRNA-seq

text:
# Question
Title: Hemoglobin subunits genes in scRNA-seq
In one scRNA-seq sample I encountered the genes: Hbb-bs, Hba-a1 and Hba-a2. These genes appear on top of the list of the highest expressed genes but the 75% percentile of cells have very low counts. Therefore the expression is highly skewed.

Reading from gene cards, they appear to be coding hemoglobin subunits and involved in oxygen transport. I think they may need to be removed before normalization, also because the analysis I'm running does not consider oxygen transport.

Has anyone encountered this situation and what do I do?

# Answer
You probably have red blood cell contamination. It's probably better to identify which cells they are. They'll probably cluster away from the cells you are interested in.
> 2 votes

# Answer
They may be B cells, did you analysis if they have b cell marker?
> 0 votes

---
Tags: gene, scrnaseq
---


# Question
Title: How to create Beta Diversity distance boxplot
I'm trying to generate boxplot from Beta diversity boxplot between two conditions from a phyloseq object.

I have a Morisita-horn distance generated by vegdist fonction :

```
vegdist<-vegdist(t(otu_table(physeq_genus_3a)), method="horn")#matrice distance Morisita-horn
vegdist<-as.matrix(vegdist)

```

Now I have a matrix and I would like to somewhat create a mean distance boxplot of the Morisita-horn index winthin two condition like the graph below.

I wanted to do a simple boxplot with ggplot but beacause it is a matrix and not a dataframe, I can't do it.

Any help or suggestions would be more than welcome.

# Answer
> 1 votes
1. You can use `as.data.frame()` or `as_tibble()` to convert/coerce your matrix to a data frame or a tibble respectively (Technically a tibble is also a data frame). From there you should be fine with `ggplot()`.
2. I would not use a "dynamite plot" (this is a box plot with error bars like the one in your question), that type of plots conceal the data! This post explains why and remedies.

---
Tags: r, statistics, microbiology, graphs, phyloseq
---


title: Converting gene symbols between species

text:
# Question
Title: Converting gene symbols between species
I have a list of genes in a mouse model

How I can find the equivalent gene symbol/name for human from them?

This is my list

```
CD45
MHC II
CD11b
Ly6C
Ly6G
F4/80
CD11c
CD38
Arg1
SiglecF
CD206
CD62L
CD103
iNOS
PD-L1
TNFa
CD64
TCRgd
Foxp3
RORgt
CD8α
Tbet
CD25
IFN-γ
CD44
CD86
CD80
PD-1
B220
NK1.1
CD19
CD4
TCR β

```

From these I only find two homologous

# Answer
> 3 votes
Easiest way is BioMart. Help video to get you started here. Use the mouse genes dataset and filter by the list of gene names, get the mouse gene name and the human gene name (listed under homologues) as attributes.

---
Tags: gene, sequence-homology
---


# Question
Title: Getting conventional gene symbol for Seurat
I have a Seurat object made by human single cells

I am mapping some genes on that but no sign of expression

When I GOOGLE for those genes I see the genes have different names

How I know Seurat uses which gene name please

This is my list

```
CD45
MHC II
CD11b
Ly6C
Ly6G
F4/80
CD11c
CD38
Arg1
SiglecF
CD206
CD62L
CD103
iNOS
PD-L1
TNFa
CD64
TCRgd
Foxp3
RORgt
CD8α
Tbet
CD25
IFN-γ
CD44
CD86
CD80
PD-1
B220
NK1.1
CD19
CD4
TCR β

```

From this list only 6 genes being mapped on the Seurat data

# Answer
You might want to use Biomart to figure out the official names for those genes based on their aliases.
> 1 votes

---
Tags: scrnaseq, gene, seurat, 10x-genomics
---