# Question
Title: What's the most efficient file format for the storage of DNA sequences?
I'd like to learn which format is most commonly used for storing the full human genome sequence (4 letters without a quality score) and why.

I assume that storing it in plain-text format would be very inefficient. I expect a binary format would be more appropriate (e.g. 2 bits per nucleotide).

Which format is most common in terms of space efficiency?

# Answer
Genomes are commonly stored as either fasta files (.fa) or twoBit (.2bit) files. Fasta files store the entire sequence as text and are thus not particularly compressed.

twoBit files store each nucleotide in two bits and contain additional metadata that indicates where there's regions containing `N` (unknown) bases.

For more information, see the documentation on the twoBit format at the UCSC genome browser.

You can convert between twoBit and fasta format using the faToTwoBit and twoBitToFa utilities.

For the human genome, you can download it in either fasta or twoBit format here: http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/
> 42 votes

# Answer
The standard formats for storing sequence data are fasta and fastq. Fasta is used if you only need the raw sequence data, fastq is used if you want to store the sequence data along with the quality information from base calling. Each of these can be compressed using gzip or another standard compression algorithm.

Typically we want to keep the quality information along with the raw sequence data, but the quality information accounts for half the storage space required. Some people have developed algorithms for lossy compression of the quality data that allow us to reduce the storage requirements.

If you are interested in storing variant calling data, the standard format for that is VCF. VCF is useful if you want to store quality information of the variant calls, genomic positions, and any annotations you might have about the position. VCFs can be compressed and indexed using bgzip and tabix. Many tools require variant data to be compressed and indexed using these tools.
> 25 votes

# Answer
The standard and the most common sequence format is FASTA for sure. You can compress it with a compressor. For the ~3GB human genome, gzip reduces the size to ~900MB, depending on the option in use.

Another often used format is UCSC's 2-bit format. This format keeps each A/C/G/T with 2 bits. As I remember, it keeps non-A/C/G/T bases and lowercases in two separate lists. These lists basically tell you that bases between offset x and y are all "N"/lowercase. The 2-bit format loses IUB codes that GRCh37 has. UCSC's hg19 differs from GRCh37 at a few bases.

BWA also produces its own 2-bit format with indexing. You can generate it separately with:

```
bwa fa2pac -f hg19.fa

```

Unlike UCSC, BWA keeps all IUB codes but loses letter cases. BWA does not provide utilities to convert its 2-bit representation to FASTA, either.

The 2-bit format typically reduces the file size down to 1/4 of its original size, unless there are too many scattered ambiguous bases. For human genome, you get a file ~784MB in size. You can compress it further with gzip, but that actually doesn't work well. A gzip'd 2-bit file is only ~5-10% smaller.

If you want to achieve an even smaller file size, you can compress the BWT of 2-bit file. This gives you a ~633MB file:

```
bwa pac2bwtgen hg19.fa.pac tmp.bwt && gzip tmp.bwt

```

A bit-aware compression algorithm may achieve an even higher compression ratio. However, such BWT-based compression prevents you from extracting subsequences. It is probably of little use in practice.
> 20 votes

# Answer
There are several things to consider when asking for "the most efficient" way to store data, it all depends on your use case. Do you just need ACGT, or are there also IUPAC codings for combinations? Do you need additional data (like quality values)? What kind of application are you using the data for (does it need to load all at once or in in chunks? Once or multiple times? Sequential or random access? etc.pp)?

E.g., most efficient for:

1. Lowest footprint on disk, without a lot of hassle: use either FASTA or 2bit, but run through standard compressor (gzip, bzip2, others). The literature you want to consult here is that of standard text compression I think. Also of interest Large Text Compression Benchmark
2. Keeping the file on disk, but ultra-fast loading small subsets into memory, being able to work in memory with character sized entities: a simple dump of the DNA as characters to disk, maybe combined with an index file to know which chromosome starts where. Then use mmap
3. Storing quality values: See papers like Compression of FASTQ and SAM Format Sequencing Data or Sequence squeeze: an open contest for sequence compression
4. Any combination of the above use cases + a lot more
> 13 votes

# Answer
I think the question is a bit ambiguous so please excuse this answer that's a bit redundant from the rest of the ones provided.

As others have mentioned, if you want to store a full genome, `FASTA` and `2bit` formats are appropriate. For some context, `hg19` is about 900Mb compressed for the `FASTA` file and about 780Mb compressed for the `2bit` file . `hg19` is a reference and is haploid so doesn't represent a "full" human genome that would normally have two alleles for the autosome (non-sex chromosomes).

A common format for representing variant information is Variant Call Format (`VCF`). The `VCF` format represent differences from a reference (`hg19`, say) that can be used to recover the original full sequence by using the reference and the differences encoded in the `VCF` file. I've seen `VCF` files in the range of 100Mb, but a reference file is still needed to recover the full genome sequence which is the range of 800Mb+, as mentioned above.

If you're considering just one "whole genome" in isolation, then the answer is pretty clear: `2bit` format is probably approaching the entropy limit of the human genome and you probably won't be able to do much better. The reason why your question is a bit ambiguous is that as soon as you start encoding more than one genome, a population of genomes, say, then you can start exploiting the redundancy of the genome as shared by the population.

For example, say you want to store two "whole genomes". You could download the `hg19` reference and download two `VCF` files which would give around 1Gb worth of data (around 800Mb for the `2bit` file and around 200Mb for both of the `VCF` files). Now you've been able to represent a "whole genome" in 500Mb instead of the 800Mb. You can see a similar argument for downloading 3 `VCF` files and more.

The minimum amount of information needed to represent a population of genomes is, as far as I know, unknown, but I would guess in the 2.5Mb-5Mb range. For example, see "Human genomes as email attachments" by Christley, Lu, Li and Xie which claims a 4Mb encoding of a genome.

Things get tricky because you have to ask what you're claiming as a "whole genome". `VCF` files are notoriously bad because older versions of the specification only store high quality differences from reference, throwing away high quality called sections. If you want to store low quality information, the encoding is now going to depend on the sequencing technology in weird ways.

Insertions, deletions, mobile insertion elements, copy number variants, other structural variants, etc. all complicate this matter further. Genome Graphs are trying to tackle at least some of these problems but the focus is on variant calling rather than efficient individual whole genome representation, though perhaps can be adapted in the future.
> 8 votes

# Answer
It is not yet standardized, but graph format has the potential for being the most space-efficient method for storing genomes. The idea is this: rather than store a genome as a linear string of sequenced nucleotides, genomes are stored as overlapping graphs, where sequence variants branch off from the reference genome, and then rejoin when the alignment continues. Basically, you start with a reference genome, and for every subsequent genome added to the graph, only the differences are stored. This could allow for an enormous gain in space efficiency.
> 7 votes

# Answer
In terms of raw storage capacity 2 bits per nucleotide, and then further compressed with standard compression techniques would be the most efficient. However, you'd still have other storage considerations. Like what to do about non-standard bases: like if you want to indicate a gap or ambiguity.

I'd also query if it is really necessary to store them as binary since it reduces the readability of the data. It is quite convenient having a whole bunch of unix and programming tools that can operate on the string level in text files.
> 4 votes

# Answer
In all seriousness, the most efficient way to store DNA sequence data is...you guessed it...in DNA. (Church, Gao, and Kasuri, 2012) and others have used DNA synthesis and sequencing as an information writing/reading mechanism.

Practical? Not yet.

Storage efficiency? Unparalleled!
> 3 votes

# Answer
Genozip is lossless compression software that has been optimized for genomic data. If the benchmarks are to be believed, Illumina FASTQ files compressed with Genozip are nearly 5× smaller than FASTQ files compressed with gzip. Equally impressive gains over gzip are seen for FASTA, VCF, and GFF3 files. See the publication in *Bioinformatics* for more details:

**Genozip: a universal extensible genomic data compressor**
> 2 votes

---
Tags: human-genome, storage, file-formats
---


# Question
Title: How to remove duplicates from a fasta file using python
I am using the following command:

```
from Bio import SeqIO
import sys
import re

fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
if len(myfile.seq) > 250:
    gene_id = myfile.id
    mylist = re.match(r"H149xcV_[^\W_]+_[^\W_]+_[^\W_])_[^\W_]+", gene_id)
    print (">"+mylist.group(1))

```

And this is providing me with duplicates of the same gene:

```
>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2
>H149xcV_Fge342_r3_h2
>H149xcV_Fge342_r3_h2

```

How can I reformat my command so that I only receive unique gene id's:

```
>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2

```

# Answer
I wouldn't use Python for this, myself. Instead you can use the FastaToTbl and TblToFasta scripts I have posted before, and pipe to standard \*nix utilities:

```
FastaToTbl file.fa | sort -u | TblToFasta > file.uniq.fa

```

Or, if you don't have GNU `sort`:

```
FastaToTbl file.fa | sort | uniq | TblToFasta > file.uniq.fa

```

Alternatively, you can also do something like this:

```
awk -v RS='>' -v ORS='>' '++a[$0]==1' file.fa | sed 's/>$//'

```

For example, given this input file:

```
$ cat file.fa
>seq1
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
>seq2
ATGGCGCGCTTAGGAGCGCTAGGACT
>seq2
ATGGCGCGCTTAGGAGCGCTAGGACT
>seq1
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
>seq1
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA

```

You can do:

```
$ awk -v RS='>' -v ORS='>' '++a[$0]==1' file.fa | sed 's/>$//'
>seq1
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
ACTTCGCAGAGGCTTCGGAGAGA
>seq2
ATGGCGCGCTTAGGAGCGCTAGGACT

```

Both of the above approaches have the advantage of not needing to store the entire file in memory.
> 2 votes

# Answer
```
from Bio import SeqIO
import sys
import re

unique = []
fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
    if len(myfile) > 250:
        gene_id = myfile.id
        if gene.id not in unique:
            unique.append(gene.id)
            mylist = re.match(r"H149xcV_[^\W_]+_[^\W_]+_[^\W_])_[^\W_]+", gene_id)
            print (">"+mylist.group(1))

```
> 1 votes

# Answer
Perhaps itertools might be an option? From the itertools page / moreitertools project https://docs.python.org/3/library/itertools.html#itertools-recipes

```
def unique_everseen(iterable, key=None):
    "List unique elements, preserving order. Remember all elements ever seen."
    # unique_everseen('AAAABBBCCDAABBB') --> A B C D
    # unique_everseen('ABBCcAD', str.lower) --> A B C D
    seen = set()
    seen_add = seen.add
    if key is None:
        for element in filterfalse(seen.__contains__, iterable):
            seen_add(element)
            yield element
    else:
        for element in iterable:
            k = key(element)
            if k not in seen:
                seen_add(k)
                yield element

```
> 1 votes

# Answer
I would use a `dict` so that the lookup time is constant. If you are creating an array with all the unique gene IDs, then checking if a gene is in that array, your runtime is O(N2). When you create an entry in the `dict`, the time to access is constant, and if the entry is already present, it will just overwrite it. Something like this should do it for you:

```
from Bio import SeqIO

genes = {}
fasta_file = sys.argv[1]
for entry in SeqIO.parse( fasta_file, "fasta" ):
    genes[ entry.id ] = entry.seq

```

Then you can write a new, unique fasta file:

```
with open('unique_genes.fasta', 'w') as out:
    for key,value in genes.items():
        out.write( '>%s\n%s\n'%(key,value) )

```
> 1 votes

# Answer
Your question isn't well formatted, for example list is used as a variable name where it should be reserved as a keyword...

