# Question
Title: What column and row naming requirements exist with Seurat (context: when loading SPLiT-Seq data)
I'm trying to use Seurat for the first time and am learning about single-cell analysis for the first time, and I'm doing so with split-seq data. (Full disclosure: I'm also a lightweight when it comes to R.)

When Read10X didn't work on the mtx and associated files output by the split-seq-pipeline no matter how much I tried to make it look exactly like a 10x dataset (based on numerous examples and descriptions), I tried converting my sparse market matrix format into a dense tsv format (see update below for a sample of the file) and I read it in using:

```
dge_data <- as.matrix(read.table("DGE.tsv", sep="\t", header=TRUE))
pbmc <- CreateSeuratObject(counts = dge_data, min.cells = 3, min.features  = 200, project = "SPLITSEQ_TT", assay = "RNA")

```

I started following chapter 9 of the scRNA-Seq course I found here. I ran into some issues which I was able to resolve by tweaking my mtx to tsv conversion script, such as including the chromosome in the row names (so that the mitochondrial QC step would work), but I ran into some errors that I assume have to do with the format of the gene names/IDs and possibly some issues stemming from cell names/IDs (which I've thus far been able to work through I think).

However, given errors like the following, and the fact that I can't seem to find specifications on these sorts of row/column name requirements (e.g. no commas and must include at least the "MT" chromosome name for QC steps using mitochondrial genes), I was wondering if I could be directed to somewhere where those specifications are laid bare. For the most part, the gene names have not been a problem, but I run afoul in this step (where it appears to complain about a comma):

```
> FeatureScatter(object = pbmc, feature1 = "hg38-4,ENSG00000074966,TXK", feature2 = "PC_1")
Error in parse(text = x) : <text>:1:7: unexpected ','
1: hg38.4,
          ^

```

Also, when I look at how CellsByIdentities behaves (while I was working through other errors/warnings), I realized I'm not sure that the column names are being processed correctly due to the underscores (followed by a number (e.g. `_0`)), because my understanding is that the `_#` portion of the split-seq cell barcode indicates the first well barcode, which would be needed to uniquely identify the cells in the table. I had retrieved the cell barcode IDs from the cell\_metadata.csv file output by split-seq. Should I modify these column names? Because I infer from this output from CellsByIdentities that it's only using the sequence portion of the cell barcodes to uniquely identify cells, which I don't think is correct:

```
> CellsByIdentities(object=pbmc, cells=colnames(pbmc))
$AAACATCGAAACATCG
[1] "AAACATCGAAACATCG_0" "AAACATCGAAACATCG_1"

$AAACATCGAACGTGAT
[1] "AAACATCGAACGTGAT_0"

$AACGTGATAAACATCG
[1] "AACGTGATAAACATCG_0" "AACGTGATAAACATCG_1"

$AACGTGATAACGTGAT

[1] "AACGTGATAACGTGAT_0" "AACGTGATAACGTGAT_1"

```

The other assumptions have to do with data size. I was running a small set of test data, and I had to adjust some parameters. So I would like to know if there's a hard lower limit on the data. And do I need to restructure my data to separate counts of genes from different species (e.g. mouse & human)? I.e. Should I create separate matrices that have only the human-identified cells and gene annotations to use it with Seurat?

I can't seem to find these specifications, presumably because Seurat was designed exclusively around 10x data and that the tool wasn't tested on data from other sources. Perhaps I'm worrying too much, but without knowing what may rely on column/row name format assumptions, I can't be sure my data is being analyzed properly.

**UPDATE**

Here is the head of the tsv file I originally tried to load when I ran into various errors that I suspect were due to the formatting of the row and column names (note, there are only 7 cells total, and it's a small but real test dataset from our core facility):

```
AAACATCGAAACATCG_0  AAACATCGAAACATCG_1  AAACATCGAACGTGAT_0  AACGTGATAAACATCG_0  AACGTGATAAACATCG_1  AACGTGATAACGTGAT_0  AACGTGATAACGTGAT_1
hg38_X,ENSG00000000003,TSPAN6   8   2   2   26  15  17  13
hg38_X,ENSG00000000005,TNMD 0   0   0   0   0   2   0
hg38_20,ENSG00000000419,DPM1    22  26  27  73  43  50  28
hg38_1,ENSG00000000457,SCYL3    6   6   5   12  7   19  8
hg38_1,ENSG00000000460,C1orf112 23  10  24  59  34  49  21
hg38_1,ENSG00000000971,CFH  101 44  129 331 151 270 122
hg38_6,ENSG00000001036,FUCA2    33  28  24  82  58  74  36
hg38_6,ENSG00000001084,GCLC 30  16  14  72  27  54  42
hg38_6,ENSG00000001167,NFYA 14  9   17  41  28  31  12

```

Yesterday, I did end up trying to change the column and row names to the following:

```
AAACATCGAAACATCG-0  AAACATCGAAACATCG-1  AAACATCGAACGTGAT-0  AACGTGATAAACATCG-0  AACGTGATAAACATCG-1  AACGTGATAACGTGAT-0  AACGTGATAACGTGAT-1
hg38-X.ENSG00000000003.TSPAN6   8   2   2   26  15  17  13
hg38-X.ENSG00000000005.TNMD 0   0   0   0   0   2   0
hg38-20.ENSG00000000419.DPM1    22  26  27  73  43  50  28
hg38-1.ENSG00000000457.SCYL3    6   6   5   12  7   19  8
hg38-1.ENSG00000000460.C1orf112 23  10  24  59  34  49  21
hg38-1.ENSG00000000971.CFH  101 44  129 331 151 270 122
hg38-6.ENSG00000001036.FUCA2    33  28  24  82  58  74  36
hg38-6.ENSG00000001084.GCLC 30  16  14  72  27  54  42
hg38-6.ENSG00000001167.NFYA 14  9   17  41  28  31  12

```

And it seems to have fixed most of the issues. The violin plots now have curvy brown swaths down the center that were not there before. And in particular, when I run CellsByIdentities manually on the column names now, I get what looks correct (as opposed to what I pasted yesterday above (which appeared to group the 7 cells in 3 pairs and one individual):

```
CellsByIdentities(object = pbmc, cells = colnames(pbmc))
$SPLITSEQ_TT
[1] "AAACATCGAAACATCG.0" "AAACATCGAAACATCG.1" "AAACATCGAACGTGAT.0" "AACGTGATAAACATCG.0" "AACGTGATAAACATCG.1" "AACGTGATAACGTGAT.0"
[7] "AACGTGATAACGTGAT.1"

```

I think that the remaining issues are likely due to the small data size. For reference as to the data size, this is the first data line from the original mtx file:

```
7 38526 194539

```

Please note also that I'm not looking for how to get this particular dataset to work. I am (as a member of the bioinformatics team that supports users of our core sequencing facility) trying to establish the requirements for the conversion script so that when we deliver SPLiT-Seq results, our output can be directly used with Seurat by the users without issues or need for further conversion. So I specifically need to know what makes row or column names valid or invalid, as opposed to what to change WRT this specific dataset. We may be supplied with annotations from various sources that may have different characters in the gene symbols, or accession numbers, or even chromosome names.

# Answer
`Read10X()` works on `mtx` files, that is not the root of your problem. The problem is that you are pushing non-10x-data to look and work like 10x-data.

You can just use `CreateSeuratObject()` and feed your count matrix to this function. I am not sure if this particular function accepts `mtx` files, could not find this in the documentation, but in the worst case you will just need to supply a plain count matrix, count matrix being genes as rows and cells as columns.

And that could be a good idea to show the first few rows and columns of your `SPLIT-Seq` output so that we can have and idea.

Regarding gene identifiers, I would stick with gene symbols, to help you in cluster annotation but also using the "MT" prefix for mitochondrial genes for example.

It looks like you are providing somewhat obscure identifiers as gene identifiers, a mix of genome build, ensemble gene id and gene symbol and apparently `Seurat` does not appreciate the `,`. While I cannot give you a definitive recipe on what can and cannot be set as gene identifiers when using Seurat, I can tell that I have worked both with gene symbols and Ensembl identifiers, separately, and both was fine with Seurat.
> 1 votes

---
Tags: scrnaseq, seurat, splitseq
---


# Question
Title: A new paper suggests the Corona Virus has "Uncanny similarity of unique inserts in the 2019-nCoV spike protein to HIV-1" - What does this mean?
Quote:

> We found 4 insertions in the spike glycoprotein (S) which are unique to the 2019-nCoV and are not present in other coronaviruses. Importantly, amino acid residues in all the 4 inserts have identity or similarity to those in the HIV-1 gp120 or HIV-1 Gag. Interestingly, despite the inserts being discontinuous on the primary amino acid sequence, 3D-modelling of the 2019-nCoV suggests that they converge to constitute the receptor binding site. The finding of **4 unique inserts** **in the 2019-nCoV, all of which have identity/similarity to amino acid residues in key structural proteins of HIV-1 is unlikely to be fortuitous in nature**.

Source: https://www.biorxiv.org/content/10.1101/2020.01.30.927871v1

I know this is a preprint paper and it's not yet peer reviewed, but can someone tell me what the implications are if this is true? Does this mean that the virus is artificially constructed?

# Answer
> 11 votes
---

**UPDATE:** The article has now been withdrawn with the following note:

> This paper has been withdrawn by its authors. They intend to revise it in response to comments received from the research community on their technical approach and their interpretation of the results. If you have any questions, please contact the corresponding author.

---

This is very odd, and will require a rigorous investigation, but my initial reaction is one of scepticism. Considering just the 1st insert, the insert sequence is `GTNGTKR`, short at just 7 amino acids. A simple BLASTP vs NR did not find perfect matches to HIV sequences, but did reveal 100% identity across the full sequence length to \>50 other short protein sequences, which could be spurious 'chance' hits of course. Many of the top 100 hits are against eukaryotic sequences, for example one is against *Pristionchus pacificus*, a type of nematode worm (see alignment below). Someone needs to do a proper peer review of this preprint before any conclusions are drawn from it.

```
>tank-1 [Pristionchus pacificus]
Sequence ID: PDM74036.1 Length: 2481
Range 1: 1474 to 1480

Score:24.0 bits(49), Expect:2477,
Method:,
Identities:7/7(100%), Positives:7/7(100%), Gaps:0/7(0%)

Query  1     GTNGTKR  7
             GTNGTKR
Sbjct  1474  GTNGTKR  1480

```

---

**UPDATE 1:** Things are moving fast. There are now 10 comments under the preprint, all of them critical of the idea that these inserts are meaningfully similar to HIV-1 sequences for the reason I outlined above (i.e. short sequences with many hits against organisms from across the tree of life) and additional points, such as small insertions being quite normal evolution for RNA viruses. There are also similar critiques on twitter and biorxiv has added the following disclaimer header to their website:

> bioRxiv is receiving many new papers on coronavirus 2019-nCoV. A reminder: these are preliminary reports that have not been peer-reviewed. They should not be regarded as conclusive, guide clinical practice/health-related behavior, or be reported in news media as established information.

I have never seen such rapid (and unanimous) post-publication peer review before. It is therefore already (just 1 day after publication) fairly clear that at least some claims in this study are completely false and there is no good evidence to support the claim that the 2019-nCoV have acquired sequences from HIV or evidence to suggest the virus is engineered.

# Answer
> 7 votes
Normally "inserts" used in the manuscript are "indels" in protein alignments, short for insertions and deletions.

What I think has happened is a group investigating indels in HIV env noticed indels in 2019-nCov. Essentially I think the correlation is spurious - but I haven't test it, but the area of research in understanding indels is certainly valid and important.

What is certain is that indels induce a large structural change to a protein structure and any Gibbs free-energy style calculation will identify this.

**Vaccine** The spike protein will be the primary candidate to make a 2019-nCov vaccine and that is a very important reason why the sequence was rapidly released. So it is an important protein and the structural changes indels induce mean that a SARS vaccine will probably not provide much protection against 2019-Cov, even apart from the amino acid divergence (below).

**Differences 2019-nCov vs HIV** In summary, alot. HIV env and particularly HIV gag are very different from coronaviruses, both in the mechanism of genome replication, coronavirus never leaves the cytoplasm, clinical outcomes, tissue tropism and duration of infection.

**Similarities** HIV env and the glycoprotein spike of coronaviruses are the receptor binding protein to gain entry into a cell. They are called **structural** proteins. Entry to a cell can be blocked by antibodies and these antibodies are called "neutralizing antibodies". Neutralizing antibodies are catastrophic for a virus. Other antibody responses can be effective, such as IgM, but to clear an infection just using antibodies, you need neutralising antibodies. Both HIV env and the coronavirus spike protein are subject to neutralising antibodies. HIV gag has nothing to do with HIV env, in terms of function or antibody exposure. This is why the spike protein will be the primary vaccine candidate for a subunit vaccine.

**Coincidence, law of chance** There is large variation of indels in HIV env within HIV and what the authors are inferring is there is a resemblence to that between SARS and 2019-nCov. In my opinion this is a coincidence, because they are comparing a large repertoire of HIV varients against a single indel pattern in the coronaviruses.

**Why coronavirus indels?**

That is a very good question. Generically indels in viral surface antigen genes are common, much more common in other proteins - such as those involved in virus replication (non-structural proteins). The amino acid identity between SARS and 2019-nCov is 80%, and in any virus, such as flaviviruses 80% identity means indels will be present in surface antigens between the viruses. The answer is it is not unusual in any RNA virus to see indels at a comparatively large amino acid divergence.

**What function could they serve**

I've briefly looked at indel bioinformatics between flaviviruses (Zika virus, yellow fever virus etc..) notably using envelope (E) protein sequences, and they also occur between African Zika viruses in the E-protein. E-protein being the equivalent of coronavirus spike protein, the receptor-binding protein. No-one has ascribed a function to them and that is the problem with this manuscript.

Hypthoses

* One theory is that a structural change in the protein will occur to stop antibody binding.
* Another theory is they have functional differences, such as cell tropism

Bioinformatically separating the two theories is extremely hard without wetlab experimentation.

# Answer
> 0 votes
SurpriseDog - I suspect the wording "unlikely to be fortuitous in nature" is what led you to ask "Does this mean the virus is artificially constructed?" - No!

I think the writing here is imprecise and confusing. There is no implication that these are human-generated (non-natural/engineered) insertions but rather I think they are implying that these mutations are likely to have functional implications. Whether there is any relation to HIV is questionable.

There has been some political finger pointing and threats of full-scale investigation into whether the virus escaped from a lab. Hard to say, but current thinking is that SARS-Cov-2 is a natural isolate, not engineered (see https://www.nature.com/articles/s41591-020-0820-9). It was probably transmitted from bats.

---
Tags: sequence-alignment, phylogenetics, ngs, phylogeny, covid-19
---


title: Books on bioinformatics algorithms

text:
# Question
Title: Books on bioinformatics algorithms
I'm looking for a book about bioinformatics algorithms, such as alignment, BLAST search, and variant calling.

I'm hoping reading about this subject will give me a deeper understanding of the foundations of bioinformatics, and I'm also interested in developing my own tools and algorithms for sequence and variant analysis. I studied applied mathematics and computer science at university, so I don't shy away from math-heavy books.

I've listed some candidates below.

Durbin et al. is quite old, Sung is the newest.

I am looking for books with:

* material (algorithms) on alignment, assembly, and variant calling if possible.
* New books are preferred

# Answer
> 4 votes
Here are a couple books I'd recommend:

Dan Gusfield's Algorithms on Strings, Trees, and Sequences is a deep and wide treatment of aligning, searching, and processing strings, trees, and sequences.

Warren Ewens' Statistical Methods in Bioinformatics devotes a chapter to BLAST and the math underneath it.

*Edit* \- Another book that may be useful for some mathematical underpinnings would be Michael S. Waterman's Introduction to Computational Biology.

# Answer
> 3 votes
If you're interested in variant calling, I just co-authored an O'Reilly book called Genomics in the Cloud that covers the GATK Best Practices from a scientific standpoint (germline and somatic short variants + somatic copy number alterations), as well as technical considerations like how to run it efficiently at scale (that's where the cloud part comes in).

You can read more about it on the GATK blog; there's an announcement and a second more personal blog post that I wrote to explain how I approached the topic. Happy to discuss more if you're interested.

Oh and you can preview the first few chapters for free on Amazon if you like.

# Answer
> 2 votes
Bioinformatics and Functional Genomics by Jonathan Pevsner is a good one, he goes over the biology as well as how the algorithms work, and provides real-world examples.

---
Tags: sequence-alignment, variant-calling, sequence-analysis, algorithms, literature-search
---


# Question
Title: Generate VCF from different .bam files with different chromosome names
I have two resources of .bam files. One is generated by our lab (1 sample = 1 bam). One is downloaded online (again 1 sample = 1 bam). For the downloaded samples the chromosomes are labelled: chr1, chr2, chr3 etc For our lab samples, the chromosomes are labelled: 1, 2, 3 etc. I want to generate a single VCF file of variants across all samples.

I'm using bcftools:

```
bcftools mpileup -Ov -f ref.fasta -b samples.txt | bcftools call -mv -o bamMge.vcf

```

However, I get no calls and the repeated error: **\[E::faidx\_adjust\_position\] The sequence "1" was not found**

I have two questions:

1. Is my strategy correct (i.e. bcftools mpileup)? Do I need to incorporate extra steps (NB I also have matching gVCF files)
2. How can I either (i) alter the chromosome labeling of one of the subsets of bam files, or, (ii) use a mapping file to match chromosome labels during the mpileup run?

Thanks

# Answer
The main issue here is that your bam files have different chromosome labels. (i.e. 1,2,3 vs. chr1,chr2,chr3) as you mentioned.

This suggests that the data was mapped to a slightly different version of the reference genome. As a first layer, you might want to double check that these reference chromosomes are cross compatible (e.g. both correspond to hg19 for instance). To find more information on the reference genome used, you can check the header information using "samtools view -H "

Assuming that the reference genomes are cross compatible, you can try converting it to sam, adding the missing "chr" label in a pipe, and converting it back to a bam file. However, this will also require modification of the header file which can be challenging in itself.

The easiest (and more rigorous) approach is probably just to convert the batch of bam files with the wrong chromosome labels back into fastq files. You will then want to remap them using the exact same reference genome as the other batch of bam files. After which, you should be able to perform variant calling on the combined batch of bam files.
> 3 votes

# Answer
I am not familiar with getting variants from bam files, so I will just answer the second part. From the error message, it seems that it is necessary to have consistent chromosome labels across datasets. It is also a good idea to do so because this keeps your datasets consistent and makes things a lot easier. To change the chromosomes labels, you can write your own scripts (I use Python, for example). You can also find several solutions on BioStars (https://www.biostars.org/p/98582/).
> 1 votes

---
Tags: vcf, ngs, mpileup, bcftools
---


title: Phasing problem using "samtools phase"

text:
# Question
Title: Phasing problem using "samtools phase"
so I'm having a problem using samtools phase. I'm trying to phase my .bam archives that where generated by aligning my samples with my reference genome. Following the documentation I tried the command:

`samtools phase -A my_referenceGenome.fasta my_sample_toBePhased.bam`

The -A flag stands for dropping reads with ambiguous phase. The **error** I get is: `[W::sam_read1] Parse error at line 1`. I generated .bam files with and without the SAM header but the error persists.

I feel like I am missing something basic but that isn't clearly stated anywere. Does the `my_referenceGenome.fasta` archive needs to be previously phased? If so, how can I do it?

Thank you so much!

# Answer
This is what I see in my samtools phase (v1.10)

```
Usage:   samtools phase [options] <in.bam>

```

samtools phase only takes a bam file as an input. But you are specifying a .fasta file as an input. If you want to specify the reference genome, you should be using the argument "--reference" as described here

```
      --reference FILE
           Reference sequence FASTA FILE [null]

```
> 2 votes

---
Tags: sam, fasta, samtools, phase
---


# Question
Title: What is the perferred method of optimization or energy minimization of small molecules downloaded from PubChem?
I will be docking a small library of molecules that I have downloaded from **PubChem** using **AutoDock**. I thought why not **minimize their energy** before docking using **ChemDraw**. I am thinking; what is the best method of **optimization/ energy minimization** of small molecules. **MM2**, **MMF94** or other? And why?

# Answer
> 2 votes
**The difference between different force-fields is not going to be major, it is the side steps which are.**

## When

If you are starting from a SMILES string, optimisation is a must obviously.

If you are using a 3D conformer from PubChem or even an actual sub-1 Å crystal structure from CSD, optimisation is nice for consistency.

## Which

**MMF94 is a solid choice.** RDkit offers MMF94 (Merck force field 1994) and UFF (universal force field). The latter has some extra features, such as not kicking the bucket at atoms that aren't your usual suspects. ChemDraw offers MM2 and MMF94 IIRC, the former is very old. g16 in Gaussian (paid software) uses Gaussian Electrostatic Model, which is polarisable and will solve partial charges in a nice way (which I'll get to in a minute). OpenBabel uses MMF94 under the hood. You can even optimise with macromolecular forcefields, such as AMBER and CHARMM.

Given that the AutoDock uses a flexible ligand it should not matter too much anyway.

## What else?

Autodock does not require conformers by virtue of the flexible ligand business (the sampler does rototorsions). Other programs might if they use rigid ligands (e.g. Fred) or they use a mixed approach and use a conformer library for speed (e.g. Rosetta ligand\_dock). Conformer diversity is actually a bigger issue as conformers that are only slightly contorted, often even less than kB (0.6 kcal/mol), may not be sampled.

The thing that messes things is **incorrect protonation**. OpenBabel can protonate at a given pH, RDKit cannot. Not an issue for ligands with drug screens (cf. Lipinski rule), but is for natural substrates. However, for the protein it is essential to make sure acidic and basic residue are correctly protonated and that histidine have the correct tautomer —there are tools, c.f. MOE, that can do it but by far the best approach is reading up on the mechanism.

The next thing is **dodgy Gasteiger Marsili partial charges** —for Autodock 4 (and most other full/hybrid physics-based docking programs), partial charges are important: after all PDBQT has a Q in it for charge! For `prepare_ligand4.py` you give it a mol2 file. A major consideration is what format you are using for your molecules. This sounds a very silly thing to point out, but you'd be surprised how many trip up. An mol file (or an sdf file, which is the same but multi-entry) does not have an official partial charge column (although they are often added as comments after `M END`). So whatever force-field you chose make sure that the program you use gave you partial charges are not formal charges. This does not quite apply to AutoDock Vina, which is charge independent and more reliant on machine learning derived factors.

## Verdict

Open babel is nice in that it can start with a SMILES, handles multiple formats including mol2 and PDBQT, can protonate at pH 7, uses MMF94 and does Gasteiger Marsili charges (and its confab algorithm is very good, not need here though).

---
Tags: protein-structure, docking, software-usage, 3d-structure, structural-biology
---


# Question
Title: how to use list of gene id to get cds sequence(cds fasta file have many annotation, only gene id: is same to query id)
i have a question when i want to extract cds sequence using gene id. but cds file is not just start with \>gene is, it has many other annotation. the only same is star with gene:

cds fasta:

```
>Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding
ATCGGCAT
>Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding
ATGCGGCA

```

gene\_list

```
Zm1d035916
Zm1d048

```

how to get result like

```
>Zm1d035916
ATCGGCAT
>Zm1d048
ATGCGGCA

```

# Answer
> 2 votes
If you can use `awk`, you could use this simple script:

```
awk '
BEGIN {p=0}
NR==FNR {gid[$1];next}
/^>/{
  p=0;
  for(g in gid){
    re="\\sgene:"g"\\s"
    if($0~re){p=1;break}
  }
}p' gene_list cds.fasta

```

You can also put the code inside a script file:

extract\_genes.awk

```
#!/usr/bin/awk
BEGIN {
  p = 0; # p indicates whether to print a line or not
}
# while scanning the first file, store gene ids in an array
NR==FNR {
  gid[$1];
  next
}
# check FASTA headers
/^>/ {
  p = 0;
  for(g in gid) {
    # construct regular expression to match gene id in FASTA header
    re = "\\sgene:"g"\\s";
    if ($0~re) { # a match was found!
      p = 1;
      break;
    }
  }
}
p # print line or not?

```

And call it like this:

```
awk -f extract_genes.awk gene_list cds.fasta

```

# Answer
> 1 votes
If you read it in python, they would take care of that:

```
from Bio import SeqIO
records = list(SeqIO.parse("test.fa","fasta"))

records

[SeqRecord(seq=Seq('ATCGGCAT', SingleLetterAlphabet()), id='Zm002', name='Zm002', description='Zm002 cds gene:Zm1d035916 gene_biotype:protein_coding', dbxrefs=[]),
 SeqRecord(seq=Seq('ATGCGGCA', SingleLetterAlphabet()), id='Zm001', name='Zm001', description='Zm001 cds RefGen_v4:9:153880862:153883850:-1 gene:Zm1d048 gene_biotype:protein_coding', dbxrefs=[])]

```

Get the gene and sequence.. maybe there's a easier to do this in biopython, this is like my 10th time using it:

```
def format(rec):
    des = rec.description.split()
    gene = [entry[5:] for entry in des if "gene:" in entry][0]
    return(">"+gene+"\n"+str(rec.seq))

[format(rec) for rec in records]

['>Zm1d035916\nATCGGCAT', '>Zm1d048\nATGCGGCA']

```

Write it out:

```
f=open('example.fasta','w')
f.writelines("\n".join([format(rec) for rec in records]))
f.close()

list(SeqIO.parse("example.fasta","fasta"))
[SeqRecord(seq=Seq('ATCGGCAT', SingleLetterAlphabet()), id='Zm1d035916', name='Zm1d035916', description='Zm1d035916', dbxrefs=[]),
 SeqRecord(seq=Seq('ATGCGGCA', SingleLetterAlphabet()), id='Zm1d048', name='Zm1d048', description='Zm1d048', dbxrefs=[])]

```

---
Tags: python, genome, sequence-analysis, perl
---


# Question
Title: Extract reads from bam files by their @RG
How could I extract reads from bam files by their read groups @RG ?, I've got one file containing all reads for 5 samples, the SM is appears NONE, So, I want to extract each sample reads by the @RG tags.

# Answer
> 4 votes
Typically I use samtools for operations like this. Specifically I use `samtools view` with either `-r` or `-R` flag depending on the use case.

```
-r STR
Output alignments in read group STR [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release.

-R FILE
Output alignments in read groups listed in FILE [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release.

```

# Answer
> 2 votes
This sounds more like a job for `samtools split` if you want to split out all the read groups into separate bams in one go.

http://www.htslib.org/doc/samtools-split.html

```
samtools split [options] merged.sam|merged.bam|merged.cram

Splits a file by read group.

Options:

-u FILE1
Put reads with no RG tag or an unrecognised RG tag into FILE1

-u FILE1:FILE2
As above, but assigns an RG tag as given in the header of FILE2

-f STRING
Output filename format string (see below) ["%*_%#.%."]

-v
Verbose output

```

Alternatively if you do want to proceed one read group at a time, in addition to `samtools view`, you can also do this with GATK using the PrintReads tool with the ReadGroupReadFilter as documented here:

https://gatk.broadinstitute.org/hc/en-us/articles/360041415172-ReadGroupReadFilter

The command line would look like this:

```
gatk PrintReads \
   -I input.bam \
   --read-filter ReadGroupReadFilter \
   --keep-read-group <RG:readgroup> \
   -O output.bam

```

There's also a filter called ReadGroupBlackListReadFilter that allows you to exclude a list of readgroups.

https://gatk.broadinstitute.org/hc/en-us/articles/360041851151-ReadGroupBlackListReadFilter

The PrintReads tool doc is here:

https://gatk.broadinstitute.org/hc/en-us/articles/360041849911-PrintReads

---
Tags: sam, ngs, reads, genome-sequencing, bacteria
---


# Question
Title: How to subset a VCF by chromosome and keep the header?
I would like to subset a VCF which only has chromosome 2.

The problem with using various `grep` commands, e.g.

```
grep -w '^#\|^2' my.vcf > my_new.vcf

```

or if there's a 'chr' prefix

```
grep -w '^#\|^chr2' my.vcf > my_new.vcf

```

is that this will remove the header.

I've been trying to use `vcftools` with the following command:

```
vcftools --vcf input.vcf  --chr 2 --out subset

```

which does not output a VCF `subset.vcf` as expected, but rather `subset.log`:

```
VCFtools - 0.1.15
(C) Adam Auton and Anthony Marcketta 2009

Parameters as interpreted:
    --vcf input.vcf
    --chr 2
    --out subset

After filtering, kept 2 out of 2 Individuals
After filtering, kept 80 out of a possible 720 Sites
Run Time = 0.00 seconds

```

Given that the run time is `Run Time = 0.00 seconds`, I guess there's an error here.

How can I subset a VCF and keep the header? Is there an option with bcftools instead maybe?

# Answer
This seems to work for me with `bcftools filter` and the `-r` or `-R` argument.

> ```
>    -r, --regions chr|chr:pos|chr:from-to|chr:from-[,...]
>
> ```
>
> Comma-separated list of regions, see also **-R**, **--regions-file**. Note that **-r** cannot be used in combination with **-R**.
>
> ```
>    -R, --regions-file FILE
>
> ```
>
> Regions can be specified either on command line or in a VCF, BED, or tab-delimited file (the default). The columns of the tab-delimited file are: CHROM, POS, and, optionally, POS\_TO, where positions are 1-based and inclusive. Uncompressed files are stored in memory, while bgzip-compressed and tabix-indexed region files are streamed. Note that sequence names must match exactly, "chr20" is not the same as "20". Also note that chromosome ordering in FILE will be respected, the VCF will be processed in the order in which chromosomes first appear in FILE. However, within chromosomes, the VCF will always be processed in ascending genomic coordinate order no matter what order they appear in FILE. Note that overlapping regions in FILE can result in duplicated out of order positions in the output. This option requires indexed VCF/BCF files. Note that **-R** cannot be used in combination with **-r**.

Here's an example:

```
$ bcftools filter vcf_nocomp_merge_geno98.vcf.gz -r 4 | head -n 38 | colrm 100 1000000
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20171222
##source=PLINKv1.90
##contig=<ID=1,length=249167691>
##contig=<ID=2,length=242695901>
##contig=<ID=3,length=197800245>
##contig=<ID=4,length=190915651>
##contig=<ID=5,length=180666277>
##contig=<ID=6,length=170877445>
##contig=<ID=7,length=159086805>
##contig=<ID=8,length=146293415>
##contig=<ID=9,length=141018424>
##contig=<ID=10,length=135434552>
##contig=<ID=11,length=134938471>
##contig=<ID=12,length=133763353>
##contig=<ID=13,length=115045730>
##contig=<ID=14,length=107285438>
##contig=<ID=15,length=102369712>
##contig=<ID=16,length=90141356>
##contig=<ID=17,length=81006630>
##contig=<ID=18,length=78014583>
##contig=<ID=19,length=59071322>
##contig=<ID=20,length=62906515>
##contig=<ID=21,length=48077813>
##contig=<ID=22,length=51156934>
##contig=<ID=23,length=154847490>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on rea
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_filterVersion=1.2-187-g1a55e45+htslib-1.2.1-256-ga356746
##bcftools_filterCommand=filter -r 4 vcf_nocomp_merge_geno98.vcf.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  1_125   3_5555  4_7565  5_9
4   71566   rs13125929  T   C   .   PASS    PR  GT  0/1 0/1 0/0
4   87113   rs2006748   T   C   .   PASS    PR  GT  0/1 0/0 0/0
4   110646  rs11727494  C   T   .   PASS    PR  GT  0/1 0/0 1/1
4   142550  rs11735742  T   C   .   PASS    PR  GT  0/1 0/0 0/0
4   200631  rs6826124   A   G   .   PASS    PR  GT  0/1 1/1 0/1
4   221623  rs7695945   G   T   .   PASS    PR  GT  0/0 0/0 0/0

```
> 11 votes

# Answer
You can use `grep` for this. The problem is the `-w` flag which means "match entire words", but VCF header lines start with `##`, not `#`. Here are some alternative `grep` commands:

```
grep -E '^(#|2[[:space:]])' my.vcf

```

Or, if you have GNU `grep` (default on Linux):

```
grep -P '^(#|2\t)' my.vcf

```

To match either `2` or `chr2`:

```
grep -E '^(#|(chr2|2)[[:space:]])' my.vcf

```

or

```
grep -P '^(#|(chr){0,1}2\t)' my.vcf

```

Alternatively, you can use `awk`:

```
awk '/^#/ || $1=="2" || $2=="chr2"' my.vcf

```

Finally, if you want to use `vcftools`, you need to also use the `--recode` option:

> --recode
>
> --recode-bcf
>
> These options are used to generate a new file in either VCF or BCF from the input VCF or BCF file after applying the filtering options specified by the user. The output file has the suffix ".recode.vcf" or ".recode.bcf". By default, the INFO fields are removed from the output file, as the INFO values may be invalidated by the recoding

(e.g. the total depth may need to be recalculated if individuals are removed). This behavior may be overriden by the following options. By default, BCF files are written out as BGZF compressed files.

But that will remove the INFO fields as explained above, so just use normal text-parsing tools like the ones I suggested.
> 4 votes

# Answer
You can do this more easily with the GATK tool SelectVariants, which allows you to subset variants based on various criteria. Here's the doc page:

https://gatk.broadinstitute.org/hc/en-us/articles/360041849831-SelectVariants

For this case (subsetting by chromosome) you would simply do:

```
gatk SelectVariants \
     -R reference.fasta \
     -V input.vcf \
     -L chr2 \
     -O output.chr2.vcf

```

If you wanted to do this for several chromosomes, or regions within chromosomes, you would give the tool a file of intervals with `-L` instead of the chromosome name. Suitable interval file formats are described here:

https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists
> 3 votes

---
Tags: vcf, vcftools
---


# Question
Title: How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?
*This is a problem I have to solve frequently, and I'd be interested in knowing what other methods people use to solve the same problem.*

About twice a year, I get asked to determine variants from Illumina reads, usually from either mouse or human. These are things that have good reference genomes (e.g. human genome from ensembl, mouse genome from ensembl), and are frequently tested out on various bioinformatics tools.

I have at my disposal three computers, in order of usage preference: a laptop with 4 processing threads and 8GB of memory, a desktop with 12 processing threads and 64GB of memory, and a server with 24 processing threads and 256GB of memory. Occasionally I get access to better computers when I need to do lots of processing, or to do things quickly, but these comprise my basic bioinformatics team.

I've had whole-genome projects that involved processing data from 4 individuals, up to projects with about 120 individuals, with each individual having about 10-40M paired-end reads (100-125bp each). Beyond this level, I suspect other researchers will typically have their own bioinformatics team and/or capabilities, which is why I don't get any of the larger-scale jobs that I hear about from salaried bioinformaticians.

What is the standard approach that you would use to process these reads into reference-mapped variant data (i.e. gzipped VCF) for use in downstream analysis tools (e.g. VEP)?

# Answer
> 9 votes
## 1. Adapter Trimming

One of the first things I do after encountering a set of reads is to remove the adapter sequences from the start and end of reads. Most basecalling software includes some amount of built-in adapter trimming, but it is almost always the case that some adapter sequence will remain. Removing adapters is helpful for mapping because it reduces the effort the mapper needs to go through to make a match (which may help some borderline reads to be mapped properly).

My preferred adapter trimmer is Trimmomatic. It has the ability to search for palindromic adapter read through, and includes an adaptive sliding window option. Trimmomatic is threaded, and seems to work reasonably fast when run through files one at a time in threaded mode. Here's an example script that I might run to trim reads from some samples. This puts trimmed FASTQ files into a `trimmed` subdirectory:

```
TRIMPATH=/data/all/programs/trimmomatic/Trimmomatic-0.36;
JARPATH="${TRIMPATH}/trimmomatic-0.36.jar";
ADAPTfile="${TRIMPATH}/adapters/TruSeq3-PE-2.fa";
CMDS=" ILLUMINACLIP:${ADAPTfile}:2:30:10:5:true LEADING:3 TRAILING:3  SLIDINGWINDOW:10:20 MINLEN:40 ";

mkdir -p trimmed;

for INFILE1 in *_1.fastq.gz; do
  base=$(basename "$INFILE1" _1.fastq.gz);
  echo ${base};
  INFILE2="${base}_2.fastq.gz";
  OUTFILE_P1="trimmed/${base}_P1.fastq.gz";
  OUTFILE_P2="trimmed/${base}_P2.fastq.gz";
  OUTFILE_U1="trimmed/${base}_U1.fastq.gz";
  OUTFILE_U2="trimmed/${base}_U2.fastq.gz";
  java -jar "${JARPATH}" PE -threads 20 -phred33 \
        "${INFILE1}"  "${INFILE2}" \
        "${OUTFILE_P1}" "${OUTFILE_U1}" "${OUTFILE_P2}"  "${OUTFILE_U2}" \
        "${CMDS}";
done;

```

Trimmomatic is somewhat sensitive to bad input data, and will loudly fail (i.e. stop running) if it sees quality strings of different length to sequence strings:

```
Exception in thread "Thread-1" java.lang.RuntimeException: Sequence and quality length don't match: 'TACATGGCCCTGAAATGACTTTCACCCAGGCAACCAGTGCCCCCTGTATAGACACATGCCTTGGGCGCTCCCCACCCTTCCTCGCGTGGCCACACCTCTGT' vs '-AAFFJFJJ7A-FJJJFJFJFJJFFA-FFAF-<A-<FFFJA-<-A7-F<<FFJJAJJJJJJJJJ--<--7-7AJ7<AAJA--J7<ACTGCTGTGGGGCACCCAGCCCCCCAGATAGCCTGGCAGAAGGATGGGGGCACAGACTTCCCAGCTGCACGGGAGAGAC'
        at org.usadellab.trimmomatic.fastq.FastqRecord.<init>(FastqRecord.java:25)
        at org.usadellab.trimmomatic.fastq.FastqParser.parseOne(FastqParser.java:89)
        at org.usadellab.trimmomatic.fastq.FastqParser.next(FastqParser.java:179)
        at org.usadellab.trimmomatic.threading.ParserWorker.run(ParserWorker.java:42)
        at java.lang.Thread.run(Thread.java:745)

```

## 2. Read Mapping

Read mapping finds the most likely location for a read within a target genome. There are extremely fast approximate mappers available, but these don't yet work for variant calling, and an exact mapping approach is necessary. My current preferred mapper is HISAT2. I use this instead of BWA due to the double-read issue, and because of the local variant-aware mapping. HISAT2 is made by the same computing group as Bowtie2 and Tophat2 (JHUCCB), and is their recommended tool for replacing those other programs.

> HISAT2 is a successor to both HISAT and TopHat2. We recommend that HISAT and TopHat2 users switch to HISAT2.

The HISAT2 page includes a link to a genomic index file for the human genome, including SNP variants, which I use when mapping human reads. HISAT2 is also threaded, so I run it on the files one at a time and pipe through samtools to create a sorted BAM file:

```
mkdir mapped
for r1 in trimmed/*_P1.fastq.gz;
  do base=$(basename "${x}" _P1.fastq.gz);
  r2="trimmed/${base}_P2.fastq.gz";
  echo ${base};
  hisat2 -p 20 -t -x /data/all/genomes/hsap/hisat_grch38_snp/genome_snp -1 \
    "${r1}" -2 "${r2}" 2>"mapped/hisat2_${y}_vs_grch38_stderr.txt" | \
    samtools sort > "mapped/hisat2_${y}_vs_grch38.bam";
done

```

## 3. Variant Calling

I currently use samtools/bcftools to do this, but would be interested in other opinions as we recently had a grant reviewer response that samtools was a dinosaur program and better approaches were available. Samtools doesn't currently work with threads for variant calling, so I save the commands to a text file and then run them through GNU Parallel. This requires the genome files to be first downloaded and indexed via `samtools faidx`. This produces a set of gzipped VCF files in the `variants` directory:

```
mkdir -p variants;
(for x in mapped/*.bam;
   do echo samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa "${x}" \| \
   bcftools call --ploidy GRCh38 -v -m -O z -o variants/$(basename "${x}" .bam).vcf.gz; done) > call_jobs.txt
cat call_jobs.txt | parallel -u;

```

As a slow, but more accurate alternative, variants can be called for all samples at the same time:

```
samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa mapped/*.bam | \
  bcftools call --ploidy GRCh38 -v -m -O z -o variants/hisat2_allCalled_vs_grch38.vcf.gz

```

If a faster processing is desired, then the variant calling can be parallelised across chromosomes using `parallel` and merged afterwards using `bcftools norm`. The implementation of this is left as an exercise to the reader.

*Aside: I see that GATK's best practise guidelines call per-sample with HaplotypeCaller, then process those VCF files as aggregated data. @Geraldine-vdauwera says this two-step joint calling workflow was designed to be scalable for any size cohort, so you can apply it pretty much the same way whether you have two samples or twenty thousand. Once you go higher than that there are some additional tweaks to keep scaling up.*

# Answer
> 5 votes
The GATK is probably the most widely used program for this. There are several workflows you can use depending on what types of variants you're interested in (germline or somatic, short variants, copy number, structural). The most common use case is Germline short variant discovery (SNPs and Indels), for which the overall workflow is described here:

https://gatk.broadinstitute.org/hc/en-us/articles/360035535932

The GATK team at the Broad Institute now provides workflow scripts written in WDL (Worfklow Description Language) that can be run on pretty much any platform with the Cromwell engine. The workflows include data pre-processing, including mapping (alignment) and some important cleanup operations, plus variant calling and filtering.

One of the advantages of using those workflows is that the Broad uses them in production for all the data they sequence, which nowadays amounts to abut a new human whole genome every five minutes. That means there's a lot of pressure to make sure the tools are reliable and efficient because any interruption of the pipelining operations would be a big problem. They're also highly standardized so that the outputs will be compatible with other datasets and tools.

Disclosure: I work at the Broad Institute and wrote a book about using GATK (https://oreil.ly/genomics-cloud).

---
Tags: variant-calling, read-mapping, human-genome
---


# Question
Title: Meaning of the FORMAT fields of the VCF file coming from GIAB project
After reading the GIAB paper in https://www.biorxiv.org/content/early/2018/05/25/281006 and its Figure 1, I am still having trouble understanding the data inside the GIAB VCF file for HG001 (HG001\_GRCh38\_GIAB\_highconf\_CG-IllFB-IllGATKHC-Ion-10X-SOLID\_CHROM1-X\_v.3.3.2\_all.vcf). In particular, I need to understand the information that led to a FILTER value of PASS for some variants.

The first case that I would like to discuss is the following variant:

```
chr1    110176467   .   A   G   50  PASS    platforms=3;platformnames=Illumina,CG,10X;datasets=3;datasetnames=HiSeqPE300x,CGnormal,10XChromium;callsets=4;callsetnames=HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo;datasetsmissingcall=IonExome,SolidPE50x50bp,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable    GT:PS:DP:ADALL:AD:GQ    0/1:.:512:114,106:145,148:760

```

Translating the information in the line to English, this is what I get:

* We considered data from the following sequencing technologies: Illumina,CG,10X,IonExome,SolidPE50x50bp, and SolidSE75bp.
* Of those, some of them (IonExome,SolidPE50x50bp,SolidSE75bp) were missing the call, so we discarded them.
* Therefore, the data for the variant is going to come from 3 platforms (Illumina,CG,10X), 3 datasets (HiSeqPE300x,CGnormal,10XChromium), and 4 callsets (HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo), because we analyzed the data from HiSeqPE300x with GATK and also with Freebayes.
* However, the callsets from CS\_HiSeqPE300xGATK,CS\_CGnormal, and CS\_HiSeqPE300xfreebayes had this call in a region with low coverage of high MQ reads.
* Across all the platforms, DP=512
* ADALL (Values 114,106) is coming from “all" the datasets. Are those “all” the set (HiSeqPE300x + Cgnormal + 10XChromium), or the set (Illumina + CG + 10X + IonExome + SolidPE50x50bp + SolidSE75bp).
* AD (Values 145,148) is coming from the same “all” datasets as ADALL, but they are unfiltered, so 145\>114 and 148\>106.

All this make sense, but I would like to know what case of Figure 1c was applied to arrive at the PASS value.

The second variant that I would like to mention is this one:

```
chr1    5705293 .   T   C   50  PASS    platforms=3;platformnames=Illumina,CG,10X;datasets=3;datasetnames=HiSeqPE300x,CGnormal,10XChromium;callsets=4;callsetnames=HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo;datasetsmissingcall=IonExome,SolidPE50x50bp,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable;filt=CS_HiSeqPE300xGATK_filt   GT:PS:DP:ADALL:AD:GQ    0/1:.:599:121,139:49,28:627

```

Again, my translation would be:

* We considered data from the following sequencing technologies: Illumina,CG,10X,IonExome,SolidPE50x50bp, and SolidSE75bp.
* Of those, some of them (IonExome,SolidPE50x50bp,SolidSE75bp) were missing the call, so we discarded them.
* Therefore, the data for the variant is going to come from 3 platforms (Illumina,CG,10X), 3 datasets (HiSeqPE300x,CGnormal,10XChromium), and 4 callsets (HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo), because we analyzed the data from HiSeqPE300x with GATK and also with Freebayes.
* However, the callsets from CS\_HiSeqPE300xGATK,CS\_CGnormal, and CS\_HiSeqPE300xfreebayes had this call in a region with low coverage of high MQ reads.
* Also, the data from CS\_HiSeqPE300xGATK\_filt was filtered/discarded
* Across all the platforms, DP=599
* ADALL (Values 121,139) is coming from “all" the datasets. I assume (Cgnormal + 10XChromium), because CS\_HiSeqPE300xGATK\_filt was filtered.
* AD (Values 49,28) are unfiltered.

But notice that 49\<121 and 28\<139. How come that there are less unfiltered reads (AD)? Should not be the case that AD \> ADALL for all instances?

I don’t know how to reconcile the two variants. Aren’t they providing opposite information in the AD and ADALL fields? Which specific datasets and callsets are involved in getting AD and ADALL for each variant?

I would also appreciate answers to a couple of other questions:

* Why are these two variants a PASS if there are 3 “callable” datasets in regions of low coverage?
* Did the variants follow an arbitration process? I guess not, if the INFO flag “arbitrated” is missing.

For reference, here are the descriptions in the header of the VCF file that I have been using:

```
##INFO=<ID=platforms,Number=1,Type=Integer,Description="Number of different platforms for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=platformnames,Number=.,Type=String,Description="Names of platforms for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=datasets,Number=1,Type=Integer,Description="Number of different datasets for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=datasetnames,Number=.,Type=String,Description="Names of datasets for which at least one callset called this genotype, whether filtered or not">
##FORMAT=<ID=ADALL,Number=R,Type=Integer,Description="Net allele depths across all datasets">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Net allele depths across all unfiltered datasets with called genotype">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Total read depth summed across all datasets, excluding MQ0 reads">
##INFO=<ID=datasetsmissingcall,Number=.,Type=Integer,Description="Names of datasets that are missing a call or have an incorrect call at this location, and the high-confidence call is a variant">
##INFO=<ID=callable,Number=.,Type=String,Description="List of callsets that had this call in a region with low coverage of high MQ reads.">
##INFO=<ID=filt,Number=.,Type=String,Description="List of callsets that had this call filtered.">
##INFO=<ID=arbitrated,Number=1,Type=String,Description="TRUE if callsets had discordant calls so that arbitration was needed.">

```

# Answer
> 2 votes
The filtering criteria should be defined in the header of the VCF file. Can you include the header from your VCF file?

VCF Format and example header:

```
FILTERs that have been applied to the data should be described as follows:

##FILTER=<ID=ID,Description=”description”>

```

header:

```
##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">

```

# Answer
> 0 votes
For your second question, I don't think that variant matches any of the examples shown in Figure 1c. Note that that panel's caption says it shows four examples of arbitration. Per you last sub-question, I don't think this variant was subjected to arbitration.

For your second question, I believe the subtlety is in the field descriptions:

ADALL = Net allele depths across all datasets

AD = Net allele depths across all unfiltered datasets **with called genotype**

My understanding is that ADALL is a count of allele depths that includes samples where the variant was filtered out or the genotype was not concordant; something to that effect. So generally I would expect ADALL \> AD.

I'm honestly not sure exactly what they mean by "unfiltered datasets". That could mean "counting all reads without filtering on quality" or "callsets to which variant filtering was not yet applied". Maybe there's another subtlety there that explains why the two variants have opposite AD vs ADALL proportions.

For your additional questions:

> Why are these two variants a PASS if there are 3 “callable” datasets in regions of low coverage?

According to the paper, having the variant in low-coverage regions isn't disqualifying if the genotype calls are concordant, which makes sense to me. Low coverage makes me more suspicious, but if multiple callsets are in agreement, I'm still willing to consider the call good (unless of course they all exhibit the same error mode, which is the main potential weakness of the consensus approach used by GiaB).

> Did the variants follow an arbitration process? I guess not, if the INFO flag “arbitrated” is missing.

Right, if that flag is not present, it means arbitration was not needed/applied.

---
Tags: vcf, ngs, file-formats, genomics, human-genome
---


# Question
Title: How can I count the number of reads that support a variant in a bam file?
I am calling variants from a human sample using bwa mem to align the reads and gatk to call the variants. I'm trying to understand why a specific variant was not called in my sample. I have checked the bam alignments in a GUI viewer and I can see that there are reads supporting the missing variant. It looks like the issue is a low allelic balance, with far more reads supporting the reference than the alternate allele but I want to get the actual numbers.

So, given a specific variant like this:

```
chr22   425236  C T

```

How can I count the number of reads in my `sample.bam` file that support that variant and the number that don't on Linux?

# Answer
If you don't mind a bit of manual counting, then `samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam` will produce output where you can count the bases for that position. You could, of course, use the command line to do most of that automatically:

```
samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam |
    cut -f 5 | tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c

```

That'll give you a count of how many of each base were seen.
> 8 votes

# Answer
ASCIIGenome (I'm the author) has a command, filterVariantReads, designed to inspect reads having a variant at a position or range. It would go along these lines:

```
ASCIIGenome -fa genome.fa aln.bam

```

Then go to the region of interest and use:

```
goto chr9:4917981-4918161
filterVariantReads -r 4918011

```

From this:

You get:

You can also print to screen or save to file the underlying sam records with:

```
print > reads.sam

```

For several regions, the whole process can be scripted for automation.

Hope this helps!
> 4 votes

# Answer
The variant record in the VCF produced by GATK should include information about depth -- both total read depth (DP) and allelic depth (AD), ie depth per allele.

Keep in mind however that the caller makes its decision based on more than just allelic depth. There are some additional statistics involved, including the quality of the individual bases. Also, if any filtering was applied, your variant may have been excluded due to other statistics, like for example if there is evidence of strand bias -- eg if all the reads supporting your variant are in the same orientation, that is a classic red flag suggesting it's an artifact.
> 0 votes

---
Tags: sam, ngs, variant-calling
---


# Question
Title: Is there a read mapper that allows input in SAM format instead of FASTA/FASTQ?
When we get raw sequenced reads in FASTQ format, we usually need to do some pre-processing (QC, demultiplexing, etc.) prior to sequence alignment. The only way to register the results of this preprocessing is to add metadata to each read. This is best done by using the SAM format which has a well-defined system of metadata tags.

However, it seems that not a single read mapper (except maybe for *bwa aln* with option *-b*) accepts SAM (or one of the binary equivalents BAM/CRAM) as input to make use of the metadata that has been added during preprocessing. This implies that one has to convert the unaligned reads back to FASTQ format.

In any case, whether one uses FASTQ throughout or converts back to FASTQ, a common annoyance with FASTQ is that metadata is registered as part of the query name - this is not standardized or even well-defined. After alignment, we again have to do some text processing to get the metadata out of the query name and into proper metadata tags defined by SAM. This back and forth hacking is highly unfortunate and error prone. I have seen this, for **example**, in some single-cell sequencing projects where we have highly multiplexed sequencing runs with reads that need to be QC'd and classified. The new SAM standard from May 2018 does indeed define proper tags for handling cell-multiplexed runs. SAM is the de-facto standard to store both unaligned and aligned reads for years to come so it should be standard for read mappers to accept it as input.

**I would be grateful for some hints and if we could collect here some read mappers that do accept SAM/BAM/CRAM as input format.**

# Answer
The common solution for scRNA-seq is to put cell barcodes and such in read headers and then post-process things with UMItools.

But regarding your actual question, STAR can accept SAM/BAM as input with the `--readFilesType SAM PE` option (it's `SAM` for both that format and BAM). You can swap `SE` for `PE` if your data is single-end (or is effectively that way, as is the case in many scRNA-seq experiments). Note that in the case of PE experiments that mates must be next to each other. I really don't know what this does with auxiliary tags though, so you would need to see if they're lost before you start processing a lot of data).
> 5 votes

# Answer
With bwa-mem or minimap2, the recommended way is

```
samtools fastq -T BC,RX name-grouped.sam | bwa mem -C -p -

```

This passes the `BC` and `RX` tags through bwa-mem and copies them to the SAM output. If you are running paired-end alignment and your BAM is coordinate sorted, you have to run `samtools collate` to group reads by name:

```
samtools collate -Ol0 in.bam tmp | samtools fastq -T BC,RX - | bwa mem -pC -

```

Here `bwa mem -C` or `minimap2 -y` copies FASTQ comments to SAM output. This strategy is much easier to deal with than SAM/BAM. A caveat is that it doesn't work well for read groups because each read group requires a separate header line. You have to manually modify the header.

I am not sure which short-read RNA-seq mappers have implemented a similar feature. You might consider to send a feature request.
> 5 votes

# Answer
At the Broad Institute, we do something similar to what is described by user172818. We use uBAM (unmapped BAM) to store unmapped reads, because as you say the SAM format has a more appropriate way to encode metadata. However, we use a Picard tool called SamToFastq to do the initial conversion then pipe the output of that to BWA mem. Then we use an additional step to reattach metadata, read group info and so on, using a Picard tool called MergeBamAlignments.

You can see the full implementation here (starts as line 71):

https://github.com/gatk-workflows/gatk4-genome-processing-pipeline/blob/master/tasks/Alignment.wdl
> 1 votes

# Answer
SNAP is a hash-index based genomic aligner that can take SAM and BAM as input. It has been used as an aligner for RNA-seq (see SNAPR) and BS-seq (see BS-Seeker3). SNAPR seems to accept SAM/BAM but it is not clear whether BS-Seeker3 wraps this functionality as well.

I tried SNAP and it seems to do the job well, i.e. it only adds tags during alignment and does not remove metadata that has been deposited prior to alignment, neither from the SAM header nor from individual reads.

**NOTE** I post this answer to my own question after I became aware of SNAP.
> 0 votes

---
Tags: sam, ngs, read-mapping, data-preprocessing
---


# Question
Title: Picard validation error regarding bin field of BAM file
I am sure I set the path right but whatever I am trying the command not working

Any help please?

```
[fi1d18@cyan02 fi1d18]$ picard MarkDuplicates \
    I=DNA.dupmarked.bam \
    O=DNA.dupmarked1.bam \
    M= marked-dup-metrics.txt
[Thu Mar 07 17:33:42 GMT 2019] picard.sam.markduplicates.MarkDuplicates INPUT=[/temp/hgig/fi1d18/1631_WTSI-OESO_005_a_DNA/mapped_sample/HUMAN_1000Genomes_hs37d5_genomic_WTSI-OESO_005_a_DNA.dupmarked.bam] OUTPUT=/temp/hgig/fi1d18/1631_WTSI-OESO_005_a_DNA/mapped_sample/HUMAN_1000Genomes_hs37d5_genomic_WTSI-OESO_005_a_DNA.dupmarked1.bam METRICS_FILE=marked-dup-metrics.txt MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 SORTING_COLLECTION_SIZE_RATIO=0.25 REMOVE_SEQUENCING_DUPLICATES=false TAGGING_POLICY=DontTag REMOVE_DUPLICATES=false ASSUME_SORTED=false DUPLICATE_SCORING_STRATEGY=SUM_OF_BASE_QUALITIES PROGRAM_RECORD_ID=MarkDuplicates PROGRAM_GROUP_NAME=MarkDuplicates READ_NAME_REGEX=<optimized capture of last three ':' separated fields as numeric values> OPTICAL_DUPLICATE_PIXEL_DISTANCE=100 VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json
[Thu Mar 07 17:33:42 GMT 2019] Executing as fi1d18@cyan02 on Linux 2.6.32-754.el6.x86_64 amd64; Java HotSpot(TM) 64-Bit Server VM 1.8.0_51-b16; Picard version: 2.8.3-SNAPSHOT
INFO 2019-03-07 17:33:42 MarkDuplicates Start of doWork freeMemory: 2012347496; totalMemory: 2027945984; maxMemory: 3817865216
INFO 2019-03-07 17:33:42 MarkDuplicates Reading input file and constructing read end information.
INFO 2019-03-07 17:33:42 MarkDuplicates Will retain up to 14684096 data points before spilling to disk.
WARNING: BAM index file /temp/hgig/fi1d18/1631_WTSI-OESO_005_a_DNA/mapped_sample/HUMAN_1000Genomes_hs37d5_genomic_WTSI-OESO_005_a_DNA.dupmarked.bam.bai is older than BAM /temp/hgig/fi1d18/1631_WTSI-OESO_005_a_DNA/mapped_sample/HUMAN_1000Genomes_hs37d5_genomic_WTSI-OESO_005_a_DNA.dupmarked.bam
[Thu Mar 07 17:33:42 GMT 2019] picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 0.01 minutes.
Runtime.totalMemory()=2027945984
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
Exception in thread "main" htsjdk.samtools.SAMFormatException: SAM validation error: ERROR: Record 3752, Read name HX3_22030:3:2114:23155:23319, bin field of BAM record does not equal value computed based on alignment start and end, and length of sequence to which read is aligned
at htsjdk.samtools.SAMUtils.processValidationErrors(SAMUtils.java:448)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.advance(BAMFileReader.java:665)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:650)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:620)
at htsjdk.samtools.SamReader$AssertingIterator.next(SamReader.java:569)
at htsjdk.samtools.SamReader$AssertingIterator.next(SamReader.java:543)
at picard.sam.markduplicates.MarkDuplicates.buildSortedReadEndLists(MarkDuplicates.java:438)
at picard.sam.markduplicates.MarkDuplicates.doWork(MarkDuplicates.java:222)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:205)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:94)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:104)
[fi1d18@cyan02 fi1d18]$

```

The error message is:

> Exception in thread "main" htsjdk.samtools.SAMFormatException: validation error: ERROR: Record 3752, Read name HX3\_22030:3:2114:23155:23319, bin field of BAM record does not equal value computed based on alignment start and end, and length of sequence to which read is aligned
>
> EDITED

```
[fi1d18@cyan02 sambamba]$ samtools view -H /temp/hgig/fi1d18/1631_WTSI-OESO_005_a_DNA/mapped_sample/HUMAN_1000Genomes_hs37d5_genomic_WTSI-OESO_005_a_DNA.dupmarked.bam
[W::bam_hdr_read] EOF marker is absent. The input is probably truncated.
@HD     VN:1.5  SO:coordinate
@SQ     SN:1    LN:249250621    AS:NCBI37       M5:1b22b98cdeb4a9304cb5d48026a85128     SP:Human  UR:/lustre/scratch117/core/sciops_repository/references/Homo_sapiens/1000Genomes_hs37d5/all/fasta/hs37d5.fa
@SQ     SN:2    LN:243199373    AS:NCBI37       M5:a0d9851da00400dec1098a9255ac712e     SP:Human  UR:/lustre/scratch117/core/sciops_repository/references/Homo_sapiens/1000Genomes_h

```

# Answer
> 0 votes
We've seen this error before on the GATK forum; if you google "bin field of BAM record does not equal value computed based on alignment start and end" you will find this forum thread:

https://gatkforums.broadinstitute.org/gatk/discussion/4290/sam-bin-field-error-for-the-gatk-run

This part of the last post should help you solve this:

> The only way I can think of to fix this is to fix the BAM. In order to do that, you would need to convert to SAM format, which does not have a field for indexing bin, and then back to BAM. That would be the safest thing, followed by recreating the BAI. So you can do this with GATK tools using the following commands.
>
> ```
> gatk SamFormatConverter - to convert Bam to Sam
>
> gatk SamFormatConverter - to convert Sam back to Bam
>
> gatk BuildBamIndex - to create the Bam index file.
>
> ```
>
> Also, please make sure your Picardtools and GATK versions are up to date to the most recent versions.

---
Tags: sam, wgs, picard
---


# Question
Title: Help with the definitions of fields in a VCF output by Strelka
In Strelka vcf in INFO column we have these

```
Depth="The number of reads covering the variant position including duplicates, supplementary records and reads that fall below minimum base and mapping quality thresholds."

ReadCount="The number of reads covering the variant position excluding duplicates, supplementary records and reads that fall below minimum base and mapping quality thresholds.">
ReadCountControl="The number of reads covering the variant position in the control excluding duplicates, supplementary records and reads that fall below minimum base and mapping quality thresholds.">

```

And in FORMAT column we have

```
AU="Number of 'A' alleles used in tiers 1,2">
CU="Number of 'C' alleles used in tiers 1,2">
DP="Read depth for tier1 (used+filtered)">
GU="Number of 'G' alleles used in tiers 1,2">
TU="Number of 'T' alleles used in tiers 1,2">

```

What is the difference of **ReadCount** in INFO column and **DP** in FORMAT column?

What is the difference of **number of allele** versus n**umber of reads covering the variant position** ?

Please help for clarification if you have any intuition on this definitions

Thanks

# Answer
> What is the difference of ReadCount in INFO column and DP in FORMAT column?

I'm not super familiar with Strelka but I believe it follows the VCF spec, so DP in the FORMAT column should give the depth of coverage (ie read count) per sample. ReadCount in INFO is not in the VCF spec but based on the definition it sounds like the depth of coverage aggregated over all samples.

> What is the difference of number of allele versus number of reads covering the variant position ?

I believe the "number of alleles used" refers to counts of each base that are used in the indel error modeling process as described in the preprint here:

https://www.biorxiv.org/content/10.1101/192872v2.full

So it sounds like subsets of the number of reads covering the variant position.
> 2 votes

---
Tags: vcf, variant-calling, wgs, snv
---


title: Parsing .vcf file for this information

text:
# Question
Title: Parsing .vcf file for this information
I have a .vcf file

https://www.dropbox.com/s/8v73nppwg3a1tnd/LP2000109-DNA\_A01\_vs\_LP2000103-DNA\_A01.SVannotated.txt?dl=0

with this header

```
##startTime=Fri Mar 29 16:46:32 2019
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  NORMAL  TUMOR
1   54586   .   T   C   .   PASS    DP=39;MQ=50.55;MQ0=0;NT=ref;QSS=48;QSS_NT=48;ReadPosRankSum=1.92;SGT=TT->CT;SNVSB=0.00;SOMATIC;SomaticEVS=10.83;TQSS=1;TQSS_NT=1    AU:CU:DP:FDP:GU:SDP:SUBDP:TU    0,0:0,0:20:0:0,0:0:0:20,20  0,0:6,6:18:0:0,0:0:0:12,13
1   103241  .   C   T   .   PASS    DP=120;MQ=24.94;MQ0=35;NT=ref;QSS=47;QSS_NT=47;ReadPosRankSum=2.09;SGT=CC->CT;SNVSB=0.00;SOMATIC;SomaticEVS=9.44;TQSS=2;TQSS_NT=2   AU:CU:DP:FDP:GU:SDP:SUBDP:TU    0,1:32,47:33:1:0,0:0:0:0,5  0,

```

The "DP" field in the vcf shows the depth of the individual samples; So in this file, the first locus has the following format fields:

```
AU:CU:DP:FDP:GU:SDP:SUBDP:TU    0,0:0,0:20:0:0,0:0:0:20,20  0,0:6,6:18:0:0,0:0:0:12,13

```

So according to this (DP field of normal and tumor samples), normal sample has a depth of 20 and tumor sample has a depth of 18.

So how I could extract the read depth for all loci as described for the first position? The desired output would be like this \[note that the VCF is taken from my own data but the table is my desired format that I don't know how to get that from my own data. chr have been added manually because my reference genome is hg19\]:

```
  Sample   Type   CHROM       POS    REF  ALT       Tumor_Depth Normal_Depth
  CHC2432T SNV    chr1  102961055    G    A                  64           62
  CHC2432T SNV    chr1  105492588    A    T                  66           73
  CHC2432T SNV    chr1  108628724    C    T                  45           54
  CHC2432T SNV    chr1  109692113    G    T                  53           29
  CHC2432T SNV    chr1  109692114    G    T                  53           31
  CHC2432T SNV    chr1  120676701    T    C                  48           87

```

# Answer
> 1 votes
To extract the DP fields from a VCF file, you could use a tool like `bcftools query`:

> Extracts fields from VCF or BCF files and outputs them in user-defined format.

You could start from something like this:

```
bcftools query -Hf 'CHC2432T\t%TYPE\t%CHROM\t%POS\t%REF\t%ALT[\t%DP]\n' file.vcf

```

# Answer
> 1 votes
You can do the extraction part with the GATK tool VariantsToTable, as described here:

https://gatk.broadinstitute.org/hc/en-us/articles/360041414592-VariantsToTable

The usage example from that doc:

```
gatk VariantsToTable \
    -V input.vcf \
    -F CHROM -F POS -F TYPE -GF AD \
    -O output.table

```

would produce a file that looks like:

```
 CHROM  POS        TYPE   HSCX1010N.AD  HSCX1010T.AD
 1      31782997   SNP    77,0          53,4
 1      40125052   SNP    97,0          92,7
 1      65068538   SNP    49,0          35,4
 1      111146235  SNP    69,1          77,4

```

So you might still need to reorder columns etc but that should allow you to at least get the values out in a tabular format that will be easier to work with.

---
Tags: vcf, wgs, snv
---


# Question
Title: What is the meaning of the "*" character in bwa fastmap's output?
I am mapping kmers back to a few bacterial genomes using `bwa fastmap`:

```
bwa fastmap -l 9 ref.fasta kmers.fasta > out.fastmap
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[main] Version: 0.7.17-r1188
[main] CMD: bwa fastmap -l 9 ref.fasta kmers.fasta
[main] Real time: 0.168 sec; CPU: 0.157 sec

```

The output looks like this:

```
SQ      AAAAGTAGTAAGCAGGAAGACAACACGGTTG 31
EM      0       11      1       contig1:+2041368
EM      3       13      2       contig1:+1875252        contig1:-2474744
EM      4       15      1       contig1:+3779779
EM      5       16      1       contig1:-3253348
EM      6       17      2       contig1:+938710 contig1:+1066682
EM      7       20      1       contig1:-4912797
EM      10      21      4       contig1:+4803473        contig1:-2830252        contig1:+4495150        contig1:-4907283
EM      11      23      1       contig1:-3148132
EM      13      24      1       contig1:-4651172
EM      14      25      3       contig1:-2142223        contig1:+4994474        contig1:+873066
EM      16      26      4       contig1:+156775 contig1:+27749  contig1:+2207492        contig1:-1340811
EM      17      29      1       contig1:+3523989
EM      19      30      24      *

EM      20      31      3       contig1:-2533354        contig1:-208660 contig1:-1080177
//

```

My question is: what is the meaning of the `*` char in the row before last? Why is `bwa fastmap` not reporting all the hits as in the other lines? I could not find a page that explains this output anywhere.

# Answer
> 3 votes
Just starting to work with bwa fastmap and indeed the documentation is very sparse... This is what we get from running the program without arguments:

```
bwa fastmap

Usage:   bwa fastmap [options] <idxbase> <in.fq>

Options: -l INT    min SMEM length to output [17]
         -w INT    max interval size to find coordiantes [20]
         -i INT    min SMEM interval size [1]
         -L INT    max MEM length [2147483647]
         -I INT    stop if MEM is longer than -l with a size less than INT [0]

```

Which is a bit cryptic as is...

The output does not really follow SAM specification and I'm pretty sure the '\*' means does not mean unmapped (in fact an unmapped Exact Match doesn't make any sense does it?)

In fact, the '\*' means that there are too many matches to print the reference id and position. Looking quickly into the source code this is controlled by the -w arguments.

It defaults to 20, which is why you don't get an output for your 24 matches.

# Answer
> 3 votes
It means the read is unmapped. No origin in the reference could be found that was sufficiently similar to the read to call it a proper alignment.

In the SAM specification, `*` is always used if the information for that field is not available. Therefore, if no mapping information is available (=unmapped) one sets `*`.

---
Tags: sequence-alignment, bwa
---


# Question
Title: How to downsample some of the samples in RNA-seq data?
I have 40 samples and these are into two groups. I would like to perform a differential analysis between two groups. The library size of the samples is very low. But there are two samples (GroupA\_12 and GroupB\_2) in which the number of counts is too large compared to other samples. Please check the library sizes below:

I wanted to downsample these two cases and then would like to perform differential analysis. May I know how to downsample those two samples.

Any help is appreciated. thanq.

# Answer
> 4 votes
You have a few options:

1. Downsample the fastq files and rerun the entire analysis. You can do this with `seqtk sample`.
2. Downsample the BAM files, which you can do with `samtools view -s`.
3. Divide all of the counts in the counts files by some factor and round that to an integer.

I personally prefer option 2, since it's quick and doesn't usually have any side-effects unless you use two-pass alignment or did the quantification with something like salmon. Option 3 is always tempting, but I have to say I've never been convinced that this consistently produces similar enough results to options 1 or 2.

---
Tags: r, rna-seq, differential-expression, gene-expression
---


# Question
Title: Can I automate CLUSTALO and output alignment sequence identity?
I've detected homology between targets of ligands in drugbank and proteins in the proteome of a pathogen. I've parsed the output very rudimentary and calculated my query coverage. This exists in an excel file (see attached)

I'd like to further filter these hits. For query coverages \> 70% alone I have 1500+ hits. Some redundancy exists and some hits are for drugbank targets for which no data has been deposited in the PDB.

So there are ways of further filtering these hits, but I'd like to do pairwise alignments first before I do anything else, as I think this is the most rational filter I can think of.

I do not want to filter based on lack of PDB data necessarily, as when I check my high query coverage hits, I'm getting predictions that are actually backed up by in vitro data. This includes high coverage hits lacking x-ray data

My question:

I have CLUSTALO installed on my desktop: How do programmatically send it jobs if I have the Drugbank sequences and the pathogen sequences in the same file? I would also like to parse out the pathogen IDs from the excel file. There’s 1500 and they are a bit scattered throughout the columns.

```
0.980634         -                    572         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-263             874.5          0.0            1        1      6.6e-267  1.4e-263   874.2      0.0       5         562    15     570    12     572    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005678  PE=3             SV=1

0.950704         -                    580         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-258             858.0          0.0            1        1      6.8e-262  1.5e-258   857.7      0.0       6         546    10     560    7      578    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005014  PE=3             SV=1

0.950704         -                    691         sp|P06672|PDC_ZYMMO                         -                     568                               8.2e-162             538.3          0.1            1        1      2.4e-164  5.1e-161   535.6      0.1       1         541    89     650    89     658    0.89   Acetolactate                            synthase                  OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_006830  PE=3             SV=1

0.492958         -                    585         sp|P06672|PDC_ZYMMO                         -                     568                               8e-87                290.7          0.0            1        2      3.3e-83   7e-80      267.8      0.0       3         283    2      287    1      292    0.94   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_001897  PE=3             SV=1

0.107394         -                    585         sp|P06672|PDC_ZYMMO                         -                     568                               8e-87                290.7          0.0            2        2      1.7e-08   3.6e-05    21.1       0.0       413       474    436    505    376    583    0.72   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_001897  PE=3             SV=1

```

# Answer
So, as the first reply said, there's sort of multiple questions here.

For filtering things out from a table, I would absolutely use pandas. https://swcarpentry.github.io/python-novice-gapminder/08-data-frames/ Pandas is a way of managing dataframes in python. It looks like what you want from your table is to remove certain rows, if the content of those rows at a certain column matches certain IDs from a list (pathogen IDs). That's the sort of thing the pandas library is good at.

For managing workflows, I use either bash or snakemake. Snakemake is a python thing. It has some advantages, such as being able to start from an intermediary file (for example, if a workflow is "sequences -\> multiple alignment -\> phylogenetic tree", it can easily start from the multiple alignment). Also deals with branched workflows quite well. Can be easily installed with conda (anaconda is absolutely a lifesaver for bioinformatics work imo). I still use bash for simple workflows, because it's slightly simpler/faster.

For potential filtering - If you want to keep the top hits as they are, but filter the lower-scoring hits based on whether they have a PDB ID or not, I'd personally look into Entrez queries. They're a way of making your search more specific. Here they're shown for WebBlast, but they work quite well for command line blast too. https://www.youtube.com/watch?v=bxx5uaKjMa8
> 1 votes

# Answer
The question is still in numerous parts,

* parsing a blast output
* obtaining sequences
* feeding this as input into clustalo

If I was automating this I would parse the Blast genbank codes \[regex either Perl one-liner (scripts aren't trendy ;-( ), or Pythons "re"\] and feed them into NCBI's `efetch` to hook out the sequence data. Then I would create a subprocess that feed the files into clustalo. In Perl you simply issue the linux command wrapped in "\`" either side, in Python use

```
import os
os.system('clustalo -in file.in -out file.out')

```

Ensuing your 'bin' is correctly configured to execute clustalo

or

```
import subprocess # use popen

```

To my knowledge there is no wrapper for `clustalo` unlike several other phylogenetics programs in BioPython, which enormously simplifies pipelines. Others might recommend using shell rather than e.g. Python.

BTW I personally use `Muscle` having been a long term user of clustalo, thats a long discussion.

---

To be honest you should look at the power of NBCI's online blast utility. The download facilities are certainly powerful particularly if you are pulling down e.g. 1500 sequences from blast hits, because it automatically can splice the hits of interest from their genes/genomes. However, its context specific.
> 0 votes

---
Tags: ngs, python, protein-structure, pdb, protein-protein-interaction
---


# Question
Title: Visualizing protein-ligand interactions with LIGPLOT or any other suitable software?
I have used DOCK software to dock 10 ligands into the protein structure.

I want to generate a pdb file that contains all 10 binding modes of the ligands together with the protein. So all together in one pdb file.

Now, what I want to know is when I give this file to LIGPLOT, is it possible to generate something like:

%polar interactions XX %non-polar interactions XX %hydrogen interactions XX

Is there a way to do it with LIGPLOT, if now is there some other software that can be recommended for this type of assessment?

I have been having hard time, trying to find a way or software that can perform similar type of calculation and give me similar type of output that I want.

Best

# Answer
> 1 votes
You can visualize with Discovery Studio Visualizer and Schrodinger. Arpeggio and PLIP can generate output files of the different interactions.

---
Tags: proteins, protein-structure, docking
---


title: Assemble reads for a specific gene

text:
# Question
Title: Assemble reads for a specific gene
I have a lot of unassembled sequencing data and I want to build phylogeny for some genes from these data. I can retrieved the sequene by mapping reads to a reference sequence but if the two species are too divergent, reads will not match and I won't get the full CDS. Is there a better way to deal with that? Maybe I can do that to retrieve the well conserved part and then extend this region with reads that overlap (as if it was gap-filling)? What do you think of this approach? Do you know a software for such extension process? Thanks

# Answer
That's a good question. I remember seeing a poster about making trees directly from reads, but unfortunately, I can not find any other record than this tweet. Perhaps you can try to email David, if he would be keen to share the implementation.

While searching for any information about David's work, I run into a GitHub repository called reads2trees. It's quite old and it's not super popular, so hard to say if it can be useful. But the name matches the problem description quite well.

I also remember there were some attempts (An example paper) to get phylogeny directly from kmers, but I have never tried that on my own. I also think these approached were focused on microbial phylogienies. So again, not sure if it's going to be useful for you.

That's all I got.
> 1 votes

---
Tags: gene, assembly, read-mapping, phylogeny
---


# Question
Title: How to generate the json files for Cromwell workflow execution?
I am starting to use Cromwell WDL scripts for running GATK workflows on a slurm cluster. When using the paired-fastq-to-unmapped-bam script I need to define the fastq sepcifics in the appropriate json file. This needs to be written separately for each pair of fastqs. My question is: **Is there a way of generating the json files on a batch of fastq-files by scripting?** Of course, I could start writing scripts myself, but since I assume that GATK is used in a batch environment, I have the hope that something already exists.

# Answer
Unfortunately generating inputs for Cromwell is not something that is handled by Cromwell itself or by GATK, and generally needs to be scripted to fit your custom needs.

Despite this being a bit arduous, there is a very good reason for this. WDL itself is a highly polymorphic workflow language, where inputs can be vastly different between each workflow. Cromwell *ONLY* understands how to parse the workflow and the inputs, but does not know the context of what it is running (ie this is a paired mapping job), nor does GATK know that it is running within a WDL.

Additionally, I've encountered many different cases where there are not just 2 fastq files but, due to older sequencing technologies, 10 or more for the same sample. Unless GATK/Cromwell knows how your specific samples are arranged there is no easy way for them to generate your inputs.

There are tools that would possibly help you simplify scripting better. WOMTool developed by the Broad has an `inputs` command which will take an arbitrary workflow and generate a `JSON` template of the inputs, which you can then use in your own scripts.

As an example, imagine you have this WDL workflow in a file called `my.wdl`:

```
version 1.0
workflow foo {
  input {
    Array[Pair[File,File]] fastqFiles
  }
}

```

When you run WOMtool's `inputs` command, it produces a template JSON listing all inputs and specifying the variable types expected for each:

```
$ java -jar womtool.jar inputs my.wdl
{
  "foo.fastqFiles": "Array[Pair[File, File]]"
}

```

You can use this template to produce actual JSONs for each of your samples.

Note that you can also get the template written to a file as follows:

```
$ java -jar womtool.jar inputs my.wdl > my.inputs.json

```
> 2 votes

---
Tags: gatk, wdl
---


# Question
Title: How to map gene names to ensembl Ids in R
I am trying to convert a list of gene names into Ensembl IDs but I'm having a bit of trouble.

I've tried using the EnsDb.Hsapiens.v86 package but I always get LRG-ID values returned as well, even when using the GeneIdFilter "startsWith" thing.
My code looks like this:

```
targets <- as.character(targetCSV[, 2])

targetEnsemble <- ensembldb::select(
      EnsDb.Hsapiens.v86,
      filter = GeneIdFilter("ENS", "startsWith"),
      keys = targets,
      keytype = "SYMBOL",
      columns = c("SYMBOL","GENEID"))

```

And for example I get back:

```
    SYMBOL          GENEID
1     CYR61 ENSG00000142871
2      CTGF ENSG00000118523
3    AMOTL2 ENSG00000114019
4    ANKRD1 ENSG00000148677
5    ANKRD1         LRG_379
6    IGFBP3 ENSG00000146674
7        F3 ENSG00000117525

```

I also tried using org.Hs.eg.db package but get a different problem where some of my genes return NA.

My code:

```
  mapIds(x = org.Hs.eg.db,
     keys = targets,
     "ENSEMBL",
     "SYMBOL",
     fuzzy = TRUE,
     multiVals = first)

```

And this works on some but not all the genes for some reason.

```
    $CYR61
[1] NA

$CTGF
[1] NA

$AMOTL2
[1] "ENSG00000114019"

$ANKRD1
[1] "ENSG00000148677"

$IGFBP3
[1] "ENSG00000146674"

$F3
[1] "ENSG00000117525"

```

If anyone can solve either of these problems that would be great, thanks in advance!

# Answer
You can pull a full list of Ensembl gene IDs and MGI symbols from biomaRt in R. From there it should be easy enough to filter for the genes you want:

```
library(biomaRt)
hsapiens_genes <- getBM(attributes = c("ensembl_gene_id",
                                       "hgnc_symbol"),
                        mart = useMart("ensembl", dataset = "hsapiens_gene_ensembl"))

```

You do not find some genes in your example because you were not using the official gene name. For example CTGF is offocially called CCN2 and with this name you will find the appropriate Ensembl gene ID in the biomaRt table. This is why I personally prefer to use gene name that are like `ENSG00000118523|CCN2` in my work. That way you always have a unique identifier (the Ensembl ID) that is unambiguous and you do not run into trouble that for one gene many more or less official gene names exist.
> 3 votes

---
Tags: r, gene, ensembl
---


title: DESeq2 compare within a condition

text:
# Question
Title: DESeq2 compare within a condition
This might be a really stupid question, but I can't figure out how to do this. I've read through the DESeq2 vignette and manual pages but couldn't find an answer.

I have a bunch of samples split up into different conditions (eg. celltypes and disease state). I would like to make a comparison between the two possible values within one condition, but only for the samples with a specific value of the other condition.

For example, with the following samplesheet:

```
patient | phenotype | type
--------+-----------+-----
1       | healthy   | A
1       | healthy   | B
1       | sick      | A
1       | sick      | B
2       | healthy   | A
2       | sick      | A
2       | sick      | B

```

I would like to compare "healthy" vs "sick", but only for type "A".

Currently, I have the following code, but this will also include the type "B" sample:

```
dds <- DESeqDataSetFromMatrix(countData=counts, colData=design, design = ~ patient + phenotype + type)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
res <- results(dds, contrast=c("phenotype", "healthy", "sick"))

```

Any idea how to accomplish this?

# Answer
> 1 votes
Welcome to the world of all possible pairwise comparisons in DESeq2 !

So the easiest way is to create another group in your data frame, first I simulate some sensible counts:

```
counts = counts(makeExampleDESeqDataSet(m=48))
design = expand.grid(id=1:3,
phenotype=rep(c("sick","healthy"),each=2),
type = rep(c("A","B"),2))

```

And you make a group that is a combination of phenotype and type:

```
design$group = paste(design$phenotype,design$type,sep="_")

```

Run DESeq2, including id:

```
dds <- DESeqDataSetFromMatrix(countData=counts, colData=design,
design = ~ id + group)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
res <- results(dds, contrast=c("group", "healthy_A", "sick_A"))

```

So this compares only healthy\_A with sick\_A.. hopefully this is what you need... The other option (normal statistical) is to fit a model with interaction, but I think this is easier to explain and interpret.

Another question you might have is why not do two separate analysis for A and B, I think it's sometimes better if you model the variance with all available data, given the limited number of replicates.

---
Tags: differential-expression, deseq2
---


# Question
Title: After KEGG and GO analysis, how to make tables+phylogenetic trees
hope everyone is ok. I used Trinity to do a de novo transcriptome assembly, then blastp/blastx and then used Blast2GO software to do KEGG and GO analysis. So i got some txt files with header :

for GO

`Sequence name Sequence desc. Sequence length Hit desc. Hit ACC E-Value Similarity Bit-Score Alignment length Positives`

for KEGG

`Seq Pathway Pathway ID #Enzs in Pathway Enzyme`

All well so far. Then i was asked to do 3 tasks.

1. Tables with terpene synthases, (monoterpenes, sesquiterpenes, diterpenes and triterpenes) and corresponding phylogenetic trees including other known terpene synthases from all categories from other plant species.
2. Tables with Cytochrome P450s enzymes implicated in terpene biosynthesis (categorization per clades and clans) and corresponding phylogenetic trees including other known terpene CYPs from all clades and clans from other plant species.
3. From the graphs corresponding to the different KEGG and GO ontology analyses, export the data for the GO analysis for all three parameters (categories, names, blasts, total number of contigs and percentage) to an excel file to be able to transform them in other formats for future publication.

I have no idea how to proceed. Any guidance? Thanks

# Answer
I'm going to answer 1+2 because 3 is hard to understand and maybe trivial?

Each of 1+2 seems to be have these sub-tasks:

1. identify the contigs in your assembly that have the genes of interest.
2. pull out the gene sequences from the contigs (you should have intermediate tblastx output somewhere with coordinates, if you don't you can just re-blast those contigs).
3. possibly translate those in silico to protein sequences.
4. tblastx (gene) or blastp (protein) those gene/protein sequences against whatever db is relevant.
5. pull sequences for the hits above some threshold for each gene from the db.
6. perform a multiple sequence alignment for each gene against the hits.
7. plug the MSA into whatever phylogenetic tool you like. I used phyml a lot back in the day, don't know what the kids like now.
8. make a visualization of the tree, probably?
> 1 votes

---
Tags: phylogenetics, assembly, phylogeny, kegg, go
---


# Question
Title: CCP4 file to a Python 3 numpy array or similar workaround
I would like to merge together several ccp4 formatted density maps (and do a few minor things).

So ideally I would like to open the ccp4/mrc files as numpy arrays in Python 3 and save the array as a ccp4 file. CCTBX can open ccp4 maps as numpy arrays, but I really would like to avoid using Python 2.7 for various reasons.

This is to be used for constraining energy minimisation in PyRosetta using the `pyrosetta.rosetta.protocols.electron_density.SetupForDensityScoringMover` mover/sampler. But its `pyrosetta.rosetta.core.scoring.electron_density.getDensityMap` data loader only likes ccp4/mrc files (density maps) and not mtz files (structure factor maps).

From what I can tell there is no way to fudge a ccp4 map from a mtz map, correct? The fact that only ePDB provides ccp4 maps makes me want to use mtz files.

CCP4's Blend can merge datasets but uses only mtz files...

# Answer
> 4 votes
You can use Gemmi.

```
import gemmi
ccp4_map = gemmi.read_ccp4_map('my.map')
ccp4_map.setup()  # optional
np_array = numpy.array(ccp4_map.grid, copy=False)
# ...
ccp4_map.write_ccp4_map('out.ccp4')

```

The setup() call expands data to the whole unit cell.
The numpy array above accesses the data through Python's *Buffer Protocol* (it does not copy the data). So changing values in numpy array changes gemmi.Ccp4Map.

You can also use gemmi to read mtz file and transform map coefficients to a real-space map, see https://gemmi.readthedocs.io/en/latest/hkl.html.

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