```
print (">"+list.group(1))

```

However, the answer you're after is going to be something like:

```
print (">"+set(name_of_list.group(1)))

```

Anytime you want to get the unique elements of a list, just wrap it in set()
> 0 votes

# Answer
I use the following code. It doesn't only remove repeated sequences but similar sequences (with a quite simplistic approach). I use it just sporadically and I guess it can be largely depurated. Remove the debug prints if these annoy you.

```
import sys
import requests
import shlex, subprocess
import random
from Bio import SeqIO
import editdistance

if len(sys.argv)<3:
    print("Uso: "+sys.argv[0]+" input_fasta output_fasta")
    sys.exit(0)

savefile = sys.argv[2]
content_file = sys.argv[1]

class Sequence:
    id=""
    sequence=""

def are_similar(seq1, seq2):
    l1 = len(seq1.sequence)
    l2 = len(seq2.sequence)
    if abs(l1-l2) < l1/40:
        d = editdistance.eval(seq1.sequence[0:250], seq2.sequence[0:250])
        #d = editdistance.eval(seq1.sequence, seq2.sequence)
        print(seq1.id+" "+seq2.id+" -> "+str(d))
        #if d < l1/20:
        if d < 50:
            return True
    return False

def find_repeats(sequences, seq, repeated):
    b=False
    for seq2 in sequences:
        if(b and not seq2.id in repeated):
            if are_similar(seq, seq2):
                repeated.add(seq2.id)
        if seq2.id==seq.id:
            b=True

sequences=[]

res=""

for record in SeqIO.parse(content_file, "fasta"):
    seq=Sequence()
    seq.id=record.id
    seq.sequence=record.seq
    sequences.append(seq)

repeated_sequences=set()
for seq in sequences:
    find_repeats(sequences, seq, repeated_sequences)

print(repeated_sequences)

def mySortFunc(e):
    return len(e.sequence)

sequences.sort(key=mySortFunc)

for seq in sequences:
    b= seq.id in repeated_sequences
    if not b:
        res+=">"+seq.id+"\n"
        res+=str(seq.sequence)+"\n"

output_file = open(savefile, "w")

output_file.write(res)

output_file.close()

```

EDIT: I see that's not exactly what OP is asking ...
> 0 votes

---
Tags: fasta, python
---


# Question
Title: Where to get BED files for regions that contain binding site motifs for specific transcription factors?
Quick question, is there a place where one can download BED files for the binding motifs of certain transcription factors. In my case, I'm looking for a BED file for the genomic regions (enhancers) where TET2 binds to genome-wide. I'm looking to examine the DNA methylation levels in those regions.

Thank You!

# Answer
The standard answer would be to search JASPAR, to get a curated dataset. However, there does not appear to be anything coming up for TET2 in the database.

If your TET2 gene is what you want specifically, you can find papers that perform ChIP-seq and infer binding regions. I believe that this Supp Table 3 contains information equivalent to a BED file for TET2.
> 1 votes

---
Tags: phylogenetics, bed, motifs, transcription-factors
---


# Question
Title: What is the best way to compare treated/control samples?
I have two populations of mice, say RED and GREEN; there are 5 pairs treatment/control in each.

The table is the following:

```
gene_name  p-value(treatVSctrl) RED1_ctrl RED1_treat RED2_ctrl RED2_treat ... GRN1_ctrl GRN1_treat etc...
gene1      1.3                  4         10         5         10             6.2       9.75
gene2      2.8                  3         6          6         12             3         6
gene3      4.4                  5         12         4         8              7         10

```

I need to compare two groups to have an idea of the expression of different genes after treatment. I also need to find a difference between RED and GREEN group in terms of genes expression

The thing is I cannot compare mice one by one just measuring the fold change, I need to consider p-value. What is the most correct way to perform this analysis?

# Answer
> 3 votes
I don't really understand the table, particularly the 'p-value' in the header. More generally I don't really understand the experiment nor the data.

As a general rule the Wilcoxon would be a reasonable guess of the test using raw data. This is non-parametric version of the T-test without the rigid assumptions of parametrics. Control mice go in one column and subjects in the other. Scipy.stats has the test, its certainly in R (which I don't really use).

```
print (scipy.stats.wilcoxon(x = [1,2,3,4,5], y=[1,3,5,6,7], zero_method='wilcox', correction=False, alternative='two-sided', mode='auto', *, axis=0, nan_policy='propagate'))

```

There you go, filled out with some example data.

---
Tags: differential-expression, p-values
---


# Question
Title: Why is bulk RNA sequencing reflecting AVERAGE expression but not TOTAL expression of all cells?
When I am reading papers that compares bulk RNA sequencing and single-cell RNA sequencing, we often see papers describe bulk RNA seq measures the average cell expression.

For example, in this paper Single-Cell RNA-Seq Technologies and Related Computational Data Analysis

> bulk RNA-seq mainly reflects the averaged gene expression across thousands of cells

And this paper A practical guide to single-cell RNA-sequencing for biomedical research and clinical applications

> the averaging that occurs in pooling large numbers of cells does not allow detailed assessment of the fundamental biological unit—the cell—or the individual nuclei that package the genome

And this paper Differential expression analyses for single-cell RNA-Seq: old questions on new data

> Standard RNA-Seq experiments need millions of cells for sequencing \[2,3\], and therefore can only get averaged measurements of gene expressions of the cells sequenced.

My question is, when we perform the bulk RNA seq and calculates parameters like TPM, we don't actually divide the total number of mRNA transcripts with the number of cells (I think sometimes we don't even have an accurate number on the number of cells that we used in bulk RNA seq), we simply normalise it with the transcript length and count.

Therefore, if we **DID NOT** divide the number of transcripts with the total cell number, how are we measuring the average of expression, but not the total expression?

Or in other words, what is the "average expression" referring to? Average to what?

# Answer
What these statements refer to is that bulk RNA-seq will overlook the different transcriptional cell types, hence the value of single cell RNA-seq in unraveling heterogeneous cell populations/tissues.

I think "total" is also avoided because total RNA has another meaning related to capturing all the RNA transcripts in a cell (e.g. including non polyA transcripts), which is still not the case for most RNA-seq experiments.

Also note that RNA-seq experiments generally (without spikins) are only *relative* and not giving *absolute* quantifications of the amount of RNA in a cell or groups of cells.
> 2 votes

---
Tags: rna-seq, scrnaseq, transcriptome, rna
---


# Question
Title: GATK variant filtration using "SelectVariants" and use of JEXL queries
How to do variants selection in some corner cases using GATK and JEXL expressions?

I am following the guidelines given in this links for variant selection for some specific cases are not as much detailed.

I just had to do some filtering of my own data and tested several methods (to see what the outcome would be). I am sharing those here for others and myself, if it is ever needed again.

# Answer
I am following the guidelines given in this links for variant selection for some specific cases that are not as much detailed in the GATK site.

I just had to do some filtering of my own data and tested several methods (to see what the outcome would be). I am sharing those here for others and myself, if it is ever needed again.

## **Tutorials .... !!!**

**\- method to call the sites that are all homozygous Variants (1/1, 2/2 genotype)**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomVarCount() == 6'
# outcome: no sites will have ambigous (./.) and/or homRef (0/0). All other homVar are allowed

```

**\- in the above script if we want the site where at least 4 samples are homozygous ref or variants**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomVarCount() >= 4'
# outcome: atleast 4 samples are homVar, rest can be any genotype (./., 0/0, 0/1, etc.)

```

**\- at least 1 sample is homVar (1/1, 2/2 or 3/3)**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomVarCount() != 0'

```

**\- all the samples are homRef (0/0 genotypes)**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomRefCount() == 6'

```

## Note: 6 chr = 3 samples (when samples are diploid)

**\- sites that are all HomVar or nocalls (./.)**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomVarCount() == vc.getCalledChrCount()/2'
# outcome: no 0/0, 0/1
# Note: samples with GT = ./. are not counted as called chr

```

## Useful after samples are sub-sampled from multisample vcfs

**\- sites where no genotypes were called for all the samples in the selected vcf** **\- or say, select vcf sites that are not called in any sample (i.e don't have any genotypes, only contain ./.)**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getCalledChrCount() == 0'
# outcome = vcf sites (lines) that are all no calls (GT = ./.) on all samples are selected

```

**\- samples are either homRef (0/0) or no call (./.)**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomRefCount() == vc.getCalledChrCount()/2'
# outcome: vcf sites where all the samples are GT = ./. (no call) or GT = 0/0

```

**\- site where 6 chromosomes (exact) or 3 samples were called for diploid genome**

```
java -jar /home/everestial007/GenomeAnalysisTK-3.7/GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getCalledChrCount() == 6'
# Note of Caution: 6 chr = 3 samples (when samples are diploid)

```

**\- site where at least 3 samples (diploid) are called**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getCalledChrCount() >=6'

```

**\- if I have 6-sample vcf**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getCalledChrCount() ==12'
# 12 haploid chr = 6 samples (diploid)
# outcome: sites that were called in all samples, there will be no GT = ./.

```

**\- select vcf sites where no sample has hetVar (no GT = 0/1, 0/2), only GT = ./., 0/0, 1/1, 2/2**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHetCount() == 0'

```

**\- select vcf sites where no sample has hetVar (no Gt = 0/1, 0/2), and at least one sample is HomVar (GT = 1/1)**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHetCount() == 0 && vc.getHomVarCount() >= 1'

```

**\- site where at least 1 sample is homVar (GT = 1/1 and/or 2/2 etc.)**

```
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomVarCount() >= 1'

```

**\- where a site is only homRef for a particular sample and is only SNP**

```
java -jar GenomeAnalysisTK.jar SelectVariants -R lyrata_genome.fa -V phased.MySpF1.vcf -sn ms01e -select-type SNP -select 'vc.getGenotype("ms01e").isHomRef()' -O test_variants.ms01e.Chr_2_3.vcf

```

**\- where a site is only homVar for a particular sample and is only SNP**

```
java -jar GenomeAnalysisTK.jar SelectVariants -R lyrata_genome.fa -V phased.MySpF1.vcf -sn ms01e -select-type SNP -select 'vc.getGenotype("ms01e").isHomVar()' -O test_variants.ms01e.Chr_2_3.vcf

```

**\- where a site is not homRef for a particular sample and is only SNP. This will also select empty sites**

```
java -jar GenomeAnalysisTK.jar SelectVariants -R lyrata_genome.fa -V phased.MySpF1.vcf -sn ms01e -select-type SNP -select '! vc.getGenotype("ms01e").isHomRef()' -O test_variants.ms01e.Chr_2_3.vcf

```

**NOTE: following genotype call methods won't work:**

```
# to get only heterozygous sites
-select 'vc.getGenotype("ms01e").isHetVar()

# to get any variant sites
-select 'vc.getGenotype("ms01e").isVar()

# won't work if getGenotype method is combined
-select 'vc.getGenotype("ms01e").isHomVar() || vc.getGenotype("ms01e".isHomRef()'

```

**So, I figured out the following methods where count helps with selecting particular site**

**\- where a site is a heterozygous variant site for a particular sample**

```
java -jar GenomeAnalysisTK.jar SelectVariants -R lyrata_genome.fa -V phased.MySpF1.vcf -sn ms01e -select 'vc.getHetCount() == 1'  -O test_variants.ms01e.vcf

NOTE: 'vc.getGenotype("ms01e").isHetVar()' wont work so generated the above expression.

```

**\- where a site is a variant site (either homVar or hetVar) for a particular sample and is only SNP**

```
java -jar GenomeAnalysisTK.jar SelectVariants -R lyrata_genome.fa -V phased.MySpF1.vcf -sn ms01e -select-type SNP -select 'vc.getHetCount() == 1 || vc.getHomVarCount() == 1'  -O test_variants.ms01e.vcf

NOTE: 'vc.getGenotype("NA12878").isHomRef()' || vc.getGenotype("NA12878").isHomVar()' wont work together, so generated this above expression.

```

Other updates to follow. Thanks,
> 1 votes

---
Tags: variant-calling, gatk, variant-filtration
---


# Question
Title: Retrieval of metadata of a batch of 3953 sequence data from GISAID database
I need to retrieve COVID-19 metadata (such as date of collection of sequence data, location, etc.) for a batch of 3953 GISAID accessions from the GISAID database.

I am not very comfortable with GISAID. I would prefer an NCBI database solution.

Any help would be appreciated.

# Answer
You can download the full db and and put out anything you want with your favourite scripting language.

To do so, you need to request an account first. Once you log in, scroll down to the bottom of the download options. There are 2 meta-data download options, pick what you need.
> 1 votes

# Answer
The NCBI link to the metadata database is here

If I get a moment I'll talk you through it.
> 0 votes

---
Tags: sequence-analysis, public-databases, phylogenetics, covid-19, sars-cov-2
---


# Question
Title: How can I list all amino acids in a protein using PyMOL?
What is the command in PyMOL for listing all the amino acids in a specific protein, say, 1a62.pdb?

# Answer
> 2 votes
At the command, change the current directory to your target directory:

```
cd my-directory-path-containing-my-files

```

Load the file:

```
load 1a62.pdb

```

Then you can either print to screen or save to a file:

```
print(cmd.get_fastastr('all'))

save 1a62.fasta

```

---
Tags: pdb, pymol
---


title: What is the tP statistic from ANGSD?

text:
# Question
Title: What is the tP statistic from ANGSD?
We were looking for genetic diversity in the exon level of an MHC transcript and analysing the data under ANGSD (Analysis of next generation Sequencing Data) written in C++

We used ANGSD to calculate Thetas and Tajima D neutrality tests and one of the thetas was a pairwise theta, however I am unsure what this means and how it differs from Tajima's D?

http://www.popgen.dk/angsd/index.php/Thetas,Tajima,Neutrality\_tests

We then divided the tP over length of each exon to find genetic diversity in the exon level. However, as I am very interested in statistics, if someone could explain what the tP statistic means exactly it would be highly appreciated.

# Answer
While I can't explain it to you directly, you might check out the section on theta estimators in Korneliussen et al. (2013) which is cited on the above ANGSD wiki page. The equation for this estimator in this paper traces back to the original Tajima (1989).
> 1 votes

# Answer
I don't know Theta but Tajima D is a population genetics test that assumes constant sample size. When the number of segregating sites (non-conservative sites in an alignment) significantly differs from pi (Nei's p-distance between a pair of sequence) one possibility is selection given a constant population size.

What Tajima found was cool: under neutrality this relationship conforms to a beta distribution. Just to explain that again so the nucleotide differences between every pair of sequences in an alignment (pi) against the total number of possible mutations conforms to a beta distribution for a constant population size under neutrality.

Its a cool idea - like really cool - but the violations - any level of 'extinction' will skew the distribution, population movement and importantly so could a bias in sampling.

So what you get is a pair-wise matrix at 5% critical probability. If your data entirely conformed to Tajima's D it would be an amazing result, but that is not the result you are looking for. So if Tajima D gave a particular pair \<0.05 AND it was supported by other tests THEN (sounds like code) it provides a basis to consider selection as an explanation.

Again I don't know Theta but if both tests were congruent for putative selection I'd say that was reasonable. However, Tajima D alone is not a cool result, its cool theory however but not cool in practice. I hope that explains it?

Note: did any one ever follow this up in trees? Cool if they did because there will be a relationship between the balance of a tree (e.g. Bayesian 'Skyline' plots) and Tajima D.
> 0 votes

---
Tags: statistics, phylogeny, phylogenetics, exon
---


# Question
Title: Nextflow scheduler job doesn't end upon successful workflow completion
## Bug fixed

This issue appears to have been fixed in latest `Nextflow`: https://github.com/nextflow-io/nextflow/issues/2693

---

I've been using `Nextflow` for various workflows on a cluster using the `SLURM` scheduler. I've noticed that the 'parent' job does not shut down when the `Nextflow` workflow finishes successfully. I can cancel the job, but then `Nextflow` doesn't get a chance to clean things up and also doesn't write out the `trace`, `report`, etc. The parent job does end properly if the workflow ends with an error.

**Question**: Why doesn't the parent `Nextflow` job end upon successful completion?

I saw a post on the `Nextflow` GitHub where the original poster describes various acrobatics to get the parent job to exit cleanly, but this is something that `Nextflow` should handle.

Here's my `SLURM` submission script that I submit using `sbatch run_workflow.slurm`

**run\_workflow.slurm**

```
#!/bin/bash
#SBATCH --time=60:15:00                             # Time limit for the job (REQUIRED).
#SBATCH --job-name=NF_Parent_Identify_D&C_regions   # Job name
#SBATCH --ntasks=1                                  # Number of cores for the job. Same as SBATCH
#SBATCH --mem=5G                                    # Total memory requested
#SBATCH --partition=<partition>                     # Partition/queue to run the job in. (REQUIRED
#SBATCH -e slurm-%j.err                             # Error file for this job.
#SBATCH -o slurm-%j.out                             # Output file for this job.
#SBATCH -A <account>                                # Project allocation account name (REQUIRED)

module load ccs/java/jdk1.8.0_202
nextflow run identify_dark_and_camouflaged_regions.nf \
    -with-report \
    -with-trace \
    -with-dag idcr-workflow.png \
    -with-timeline \
    -resume

```

# Update 1

I'm adding example `Nextflow` output from a different workflow (I don't have the `SLURM` output from the one I mention above) where the workflow has completed, but the parent job is still running. In this case, there were errors, but I added `errorStrategy 'ignore'` to all of the processes listed:

```
executor >  slurm (876)
[06/c5672b] process > samtools_csort_proc (A-CUHS... [100%] 304 of 304, faile...
[45/26a046] process > Lsamtools_csort_proc (A-CUH... [100%] 304 of 304, faile...
[b5/0f2d7c] process > sambamba_csort_proc (A-CUHS... [100%] 304 of 304, faile...

```

So, regardless of whether the workflow completed `successfully`, `Nextflow` has nothing left to do, but it's still running.

# Update 2 (adding partial `.nextflow.log` file)

Per @Steve's point, I checked the `.nextflow.log` file and am including part of the file, here. In this case, I don't think `Nextflow` has any tasks related to `publishDir` because I'm not using that in this workflow. The job that's still running is from the `Nextflow` workflow in my GitHub repository for my answer posted in my related question Samtools sort: most efficient memory and thread settings for many samples on a cluster. Whether it's reproducible, I haven't checked. Not sure how common this issue is.

In either case, my interpretation of the `.nextflow.log` file below is that `Nextflow` thinks it's hung, too. I just don't know why, or how to exit cleanly.

**Partial view of `.nextflow.log`**

```
...

Feb-28 05:29:04.607 [Task monitor] DEBUG n.processor.TaskPollingMonitor - !! executor slurm > tasks to be completed: 1 -- submitted tasks are shown below
~> TaskHandler[jobId: 662274; id: 876; name: sambamba_csort_proc (A-CUHS-CU003128-BL-COL-49696BL1.unsorted;CPUs:2;MEM_PER_THREAD:100GB); status: RUNNING; exit: -; error: -; workDir: /path/to/Samtools_sort_optimization_test/SORT_COMPARISON/work/b5/0f2d7c5d7f81a6e62dedb08bf6694a started: 1646038174334; exited: -; ]
Feb-28 05:30:54.613 [Task monitor] DEBUG n.processor.TaskPollingMonitor - Task completed > TaskHandler[jobId: 662274; id: 876; name: sambamba_csort_proc (A-CUHS-CU003128-BL-COL-49696BL1.unsorted;CPUs:2;MEM_PER_THREAD:100GB); status: COMPLETED; exit: 0; error: -; workDir: /path/to/Samtools_sort_optimization_test/SORT_COMPARISON/work/b5/0f2d7c5d7f81a6e62dedb08bf6694a started: 1646038174334; exited: 2022-02-28T10:30:53.516015Z; ]
Feb-28 05:34:04.641 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:39:04.647 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:44:04.653 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:49:04.658 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:54:04.662 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled

< above line repeats hundreds of times >

...

```

# Answer
How do you know your workflow actually finishes successfully? One way to check the workflow termination status is to use the `workflow.onComplete` and `workflow.onError` event handlers. Consider resuming your workflow with something like the following appended to your script:

```
workflow.onComplete {

  log.info "Workflow completed at: ${workflow.complete}"
  log.info "Time taken: ${workflow.duration}"
  log.info "Execution status: ${workflow.success ? 'success' : 'failed'}"
}

workflow.onError {

  log.info "Execution halted: ${workflow.errorMessage}"
}

```

Note that the workflow event handlers may also be defined inside your `nextflow.config` if preferred.
> 2 votes

---
Tags: nextflow, slurm
---


# Question
Title: Replacing column 1 in file1 with matching ID's in file 2
I have downloaded two files:

**File1**.csv from: https://depmap.org/portal/download/data\_slicer/download?file\_path=20220227%2F3c335b7d-c18e-4788-b2e6-512d71f0a339%2Fexport.csv&name=Metabolomics.csv

**File 1** looks like

```
,2-aminoadipate,3-phosphoglycerate,alpha-glycerophosphate
ACH-000698,6.1127272,6.034197900000001,5.8968963
ACH-000489,5.5774126,5.7270453,5.1114679,6.07325

```

**File2**.csv from: enter link description here https://cog.sanger.ac.uk/cmp/download/model\_list\_20220205.csv

File2 looks like:

```
 model_id,model_name,synonyms,model_type
 SIDM01774,PK-59,,Cell Line,Adherent
 SIDM00192,SNU-1033,,Cell Line,Adherent

```

I want to replace ACH-\* ids with SIDM\* ids in File1.csv. The ACH-\* id's for SIDM\* ids are given in the "BROAD\_ID" (col.no =40) column in the File2.csv

Desired output

```
,2-aminoadipate,3-phosphoglycerate,alpha-glycerophosphate
SIDM00523,6.1127272,6.034197900000001,5.8968963
SIDM00835,5.5774126,5.7270453,5.1114679,6.07325

```

I want to do this with python script. I tried some initial steps. Please forgive any grave mistake, as I am a beginner in Python

```
 import pandas as pd

 df = pd.read_csv("File1.csv", index_col=0).T

 print(df)

 models = pd.read_csv("File2.csv", index_col="model_id", usecols=["model_id", "BROAD_ID"]).rename(columns={"model_id": "Model"})

 df = df.join(models).rename_axis(index="Model Name").reset_index()

 print(df)

 df2 = df.set_index("Model").drop(columns="Model Name")

 df2.to_csv("Meta.tsv", sep="\t")

 df2.index = df2.index + ".Demeter"
 df2.T.rename_axis(index="Gene").to_csv("Meta_2.tsv", sep="\t")

```

# Answer
I was not sure about the exact output (columns) you desired. The following code should inner merge both data frames. It writes a copy of 'model\_list' column as a first column. You could remove any unwanted columns. I am sure there are more elegant ways of doing it. Hope it helps. It's totally different approach than yours.

```
import pandas as pd

metabolics = pd.read_csv("Metabolomics.csv", sep=",", header=0)
model_list = pd.read_csv("model_list_20220205.csv", sep=",", header=0)
metabolics=metabolics.rename(columns={'Unnamed: 0':'ACH_number'}) #rename the unnamed column.
merged_df = pd.merge(metabolics, model_list, left_on='ACH_number', right_on='BROAD_ID', how='inner') #join two files

#insert model_id column at the beginning
#You will have dupicae 'model_id' column
merged_df.insert(0, "model_id", merged_df['model_id'], allow_duplicates=True)
merged_df.to_csv('results.csv', sep=",", index=False)
print(merged_df.shape)

```
> 4 votes

# Answer
The basic strategy of `merge` 'left' (which I think is default) is the correct thing to do. However, I would be concerned by lack of a key between the two data sets and basically the merging will take place on the row order. Within `merge` is `on=` then the shared `key`, in the code by @supertech there is left-on='ACH\_number', I dunno because `ACH_number` does not appear to be shared. The code needs testing.

If you assign an index it will attempt to merge on the index (which I think will therefore fail). What I think you want is the row order 1,2,3,4 in each dataframe to be the key that is merged. Therefore not specifying `key` might be the right thing to do, what can happen is it will then use the row order as the left-key. I cannot remember.

The `insert` is correct to reorder the columns to get the desired output.

Finally you want,

```
del merged_df['column_name']

```

for each column you don't want. There is a way to loop through this, but the above keeps things simple. You then do the standard csv dump.

Note, see comments on `drop` by @GOATNine, which strongly preferred.

Setting and unsetting indexes is key to manipulating pandas but it is within a specific contex.

**Summary** the command `merge` and `insert` are correct with the addition of `del` however, the code needs testing regarding merge `left-on` and assigning the index. It is much easier to print the output and visualise it.
> 2 votes

---
Tags: python, data-download, data-preprocessing
---


title: SNPs in LD in which populations?

text:
# Question
Title: SNPs in LD in which populations?
I have 3 SNPs that are in LD in 'several' populations here: https://www.ebi.ac.uk/gwas/variants/rs13078854#ld\_panel\_label

It's a pain to look at all (30) populations one by one to manually check if the 3 SNPs are in LD in a given population or not.

Is there a way to ask this question more systematically? e.g. in which common imputation populations are a given set of SNPs in LD (above a certain threshold)?

# Answer
There is a risk of a type 2 error.

My advice would be to use `Structure` if it is still available given the situation. Its not great and it is old and its a GUI and I don't know whether it has been updated (I'm not exactly selling this am I?), but it will do the job of navigating through a type 2 error.

What you do is provide the number of populations *a priori* and it will assign a Bayesian estimate of HW equilibrium for that number of populations. You manually input 1,2,3,4,.... 30 run and wait. The point at which the joint marginal plateaus is the minimum number of populations which do not exhibit LD (thus if two of your populations exhibit LD it collapse them into one). Thus if your answer is '20' it give some confidence it is a real result. If it is like 10 (optimistically) the conclusion would be in pretty reasonable shape.

There are other ways of looking at this, but that would be my starting point.

The classic calculation is LDHat, here https://github.com/auton1/LDhat, but the last update was 2015 and in this particular circumstance I'm not sure.
> 2 votes

---
Tags: snp, imputation, ld
---


title: Gold standard benchmark

text:
# Question
Title: Gold standard benchmark
This page is claimed to contain a gold standard benchmark for viral genome assembly. https://github.com/cbg-ethz/5-virus-mix

The claim is here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411778/

> Baaijens *et al*. De novo assembly of viral quasispecies using overlap graphs Genome Research 2017 May; 27: 835–848.

"*This data set was recently presented as a gold standard benchmark (Di Giallonardo et al. 2014)*"

But at (Di Giallonardo et al. 2014) I can only find that mention:

"We previously estimated the frequencies of the virus strains in the 5-virus-mix by amplifying the protease gene using single-genome amplification (SGA), the current gold standard for studying diversity of virus populations".

That doesn't seem to mean precisely that the dataset is a gold standard benchmark for genome assembly but, as I understand, that SGA is a gold standard for something else.

My question is: is it? And who decided it is? Is still valid? Is there an official list of results somewhere for this benchmark set?

I tried to contact the owner of the page but still no answer.

# Answer
Anyone can declare anything to be a "gold-standard", the term is meaningless. All if actually means is, "we think this is the test that everything else should be benchmarked against".
> 1 votes

# Answer
Thanks for raising the issue, it is certainly topical. The question is whether the data set they unit tested is reliable.

I think the issue they are trying get to is mosaic artefacts in the sequencing process and assembly thereon. This would be particularly true in a situation of *de novo* assembly which is what the authors are targeting. There are some slick methodologies for virus sequencing, which have been around for a while now, so in any case a 'gold-standard' would be outdated at present.

I do have to agree with you that the language particularly for the journal certainly could have phrased better. The data set was 2014 and genome circularisation (its a virus small genomes) to remove the otherwise high SNP error rate would have been around in 2017 when this was published.

I fully understand why the authors choose this language. A professional editor, not affiliated with the journal (without up to date technical insight), would likely recommend the same language.

I guess the question would be whether verification is achieved on other virus data sets. There is no shortage of them at present and its certainly an interesting tool.
> 1 votes

---
Tags: genome, assembly, covid-19, benchmarking, virus
---


title: samtools view command not found error

text:
# Question
Title: samtools view command not found error
When I tried to use samtools to split a bam file based on different chromosomes, I used this command:

```
samtools view input.bam -b chr21 | chr21.bam

```

However, I get error messages like this:

```
-bash: chr21.bam: command not found
[W::hts_idx_load3] The index file is older than the data file: input.bam.bai

```

How could I solve this problem? Thanks in advance~~~

The version I am using is:

```
Program: samtools (Tools for alignments in the SAM format)
Version: 1.10 (using htslib 1.10.2)

```

# Answer
> 5 votes
Try re-indexing your bam file:

```
samtools index input.bam

```

Or update the timestamp via `touch index.bam.bai` but I think re-indexing is safer. You might want to save any original index file first, e.g. via `mv input.bam.bai old_input.bam.bai` or similar

Also I think you meant `> chr21.bam` instead of piping via `| chr21.bam` if you want to save the output file

---
Tags: sam, samtools
---


# Question
Title: bcftools mpileup "Failed to open the file" questions
When I trying to use bcftools mpileup to convert bam to vcf files, I try to use:

```
bcftools mpileup -Ob -o resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta.fai input.bam

```

However, I got this error message:

```
[E::fai_build3_core] Failed to open the file resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta.fai

```

Or, I would get another error message when I use another ref file:

```
[E::fai_build_core] Format error, unexpected "c" at line 1

```

Does anyone know how to solve these problems? Thanks in advance

Also, could anyone recommend some human gene Reference? The one I am using now is from GATK pipeline (but I don't know if that is good enough).

The version I am using is:

```
Program: bcftools (Tools for variant calling and manipulating VCFs and BCFs)
License: GNU GPLv3+, due to use of the GNU Scientific Library
Version: 1.10.2 (using htslib 1.10.2)

```

# Answer
Regarding a human genome reference, have a look at CHM13/HG002 from the Telomere-to-Telomere consortium:

https://github.com/marbl/CHM13#assembly-releases

The benefit of the T2T assembly is that it is fully-assembled across all repetitive regions, which improves the mapping accuracy. The pangenome consortium is working on additional T2T assemblies from other people, trying to capture as much human \[structural\] variation as possible.

https://humanpangenome.org/

The v2.0 assembly, adding chrY from HG002, has only just been released. From the looks of it, it is not yet annotated, but there's gene annotation for all other chromosomes via the v1.1 release:

https://github.com/marbl/CHM13#downloads
> 2 votes

---
Tags: bcftools
---


# Question
Title: Identify side chain atoms in BioPandas dataframe
I am currently implementing some metrics I could use for comparing two conformations of the same protein in Python.

For example, I know I could use the RMSD of all protein ATOMS in BioPandas using this:

```
biopandas.pdb.PandasPdb.rmsd(conformation1.df['ATOM'], conformation2.df['ATOM'])

```

Are there some out-of-the-box tools in Python for measuring the RMSD between the side chains of the two conformations? When trying to do this manually, the main issue is how to identify side chain atoms, which is not clear for me. I know how to get all atoms from a residue using BioPandas:

```
aminoacids = list(PandasPdb().read_pdb(my_pdb_filename).df['ATOM'].groupby(['chain_id', 'residue_number', 'residue_name']))

```

Now `aminoacids[i]` is a dataframe containing all the atoms of the i-th aminoacid, like this:

From here, it's not clear how to separate back bone atoms from side chain atoms. It's clear that the one with atom\_number 'CA' will be the alpha-Carbon, but I am not sure how to identify the others. But once the separation is done, the side chain RMSD could be measure for example using `Bio.SVDSuperimposer.SVDSuperimposer` which works directly on the array of 3D coordinates of the atoms.

Thanks!

# Answer
As you suggest one way of solving your problem would be by selecting all atoms that don't have backbone atoms names. In a pdb file I believe backbone atoms would be named 'CA', 'HA', 'N', 'HN' or 'H', 'C' and 'O'. Beware of the N-terminal (where the hydrogens would be named either 'H1', 'H2' & 'H3'; or 'HN1', 'HN2' and 'HN3') and the C-terminal (the other oxygen would be named 'OXT').

One way of doing it is by using the handy selections from pandas dataframes. It's probably a more elegant way than making a list of dataframes.

```
from biopandas.pdb import PandasPdb

bpdf_1 = PandasPdb().fetch_pdb('1t48').df['ATOM']
bpdf_2 = PandasPdb().fetch_pdb('1t49').df['ATOM']

# Function returns a list of boolean. True when it's not a backbone atom.
sel = lambda df: df['atom_name'].str.contains('[^(CA|HA|N|C|O|HN|H)]$')

# The number of atom is not the same, so I only align the first 200
sidechain_rmsd = PandasPdb.rmsd(bpdf_1[sel][:200], bpdf_2[sel][:200])

```

The regex says something like "(not (CA or HA or …)) end of line". The end of line is important, otherwise any atom containing N, C, O or H in their name would be excluded (so most of them).

Note that with this method it's still possible to use the PandasPdb selectors. For example, to only consider heavy atoms (no hydrogen):

```
sidechain_rmsd = PandasPdb.rmsd(bpdf_1[sel][:200],
                                bpdf_2[sel][:200],
                                s='heavy')

```

**EDIT:**

Another way to do it would be:

```
sidechain_rmsd = PandasPdb.rmsd(bpdf_1[:200],
                                bpdf_2[:200],
                                s='main chain',
                                invert=True)

```

You don't need any selection in this case, **BUT** this is a bad idea if you have hydrogens, as it will also align them (including the backbone ones).
> 3 votes

---
Tags: python, proteins, protein-structure
---


title: Selecting a part of heatmap

text:
# Question
Title: Selecting a part of heatmap
I have a heatmap

I see a block of genes for a group of samples

How I can select a part of heatmap?

I mean is there any way to export a list of genes from rows of a heatmap in R?

I have tried which returns nothing

```
> library(ComplexHeatmap)

> ht = Heatmap(heat1)
> ht[1:35,1:5]

```

# Answer
You're using ComplexHeatmap, which is pretty easy to subset. Your subset operation creates a subset object of the heatmap, not of the data in the heatmap.

In case you don't want to run hclust outside of the CH, you can use `row_order` and `column_order` to re-order your data in the CH display order.

```
library(ComplexHeatmap)

ht = Heatmap(heat1)
r_o <- row_order(ht)
c_o <- column_order(ht)
heat1_reordered <- heat1[r_o, c_o]
heat1_reordered[1:35,1:5]

```

CH does say that clustering is oly "fixed" on `draw`ing, so you may want to run this before picking the row and column orders:

```
pdf(NULL)
ht <- draw(ht)
dev.off()

```
> 3 votes

---
Tags: r, heatmap
---


# Question
Title: How to convert CRAM file with 10x data in three fastq files
10x Genomics data are stored in three FASTQ files, besides the standard R1 and R2 reads, there is also a I1 file with some metadata. Sometimes however they are shipped in a single bam/cram file (e.g. the data from the Darwin Tree of Life). How can I convert the BAM or CRAM files to FASTQ?

-- edit --

Apparently, the I3 file is not 10x specific thing. It is an index generated when demultiplexing any Illumina data. See this question for a nice explanation.

# Answer
> 1 votes
All can be done with samtools. This is how the Darwin Tree of Life folks convert it:

```
samtools fastq -@4 -i \
  -1 ${sample}_S${tag}_L%03s_R1_001.fastq.gz \
  -2 ${sample}_S${tag}_L%03s_R2_001.fastq.gz \
  --i1 ${sample}_S${tag}_L%03s_I1_001.fastq.gz \
  --index-format i8 ${lane}.cram

```

# Answer
> 1 votes
You'll want to use the 10x provided bamtofastq tool to preserve the indices properly: https://github.com/10XGenomics/bamtofastq

---
Tags: sam, fastq, format-conversion, cram
---


# Question
Title: Range of values for VAF in heterozygous genotype
(VAF = variant allele frequency)

This question follows from this one.

What is the expected range of values for VAF of mutations with heterozygous genotype (i.e., 0/1 or 0|1)? Is it possible to have 0% and 100% (respectively having allele-specific depth of 0 in ALT or REF, with non-zero depth in the other allele)?

# Answer
From the link that I edited into your question (realized that the answer was there):

> VAF is the percentage of sequence reads observed matching a specific DNA variant divided by the overall coverage at that locus. Because NGS provides a near random sample, VAF is thus a surrogate measure of the proportion of DNA molecules in the original specimen carrying the variant. For constitutional genetic testing, VAF is a measure of diploid zygosity: heterozygous loci should be near 50% VAF, homozygous loci should be near 100%, and reference loci should be near 0%. Deviations from these three expected values should be considered suspicious as potential errors due to incorrect base calls or alignment. For somatic testing, contamination from normal cells and tumor heterogeneity combine to cause unpredictable VAFs. In some clinical scenarios, such as testing a patient for therapeutic resistance mutations, the desired sensitivity dictates that variants with low VAF be included in downstream analysis. This all means that in order to maintain an acceptable level of sensitivity, variant filtering for somatic variant must be highly permissive from a technical validity perspective.

This is assuming diploid humans. For other ploidies the expected proportions of heterozygosity will match the possible heterozygous genotypes (e.g. 33%/67% for triploid heterozygotes, 25%/75% for tetraploid, etc.).
> 1 votes

---
Tags: vcf, variant-calling, allele-frequency
---


# Question
Title: How Does Cox Proportional Hazard Model Describe Interaction between Covariates
I have read STDHA and the documentation for the coxph() function in R. I understand that the cox PH model explores how different covariates simultaneously impact survival. Based on the Lung data example in STDHA, a covariate can seem significantly influential in univariate analysis, but not in multivariate analysis. Following the lung data example, "age" is no longer significant in the multivariate cox model. Quite literally, this means that age does not impact survival when you account for sex and ph.ecog. However, how does the cox PH model evaluate and decide this? Under the hood, does the cox ph model stratify groups by sex, then evaluate if age is still impactful on survival (and then repeat for ph.ecog)? My intuition says no, because I don't think you can do that with continuous variables (so, in turn, how would the model stratify by age and then test if sex still impacts survival)? I don't need a deeply mathematical explanation, as I'm just a coder!

Ultimately, I want to understand if I can run coxph() on the female and male group separately, or if these results would be redundant since we have already added "sex" as a covariate when analyzing all patients.

This is a conceptual question, but please let me know if I can provide any code. If it means anything, I am using surv & survminer in R.

# Answer
The simple answer is that Cox PH is a flavor of multiple regression. Multiple regression does indeed try to subtract off the estimated effect of each covariate and look for associations of all other independent variables. Whether this works well in practice is anyone's guess. For a somewhat philosophical high-level discussion see here.

One notable issue is that "significant" is not necessarily a great way to think about statistics, according to professional statistical societies. So, like, if a p-value goes up or down a little bit that's kind of expected when you adjust for other variables, it's not really notable. For one thing, it's closely related to questions of power, driven largely by your sample size.

In simple terms, I would urge you to put all of your data in the model, find some set of independent variables and covariates that you think makes sense *before you start analyzing the data*, and then build your model accordingly. I don't think that anyone objects to e.g. plotting raw data or model residuals separately in males and females to demonstrate some point, but most associations that are worth caring about will be obvious in a combined dataset.

You *can* try to explicitly test for an interaction effect in the model. This can be interesting if there's some reason to expect it to be the case, or if there is very compelling evidence that there exists some interaction, but it has the negative of adding to model complexity and further reducing power.
> 2 votes

# Answer
"Under the hood" a Cox model doesn't stratify, consistent with your intuition. It iteratively solves a partial-likelihood vector score function as a function of regression-coefficient values.

It tries to minimize the sum, over all event times, of the difference between each covariate value for the case having the event and a risk-weighted average of the values for all at risk at that event time. The "risk-weighting" is a function of all current regression-coefficient estimates within an iteration. Iteration continues until regression-coefficient estimates are effectively constant. This nicely handles predictors of all types you might use in a multiple-regression model, including interaction terms, while also allowing for individuals that drop out of the study before experiencing the event.

There thus is nothing to be gained by doing separate models for different groups, as Maximilian Press also says in another answer. Stratified models are often used if the proportional hazards assumption doesn't hold, but those model all the data together and often are structured to provide the same regression-coefficient estimates for all strata while allowing for different baseline survival curves among strata.

I also recommend including as many predictors as possible that you might expect to be associated with outcome in your model, but with one caveat. You risk overfitting if you have fewer than 10-20 *events* per unpenalized predictor in your model (including levels beyond the first of categorical predictors and any interaction or nonlinear terms). So you might have to make some careful choices or use a penalized approach like ridge regression.

For further information, there are nearly 900 answered questions about Cox models on the Cross Validated site.
> 1 votes

---
Tags: r, statistics, statsmodels
---


# Question
Title: Extract only exon regions from GFF/GTF file with input bed regions
I have a bed file that contains a list of regions and I would like to extract only the exon regions from a GTF/GFF file that are contained in the given regions.

Is there a simple command or tool that I can use for that?

# Answer
You could use BEDOPS:

```
$ gtf2bed < annotations.gtf | grep -wF exon > exons.bed
$ sort-bed myRegions.unsorted.bed > myRegions.bed
$ bedops --element-of 100% exons.bed myRegions.bed > answer.bed

```

Or as a one-liner, which is even faster:

```
$ bedops --element-of 100% <(gtf2bed < annotations.gtf | grep -wF exon) <(sort-bed myRegions.unsorted.bed) > answer.bed

```

The `--element-of 100%` argument requires that the exon fall entirely within a region in the `myRegions` set. This is the most stringent requirement. The most relaxed requirement would be one base of overlap: `--element-of 1`. You can adjust to either extreme, or to values in between, depending on your criteria.
> 5 votes

# Answer
Bedtools `intersect` will do this:

> By far, the most common question asked of two sets of genomic features is whether or not any of the features in the two sets “overlap” with one another. This is known as feature intersection. bedtools intersect allows one to screen for overlaps between two sets of genomic features. Moreover, it allows one to have fine control as to how the intersections are reported. bedtools intersect works with both BED/GFF/VCF and BAM files as input.

https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html
> 2 votes

---
Tags: bed, gtf, exon, gff
---


title: subsetting more than one ccle sample

text:
# Question
Title: subsetting more than one ccle sample
I'm working with the CCLE dataset and I'm trying to subset the data to just 3 cell lines of interest but not sure how to it.

Currently, I have been subsetting one by one but is there an easier way?

```
ccle_panc <- ccle_dataset$data[ , ccle_dataset$colAnnotation$SAMPLE_NAME=="KP4", "", ""]

```

# Answer
Assuming you're trying to write R code, array references can use boolean values, and `in` is useful when you want to select from more than one value:

```
ccle_panc <- ccle_dataset$data[ , ccle_dataset$colAnnotation$SAMPLE_NAME %in% c("KP4", "DC1", "DC2")]

```

You can find more information about vector and array subsetting with boolean and numerical values here:

https://datacarpentry.org/R-ecology-lesson/01-intro-to-r.html#Subsetting\_vectors
> 0 votes

# Answer
Assuming that you have a dataframe/tibble situation; you can also use the tidyverse approach and the comment from @gringer:

```
ccle_dataset$data %>%
filter(SAMPLE_NAME %in% c("KP4", "DC1", "DC2"))

```
> 0 votes

---
Tags: r, rna-seq, subset
---


# Question
Title: Efficient way of mapping UniProt IDs to representative UniRef90 IDs?
I want to map ~500,000 UniProt IDs to their representative UniRef90 IDs.

What is the most efficient way of doing this?

Do we have a table with the mapping?

# Answer
You can do this directly on UniProt: https://www.uniprot.org/uploadlists/

Just paste or upload your list of UniProt IDs, and select "UniProtKB AC/ID" in the "From" field and "UniParc" in the "To" field

---

I've also written a script, pasted below, that can do this with some useful options:

```
$ uniprot_map.pl -h

uniprot_map.pl will connect to UniProt and map IDs from the input file to ehatever ID was requested.

USAGE:
       uniprot_map.pl [options] <ID_LIST>

OPTIONS:
      -h : Print this help and exit.
      -f : From ID type. This is the type of ID found in the input file.
      -F : Format. See http://www.uniprot.org/faq/28#batch_retrieval_of_entries for details.

      -l : When a requested ID maps to many output IDs, print only the first (probably the best) match.
      -L : List all recognised ID types.
      -u : Print unreviewed (trEMBL) prots as well. By default, if an ID maps to SWISSPROT and trEMBL IDs
           only the SWISSPROTs will be printed, use this flag to print trEMBL as well.
      -t : To ID. This is the ID type that the input IDs will be mapped to.
      -v : Verbose. Print progress messages to STDERR.

```

If you save the script as `uniprot_map.pl` in a directory in your PATH and make it executable, then you can use a file like this:

```
$ cat file
P31946
P62258

```

And feed it to the script with the appropriate options to get the mapping:

```
$ uniprot_map.pl -f ID -t UniParc file
P31946  UPI000013C714
P62258  UPI0000021A46
2 of 2

```

The script:

```
#!/usr/bin/perl
use strict;
use warnings;
use LWP::UserAgent;
use Getopt::Std;
my %opts;
getopts('Lsvuhf:F:t:l', \%opts);
my $base = 'http://www.uniprot.org';
my $verbose=$opts{v}||undef;
my $tool = 'mapping';
my (@NAMES,@names);
my $print_only_first=$opts{l}||undef;
my $allow_trembl=$opts{u}||undef;
list_ids() if $opts{L};
&usage() if $opts{h};

&usage() unless $ARGV[0];

####################################
# Collect the IDs we want to match #
####################################
foreach my $ar(@ARGV){
if (-e $ar){
open(A, $ar);
while(<A>){
    chomp;
    push @names,$_;
}
}
else{push @names,$ar};
}

######################################
# What are we matching FROM and TO ? #
######################################
my $from=$opts{f}||'ID';
my $to = $opts{t}||'ACC';

if($opts{t}){
$to=translate_terms($to);
}
if($opts{f}){
    $from=translate_terms($from);
}
#################
# Output format #
#################
my $format = $opts{F}||'tab';
#my @query=

my $agent = LWP::UserAgent->new;
push @{$agent->requests_redirectable}, 'POST';

my $params= {
'from' => $from,
'to' => $to,
'format' => $format,
    'query' => "@names"
};

print STDERR "Requesting...\t" if $verbose;
my $try_again=1;
my $response;
while($try_again==1){
    $response = get_response($params);
    if(defined($response->is_success)){$try_again=0}
elsif($verbose){print STDERR "Failed, retrying...\n"; }
}

if( $print_only_first)
{
if($response->is_success){
my %k;
my @a=split(/\n/,$response->content);
foreach my $line(@a){
    next if $line=~/From/;
    my @b=split(/\t/,$line);
    next if(defined($k{$b[0]}));
    $k{$b[0]}=$b[1];
}
map{print STDOUT "$_\t$k{$_}\n"}keys(%k);
}
else{
die 'Failed, got ' . $response->status_line . ' for ' . $response->request->uri . "\n";
}
}
else{
if($response->is_success){
my %map;
my @a=split(/\n/,$response->content);
#############################
    # remove "From   To"   line #
    #############################
shift(@a) if $a[0] =~/^From/i;
my %kk;
map{
    ########################################################
        # UniProt sometimes returns the SAME line twice    #
        ########################################################
    $kk{$_}++;
    #print STDOUT "$_\n" unless $kk{$_}>1;
    next if $kk{$_}>1;
    my @b=split(/\t/);
    push @{$map{$b[0]}},$b[1];
    }@a;

    foreach my $query (keys(%map)) {
    my %sw;
    foreach my $target (@{$map{$query}}) {
    #############################################
            # Does this prot have SWISSPROT IDs?        #
            #############################################
    if ($to eq 'ID')
    {
        my @b=split(/_/,$target);
        $sw{$target}++ if length($b[0])!=6;
    }
    }
    foreach my $target (@{$map{$query}}) {
    ################################################
    # If a protein has a least one SWISSPROT ID    #
    # skip printing its trEMBL ones unless -u      #
    ################################################
    if ($to eq 'ID' and not $allow_trembl and scalar(keys(%sw)) > 0) {
        print "$query\t$target\n" if defined($sw{$target});
    }
    else {
        print "$query\t$target\n";
        }
        }
    }

    ## Count how many were actually found
    my %found;
    foreach my $resp (@a){
    #print STDERR "R:$resp\t";
    foreach my $id (@names){
    #print STDERR "I:$id\n";
    $found{$id}++ if $resp=~/^$id\t/;
        }

    }
    print STDERR scalar(keys(%found)) . " of " . scalar @names . "\n";

    }
    else{
        die 'Failed, got ' . $response->status_line .
    ' for ' . $response->request->uri . "\n";
    }
}
################################################
# This will translate $to and $from values     #
# between the names used in uniprot flat files #
# and those used by the uniprot database       #
################################################
sub translate_terms{
    my $id=shift();
my %hash=(
"AC" => "ACC", "ID" => "ID", "UniParc" => "UPARC", "UniRef50" => "NF50", "UniRef90" => "NF90", "UniRef100" => "NF100", "EMBL/GenBank/DDBJ" => "EMBL_ID", "EMBL/GenBank/DDBJ CDS" => "EMBLD", "PIR" => "PIR", "UniGene" => "UNIGENE_ID", "GI number*" => "P_GI", "IPI" => "P_IPI", "RefSeq" => "P_REFSEQ_AC", "PDB" => "PDB_ID", "DisProt" => "DISPROT_ID", "HSSP" => "HSSP_ID", "DIP" => "DIP_ID", "MINT" => "MINT_ID", "MEROPS" => "MEROPS_ID", "PeroxiBase" => "PEROXIBASE_ID", "PptaseDB" => "PPTASEDB_ID", "REBASE" => "REBASE_ID", "TCDB" => "TCDB_ID", "Aarhus/Ghent-2DPAGE" => "AARHUS_GHENT_2DPAGE_ID", "ECO2DBASE" => "ECO2DBASE_ID", "World-2DPAGE" => "WORLD_2DPAGE_ID", "Ensembl" => "ENSEMBL_ID", "Ensembl Protein" => "ENSEMBL_PRO_ID", "Ensembl  Trscrpt" => "ENSEMBL_TRS_ID", "Ensembl  Genomes" => "ENSEMBLGENOME_ID", "Ensembl    Genomes" => "Protein", "Ensembl Genomes" => "Transcript", "GeneID" => "P_ENTREZGENEID", "GenomeReviews" => "GENOMEREVIEWS_ID", "KEGG" => "KEGG_ID", "TIGR" => "TIGR_ID", "UCSC" => "UCSC_ID", "VectorBase" => "VECTORBASE_ID", "AGD" => "AGD_ID", "ArachnoServer" => "ARACHNOSERVER_ID", "CGD" => "CGD", "ConoServer" => "CONOSERVER_ID", "CYGD" => "CYGD_ID", "dictyBase" => "DICTYBASE_ID", "EchoBASE" => "ECHOBASE_ID", "EcoGene" => "ECOGENE_ID", "euHCVdb" => "EUHCVDB_ID", "EuPathDB" => "EUPATHDB_ID", "FlyBase" => "FLYBASE_ID", "GeneCards" => "GENECARDS_ID", "GeneDB_Spombe" => "GENEDB_SPOMBE_ID", "GeneFarm" => "GENEFARM_ID", "GenoList" => "GENOLIST_ID", "H-InvDB" => "H_INVDB_ID", "HGNC" => "HGNC_ID", "HPA" => "HPA_ID", "LegioList" => "LEGIOLIST_ID", "Leproma" => "LEPROMA_ID", "MaizeGDB" => "MAIZEGDB_ID", "MIM" => "MIM_ID", "MGI" => "MGI_ID", "NMPDR" => "NMPDR", "Orphanet" => "ORPHANET_ID", "PharmGKB" => "PHARMGKB_ID", "PseudoCAP" => "PSEUDOCAP_ID", "RGD" => "RGD_ID", "SGD" => "SGD_ID", "TAIR" => "TAIR_ID", "TubercuList" => "TUBERCULIST_ID", "WormBase" => "WORMBASE_ID", "WormBase Trscrpt" => "WORMBASE_TRS_ID", "WormBase Protein" => "WORMBASE_PRO_ID", "Xenbase" => "XENBASE_ID", "ZFIN" => "ZFIN_ID", "eggNOG" => "EGGNOG_ID", "HOGENOM" => "HOGENOM_ID", "HOVERGEN" => "HOVERGEN_ID", "OMA" => "OMA_ID", "OrthoDB" => "ORTHODB_ID", "ProtClustDB" => "PROTCLUSTDB_ID", "BioCyc" => "BIOCYC_ID", "Reactome" => "REACTOME_ID", "CleanEx" => "CLEANEX_ID", "GermOnline" => "GERMONLINE_ID", "DrugBank" => "DRUGBANK_ID", "NextBio" => "NEXTBIO_ID"
);
defined($hash{$id}) ?
return($hash{$id}) :
return($id) ;
}

sub get_response{
    my $params=shift;
my $response = $agent->post("$base/$tool/", $params);
    while (my $wait = $response->header('Retry-After')) {
    print STDERR "Waiting ($wait)...\n";
sleep $wait;
    $response = $agent->get($response->base);
}
return($response);

}

sub usage{
    $0=~s/.+?([^\/]+)$/$1/;

    print STDERR <<EndOfHelp;

$0 will connect to UniProt and map IDs from the input file to ehatever ID was requested.

USAGE:
       $0 [options] <ID_LIST>

OPTIONS:
      -h : Print this help and exit.
      -f : From ID type. This is the type of ID found in the input file.
      -F : Format. See http://www.uniprot.org/faq/28#batch_retrieval_of_entries for details.

      -l : When a requested ID maps to many output IDs, print only the first (probably the best) match.
      -L : List all recognised ID types.
      -u : Print unreviewed (trEMBL) prots as well. By default, if an ID maps to SWISSPROT and trEMBL IDs
           only the SWISSPROTs will be printed, use this flag to print trEMBL as well.
      -t : To ID. This is the ID type that the input IDs will be mapped to.
      -v : Verbose. Print progress messages to STDERR.
EndOfHelp

    exit(0);

}
sub list_ids{
    open(HELP, "| less") ;
    print HELP <<EndOfHelp;

Database            Abbreviation
===================== UniProt =============================
UniProtKB AC/ID :   ACC+ID
UniProtKB AC    :   ACC
UniProtKB ID    :   ID
UniParc         :   UPARC
UniRef50    :   NF50
UniRef90    :   NF90
UniRef100   :   NF100

=============== Other sequence databases ==================
EMBL/GenBank/DDBJ:  EMBL_ID
EMBL/GenBank/DDBJ CDS:  EMBLD
PIR     :   PIR
UniGene     :   UNIGENE_ID
GI number*  :   P_GI
IPI     :   P_IPI
RefSeq      :   P_REFSEQ_AC
HGNC name       :       GENENAME
================== 3D structure databases =================
PDB     :   PDB_ID
DisProt     :   DISPROT_ID
HSSP        :   HSSP_ID

============ Protein-protein interaction databases ========
DIP     :   DIP_ID
MINT        :   MINT_ID

============ Protein family/group databases ===============
MEROPS      :   MEROPS_ID
PeroxiBase  :   PEROXIBASE_ID
PptaseDB    :   PPTASEDB_ID
REBASE      :   REBASE_ID
TCDB        :   TCDB_ID

================== 2D gel databases =======================
Aarhus/Ghent-2DPAGE:  AARHUS_GHENT_2DPAGE_ID
ECO2DBASE   :   ECO2DBASE_ID
World-2DPAGE    :   WORLD_2DPAGE_ID

============= Genome annotation databases =================
Ensembl         :   ENSEMBL_ID
Ensembl Protein :   ENSEMBL_PRO_ID
Ensembl Trscrpt :   ENSEMBL_TRS_ID
Ensembl Genomes :   ENSEMBLGENOME_ID
Ensembl Genomes :   Protein
Ensembl Genomes :   Transcript
GeneID          :   P_ENTREZGENEID
GenomeReviews   :   GENOMEREVIEWS_ID
KEGG        :   KEGG_ID
TIGR        :   TIGR_ID
UCSC        :   UCSC_ID
VectorBase  :   VECTORBASE_ID

============ Organism-specific gene databases =============
AGD     :   AGD_ID
ArachnoServer   :   ARACHNOSERVER_ID
CGD     :   CGD
ConoServer  :   CONOSERVER_ID
CYGD        :   CYGD_ID
dictyBase   :   DICTYBASE_ID
EchoBASE    :   ECHOBASE_ID
EcoGene     :   ECOGENE_ID
euHCVdb     :   EUHCVDB_ID
EuPathDB    :   EUPATHDB_ID
FlyBase     :   FLYBASE_ID
GeneCards   :   GENECARDS_ID
GeneDB_Spombe   :   GENEDB_SPOMBE_ID
GeneFarm    :   GENEFARM_ID
GenoList    :   GENOLIST_ID
H-InvDB     :   H_INVDB_ID
HGNC        :   HGNC_ID
HPA     :   HPA_ID
LegioList   :   LEGIOLIST_ID
Leproma     :   LEPROMA_ID
MaizeGDB    :   MAIZEGDB_ID
MIM     :   MIM_ID
MGI     :   MGI_ID
NMPDR       :   NMPDR
Orphanet    :   ORPHANET_ID
PharmGKB    :   PHARMGKB_ID
PseudoCAP   :   PSEUDOCAP_ID
RGD     :   RGD_ID
SGD     :   SGD_ID
TAIR        :   TAIR_ID
TubercuList :   TUBERCULIST_ID
WormBase    :   WORMBASE_ID
WormBase Trscrpt:   WORMBASE_TRS_ID
WormBase Protein:   WORMBASE_PRO_ID
Xenbase     :   XENBASE_ID
ZFIN        :   ZFIN_ID

================= Phylogenomic databases ==================
eggNOG      :   EGGNOG_ID
HOGENOM     :   HOGENOM_ID
HOVERGEN    :   HOVERGEN_ID
OMA     :   OMA_ID
OrthoDB     :   ORTHODB_ID
ProtClustDB :   PROTCLUSTDB_ID

=================== Enzyme and pathway ====================
BioCyc      :   BIOCYC_ID
Reactome    :   REACTOME_ID

================= Gene expression databases ===============
CleanEx     :   CLEANEX_ID
GermOnline  :   GERMONLINE_ID

========================== Other ==========================
DrugBank    :   DRUGBANK_ID
NextBio     :   NEXTBIO_ID

EndOfHelp
close(HELP);
    exit(0);
}

```

If you have too many IDs to do in a single batch, you could either simply split the file into smaller ones and run the script on each file in a loop or modify the script to do it.

Split the input:

```
$ wc -l names
$ wc -l names
500000 names
## split into 500 files of 1K lines each
$ split -l 1000 names chunk
## iterate over each chunk of names
$ for i in chunk*; do uniprot_map.pl -f ID -t UniParc "$i" >> mapped; done

```
> 2 votes

---
Tags: uniprot
---


# Question
Title: How to simulate NGS reads, controlling sequence coverage?
I have a FASTA file with 100+ sequences like this:

```
>Sequence1
GTGCCTATTGCTACTAAAA ...
>Sequence2
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTA
......

```

I also have a text file like this:

```
Sequence1 40
Sequence2 30
......

```

I would like to simulate next-generation paired-end reads for all the sequences in my FASTA file. For `Sequence1`, I would like to simulate at 40x coverage. For `Sequence2`, I would like to simulate at 30x coverage. In other words, I want to control my sequence coverage for each sequence in my simulation.

**Q:** What is the simplest way to do that? Any software I should use? Bioconductor?

# Answer
I am not aware of any software that can do this directly, but I would split the fasta file into one sequence per file, loop over them in BASH and invoke ART the sequence simulator (or another) on each sequence.

For more information about ART, please see their paper here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3278762/
> 2 votes

# Answer
I am working on a Illumina sequencing simulator for metagenomics: InSilicoSeq

It is still in alpha release and very experimental, but given a multi-fasta and an abundance file, it will generate reads from your input genomes with different coverages.

From the documentation:

```
iss generate --genomes genomes.fasta --abundance abundance_file.txt \
    --model_file HiSeq2500 --output HiSeq_reads

```

Where:

```
# multi-fasta file
>genome_A
ATGC...
>genome_B
CCGT...
...

# abundance file (total abundance must be 1!)
genome_A    0.2
genome_B    0.4
...

```

I didn't design it to work with coverage but rather abundance of the genome in a metagenome, so you might have to do a tiny bit of math ;)
> 7 votes

# Answer
The polyester bioconductor package can do this. It says it simulates RNA-seq reads, but I don't know if that's really any different from other NGS reads.

It can use a range of error and bias models, or learn them from a dataset.
> 4 votes

# Answer
This python script takes a fasta file and tsv file with counts and prints the sequences in the fasta files that many times as it is specified in the tsv file (assuming the format in the question). So if `bar.tsv` and `foo.fasta` will be your files:

```
from Bio import SeqIO

repeat = {}
for line in open("bar.tsv"):
    seq_id, coverage = line.split()
    repeat[seq_id] = int(coverage)

for seq_record in SeqIO.parse(foo.fasta, "fasta"):
    for i in range(repeat.get(seq_record.name, 0)):
        print(">",seq_record.name,"_",i,sep='')
        print(seq_record.seq)

```
> 4 votes

# Answer
The wgsim package by Heng Li (of BWA and samtools fame) is my go-to tool for simulating Illumina reads. It doesn't provide any convenient way to simulate differential coverage across different sequences, but it shouldn't be to hard to run wgsim multiple times, generating the desired level of coverage for each sequence of interest.

I would implement a Python script to slurp up your test file, and call wgsim (using the `subprocess` module) for each sequence. This will probably require you to have each sequence in a separate file. :-(
> 4 votes

# Answer
If you're ok with some randomness you can generate reads from your sequence file using a Poisson random variable. You'll need to do some math to figure out what value of lambda to use in order for the expected coverage at each base pair in your read to match what you set in your text file.

For example you have a sequence S of length 1,000, a read length of 50, and an insert size of 100. For each base b in S generate a Poisson random variable p. You will then generate p reads from base b to b+50. Then, generate the paired read starting at b+50+100.

Again, you would have to play with it to figure out what lambda to use but this would give you basically what you want, as long as you're ok with not having exactly the coverage you're targeting for each read.
> 2 votes

# Answer
Simulating NGS reads while controlling sequence coverage is now easy with RNFtools (from version 0.3.1). See the tutorial, especially section Sequence extraction.

**Environment preparation**

First, install BioConda and add the required channels. Then either install RNFtools in the default Conda environment

```
conda install rnftools

```

or create and activate a separate Conda environment (preferable)

```
conda create -n rnftools rnftools
source activate rnftools

```

**Simulation**

Assume that you have a reference file `ref.fa` and a tab-separated coverage file `coverage.tsv` (e.g., those from your example). Then the following RNFtools `Snakefile` will do the job you want:

```
import rnftools
import csv

rnftools.mishmash.sample("simulation_with_coverage_control", reads_in_tuple=1)

fa = "ref.fa"
tsv = "coverage.tsv"

with open(tsv) as f:
    table = csv.reader(f, delimiter='\t')
    for seqname, cov in table:

        rnftools.mishmash.DwgSim(
            fasta=fa,
            sequences=[seqname],
            coverage=float(cov),
            read_length_1=10, # quick test with supershort reads
            read_length_2=0,
        )

include: rnftools.include()
rule: input: rnftools.input()

```

When you save this file (`Snakefile`) and run `snakemake`, RNFtools will simulate reads using DWGsim with the coverages defined your text file, and save all the simulated reads in `simulation_with_coverage_control.fq`.

You can play with all the parameters. In particular, you can use a different simulator (e.g., Art-Illumina using `rnftools.mishmash.ArtIllumina`). See the RNFtools documentation for more information.
> 2 votes

# Answer
Another next-gen read simulation tool is gemsim. I haven't tested it, but I would be interested if anyone has had any experience with it.
> 0 votes

# Answer
Go to ART website. Download your package, binary or source:  Uncompress the .gz file, open the folder "art\_bin\_MountRainier", then "examples". Open the text file "run\_test\_examples\_illumina.sh". This file contains numerous example of samples codes. That will give you an idea on how you can use the software. Here is an example:

```
# 4) simulation of paired-end reads of 100bp with the mean fragment size 500 and standard deviation 10
#   using the built-in seperated quality profiles for A, C, G, and T
$art -ss HS20 -i ./testSeq.fa -o ./paired_end_sep -l 100 -f 10 -p -m 500 -s 10 -sp -sam
#convert both aln files to a bed file
../aln2bed.pl paired_end_sep.bed paired_end_sep1.aln paired_end_sep2.aln

```
> 0 votes

# Answer
You can split your FASTA file sequence-by-sequence using

```
split -a 6 -p '^>' your_file.fa seq_

```

and then use any existing read simulator supporting coverage (ART, DWGsim, etc.). If you want to have all the reads mixed (not ordered by the original sequence), you can use RNFtools.

**Edit 1:**

As @terdon pointed out, the previous command works on OS X only. An analogical one liner for Linux (but with a slightly different naming scheme using numbers rather than letters) can be

```
csplit -f seq_ -n 6 your_file.fa '/^>/' {*}

```

To make this command work also on OS X, one needs to install coreutils (e.g., using brew) and then use `gcsplit` instead of `csplit`.

**Edit 2:**

Once FASTA is split by sequences, the simulation becomes straightforward and many different approaches can be used. My favorite one is using GNU Parallel. Imagine that you have your coverages in a text file called `covs.txt` on separate lines and in the same order as the sequences in `your_file.fa`, e.g.,

```
40
30
...

```

Then you can simulate reads from the original sequences using DWGsim by

```
ls -1 seq_* | paste covs.txt - \
    | parallel -v --colsep '\t' dwgsim -1 100 -2 100 -C {1} {2} sim_{2}

```

and merge the obtained FASTQ files using:

```
cat sim_seq_*.bwa.read1.fastq > reads.1.fq
cat sim_seq_*.bwa.read2.fastq > reads.2.fq

```

One possible danger of this approach is that we have assumed that the number of seq\_\* files is the same as the number of lines in `covs.txt`, which might not be true (by mistake). We should check this prior to the simulation step, e.g., by:

```
[[ "$(ls -1 seq_* | wc -l)" == "$(cat covs.txt | wc -l)" ]] \
    || echo "Incorrent number of lines in covs.txt"

```

Another caveat is that the simulated reads are not in a random order (they are grouped by source sequences).
> -1 votes

---
Tags: fasta, ngs, simulated-data
---


# Question
Title: How to subset a GRanges object based on a specific genomic window of interest?
**Goal:** To keep only rows of a GRanges object if their start coordinate is within 1 Mb of the start coordinate of a gene of interest.

I am aware of some functions that come close:

* `GenomicRanges::nearestKNeighbors()` (and related `distanceToNearest()`, `nearest()`, `precede()`, `follow()`)
* `GenomicRanges::subsetByOverlaps()`
* `GenomicRanges::window()`

And I am aware of the UNIX command-line 'Bedtools' tools: `window`, `closest`, and `intersect`, however I wish to accomplish this task using R only and the above functions have not been doing exactly as I would like.

Here is an example GRanges dataset:

```
> gr <- GRanges(
        seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
        ranges = IRanges(start = seq(300000, 2000000, length=10), end = seq(700000, 6000000, length=10), names = head(letters, 10)),
        strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
        score = 1:10, GC = seq(1, 0, length=10), GeneSymbol = c(paste("gene", seq(1:10), sep="")))

> gr
GRanges object with 3 ranges and 3 metadata columns:
    seqnames          ranges strand |     score        GC  GeneSymbol
       <Rle>       <IRanges>  <Rle> | <integer> <numeric> <character>
  a     chr1   300000-700000      - |         1  1.000000       gene1
  b     chr2  488888-1288888      + |         2  0.888889       gene2
  c     chr2  677777-1877777      + |         3  0.777778       gene3
  d     chr2  866666-2466666      * |         4  0.666667       gene4
  e     chr1 1055555-3055555      * |         5  0.555556       gene5
  f     chr1 1244444-3644444      + |         6  0.444444       gene6
  g     chr3 1433333-4233333      + |         7  0.333333       gene7
  h     chr3 1622222-4822222      + |         8  0.222222       gene8
  i     chr3 1811111-5411111      - |         9  0.111111       gene9
  j     chr3 2000000-6000000      - |        10  0.000000      gene10
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

```

Ideally, by specifying the name (e.g., `"gene3"`) or start coordinate (e.g., `chr2` *AND* `677777`) of the gene of interest AND a window size (e.g., `1000000` or `1 Mb` upstream *AND* downstream) to a function, I hope to get an output like below:

```
GRanges object with 3 ranges and 3 metadata columns:
    seqnames          ranges strand |     score        GC  GeneSymbol
       <Rle>       <IRanges>  <Rle> | <integer> <numeric> <character>
  b     chr2  488888-1288888      + |         2  0.888889       gene2
  c     chr2  677777-1877777      + |         3  0.777778       gene3
  d     chr2  866666-2466666      * |         4  0.666667       gene4
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

```

Notice that only three rows from the original GRanges object remain in this example output because the start coordinates for all other genes were either (i) on different chromosomes or (ii) were *more than 1 Mb away* from the start coordinate of gene3.

---

I have tried various methods to isolate the rows within 1 Mb of my gene of interest with no luck:

```
> gr1 <- gr[grep("gene3", gr$GeneSymbol),]
> grOut <- subsetByOverlaps(gr, gr1, maxgap = 1000000, invert=T)

> grOut <- gr[mcols(gr)$GeneSymbol == "gene3" & start(gr) < 677777-1000000 & end(gr) > 1877777+1000000]

> grOut <- window(gr, GRanges("chr2", IRanges(start=677777-1000000, end=1877777+1000000)))
Error: 'start' must be a single number or NA

> grOut <- gr[IRanges(start=677777, end=1877777)]
Error: subscript contains out-of-bounds ranges

> grOut <- nearest(gr, GRanges("chr2", IRanges(start=677777-1000000, end=1877777+1000000)), select="all", ignore.strand=T)

```

The last attempt recorded above was likely the closest I've gotten to my desired output by storing the index of rows in the GRanges object falling within the range specified, but it presents this index in the form of a SortedByQueryHits object, not a GRanges object as I require. It seems simple enough a task, but I haven't solved it yet; any thoughts or solutions?

# Answer
Specifying the range of interest (+/- 1 Mb of start coordinate) as a new GRanges object and then subsetting the original GRanges object based on the overlap of genes with the query works:

```
library("GenomicRanges")

> chrPos <- as.character(seqnames(gr))[gr$GeneSymbol=="gene3"]
> startPos <- start(ranges(gr))[gr$GeneSymbol=="gene3"]
> range <- 1000000

> q <- GRanges(seqnames=chrPos,
               ranges=IRanges(start=startPos-range, end=startPos+range))
> grOut <- subsetByOverlaps(gr, q)
> grOut
GRanges object with 3 ranges and 3 metadata columns:
    seqnames         ranges strand |     score        GC  GeneSymbol
       <Rle>      <IRanges>  <Rle> | <integer> <numeric> <character>
  b     chr2 488888-1288888      + |         2  0.888889       gene2
  c     chr2 677777-1877777      + |         3  0.777778       gene3
  d     chr2 866666-2466666      * |         4  0.666667       gene4
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

```

Credit to @zx8754 for inspiring this solution.
> 4 votes

# Answer
If preferred, this can be achieved with base R functions:

```
# First store the chromosome and start position of your gene of interest in variables
> chrPos <- as.character(seqnames(gr))[gr$GeneSymbol=="gene3"]
> startPos <- start(ranges(gr))[gr$GeneSymbol=="gene3"]

# Then add a column to the GRanges object with distances from the start position of each
# gene in the GRanges object to the start position of your gene of interest
> gr$distance <- start(ranges(gr))-startPos

# Finally subset the GRanges object based on your desired window of interest (must be
# symmetrical —— e.g., +/- 1 Mb —— for the following to work)
> grOut <- gr[seqnames(gr) %in% chrPos & (abs(gr$distance) <= 1000000)]
> grOut
GRanges object with 3 ranges and 4 metadata columns:
    seqnames         ranges strand |     score        GC  GeneSymbol  distance
       <Rle>      <IRanges>  <Rle> | <integer> <numeric> <character> <integer>
  b     chr2 488888-1288888      + |         2  0.888889       gene2   -188889
  c     chr2 677777-1877777      + |         3  0.777778       gene3         0
  d     chr2 866666-2466666      * |         4  0.666667       gene4    188889
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

```

This strategy provides the added benefit of adding a 'distance-from-gene' (`gr$distance`) column to the dataset which can be useful for visualizations.
> 0 votes

---
Tags: r, filtering, genomicranges
---


# Question
Title: How do I build a glycosylated protein from PDB files manually?
I have a large protein PDB file (~3500 residues) and a couple of PDB files containing glycan structures I sourced from GLYCAM.

How can I attach the N and O-glycans onto each protein residue with software? In PyMol, I am stuck at loading up the protein PDB file and glycan file together, but I don't know how to put the two elements together. Are there any commands I can use? or tutorial resources I can take a look at? I'm open to trying other platforms.

I've tried glycan modelling webtools (GLYCAM, CHARMM-GUI) but they can't seem to process my large glycoprotein - leading to errors or extended waiting times which also end in failure.

# Answer
You can use ChimeraX.

1. Open your glycans and your protein pdb files as different models, then move the models of your glycans until they are near the residues you want (in the Toolbar: Right Mouse -\> Move Model).
2. Then use `combine` to combine one of your glycan model with the protein
3. Finally Use Tools -\> Structure Editing -\> Build Structure and its Adjust Bonds menu to add covalent bonds. You can also add or delete atoms this way.

This is a straightforward but complicated and tedious way of doing this, and it will give you unideal structures. There should be an equivalent way to proceed using PyMol. For a more automatized way I can only recommend having a look at PyRosetta. There is a Jupyter Notebook tutorial about adding glycosylation.

There may be other tools out there for doing this, like MoleculeKit, but I'm less aware of them. A nasty solution would be to cut your protein in smaller pieces, glycosylate each part with CHARMM-GUI, then reassemble them.
> 2 votes

# Answer
*answer from @coot, converted from comments (please edit this answer to improve it if you can):*

Manual model building is a common thing in structural biology. Usually, Coot is used.
> 1 votes

---
Tags: protein-structure, pymol, modelling
---


# Question
Title: Resources to learn genome assembly workflow for small genomes (like viruses)
I have sequencing data of a few samples of a DNA virus.

I'd like to learn *de novo assembly* of 'short read' data, to construct a scaffold and then count the abundance of each strain in the data.

I heard about `SPAdes` as a good choice for these kinds of very short genomes. and also `BBmap` for statistics related to contigs.

I new to *de novo* assembly could you also recommend introductory material? For example a full workflow post/paper on genome assembly.

I have some experience with RNA-seq and small RNA-seq data.

**Note**: I am seeking an overview and test data prior reading Spades and BBmap manuals. Thus I need to learn RNA-seq analysis rather than a `aligner`.

cross-posted in biostars

# Answer
I would recommend a specific tool for viral meta genomics rather than using applying a general tool into a specific technical challenge. Here are the ones discussed on Bioinfo SE for *de novo* metagenomic assembly tested/being applied to viruses. It is not an exclusive list and needs further discussion.

SAVAGE https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411778/ Baaijens et al. De novo assembly of viral quasispecies using overlap graphs Genome Research 2017 May; 27: 835–848.

S-aligner https://www.biorxiv.org/content/biorxiv/early/2021/02/02/2021.02.02.429443.full.pdf Juanjo Bermúdez s-aligner: a greedy algorithm for non-greedy de novo genome assembly 2021 bioRχiv February: 1-7

Rey http://denovoassembler.sourceforge.net Sébastien Boisvert et al Genome Biology. 2012 13:R122, December
> 1 votes

---
Tags: genome, assembly, covid-19, scaffold, virus
---


# Question
Title: How to represent a summary of disagreement between two secondary structure predictors given predictions for many proteins?
I have 6000 predicted protein structures.

I am comparing secondary structures predicted by an algorithm with the correct secondary structure assignments.

I need to show which positions in a predicted protein chain have the most differences as to the true protein chain.

How can I do that?

How to represent a summary of disagreement between two secondary structure predictors given predictions for many proteins?

I already prepared confusion matrices. I was looking for something graphical like chart or histogram, etc..

# Answer
(This comment got too long)

> positions in a predicted protein chain

You do mean *position* in SS sequence, not *coordinates* in tertiary, yes? From the sounds of it you are talking about SS —not tertiary— and want the analysis to be independent of tertiary.

Therefore, you are comparing for each protein two 1D arrays of sequential elements containing between 3 and 6 states, one for each residue. Technically, you could add another state, a chain break, and have two colossally long 1D arrays. This very much becomes a **statistics** classification topic.

The properties of these arrays are:

* The values are sequential.
* Neighbouring elements influence each other as a span of residues makes a SS element. Thus a correct classification will propagate. This will mean that all metrics will be hugely inflated.
* Assumption: the divergence is likely on the ends of a span. This could be tested and plotted!

The simplest thing would be collapse the states into a binary loop vs non loop and classify each element of the two arrays and do a **typical confusion matrix** and its metrics (see Wikipedia). This would have the issue that neighbours influence each other.

For a figure, given N states, one could simply plot a *N* × *N* heatmap of the true states vs. predicted states. This would show a hot diagonal of the correct allocations, but will show which states are more prone to misallocation. The frequency of the states may need normalising... but this would still be affected by the length of the SS spans. Therefore, a more sophisticated *model* may be needed for better normalisation.

For more advanced classifications, more advanced algorithms would be needed, but it still remains a stats problem. It might be worth asking the question again under that discussion framework or look up a relevant multi-state classifier.
> 3 votes

---
Tags: protein-structure
---


# Question
Title: Determine reference for reference-compressed SRA file
*Note: this question was also asked on Github.*

I have 241 SRA files that appear to be reference compressed. I didn't even know this was a thing until I tried to convert them to Fastq files without an internet connection. I got the "name not found while resolving..." error described at the bottom of this page when I attempted to do the conversion.

The instructions say:

> *If you have elected to prevent the toolkit from contacting NCBI, you will need to manually acquire the reference(s) here*

We have many NCBI databases mirrored locally, including RefSeq, so it's likely that we already have all the references. I have two questions:

* How do I determine the reference used to compress a reference-compressed SRA file? (Note this needs to be scripted for 241 SRA files, so manual inspection is out of the question.)
* Once I've determined the correct reference and ensured it is available locally, how to I tell `fastq-dump` where to find it?

# Answer
Run align-info to see the references used. Run prefetch \</path/to/sra/file\> to download missed references.
> 1 votes

---
Tags: file-formats, fastq, format-conversion, sratoolkit
---