text:
Biohackers Netflix - DNA to binary and video
I'm not sponsored or anything, just interested in their challenge to decipher their DNA code.
They encoded their first episode of "Biohackers" video/binary file to DNA code and said if we could decode it we can watch it (without Netflix).
Here's their page: https://biohackersnetflix.com with description and download for the DNA sequence file.
(Don't know if it's just in German or you can translate it. If questions regarding this page, ask me.)
The file is ~550MB small and contains 3.882.771 lines (not in fasta format). Every line has a length of 147 characters including primers at both ends (Illumina?).
Here are the first 5 lines:
ACACGACGCTCTTCCGATCTCTCCCAGGGACAAAGGTTCTGCATTTGCAGCAAGACTCCTGTAGTGCTGCAGATTCTCTGGTTGGATAGTACGGCGTACATTTCTGTATTGTAGCACCATGGGGTAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTTAAGGCTTCGTAACAGATATTCTATATCGTCACATTGGTCTGAAGGAAGTCGCCTATAATCGCTCCTCTGTTTTTTAAAACTGCTATGGACCCGCTGTTCGGTGGAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTCATGGTATAAGTGTTAAGGGTAATAACCACCTACCCCCCTCATTGCTCGTTTTTCCTGGAACCTTAACATTCGCAATAGCTAGCTGTTTCCTAGTAGAACCAAGGAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTAGGATGTAGTCACAGGTCATTGTCATTAACTCAACCGAGGACATAACACTAAGTCCCACTAGGCCTGGATTCTCTAACGCGGTCTCTCTATTGGGGGAAGGGGTGAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTTCTGGTAAGGCGGGTTGATATCAGTCACCTCCCTTTGAGCTAAAATACGATGGCGATTTAGTGTGAAACTAATAATGCTTGTCATACCAGCAGTACCGGATCGGGAGATCGGAAGAGCACACGTCT

I trimmed all the primers and tried to decode {A, C, G, T} considering every permutation {00, 01, 10, 11} as the obvious(?) decryption method (4! = 24 possible decodings) using python.
Then I hoped to get 1 of these 24 files loaded into VLC media player or something to be played, but it didn't work and every file seemed to be broken in the same way. I think I'm missing something here.
Can I assume that a text file containing only 0's and 1's should be playable in VLC if the DNA code is correctly decrypted?
(If I am wrong here, please tell or move me.)
//Edit: I converted all 24 files to ASCII to see if there's some kind of "video-like header". (All videos have some sort of description in their first lines if opened in text editor?) But there's just gibberish.
//Edit: I saw that every 84th sequence position has a "T", which is kind of weird.
So I tried to run my script again with these T's removed, but still no solution.
//Edit: I searched for "AVI", "264", "codec" and some other strings in every video file I created and hexdumped. Nothing found.
For clarification: I translated the DNA into every 24 binary and then into their ASCII representation following the 19 upvotes answer: https://stackoverflow.com/questions/7290943/write-a-string-of-1s-and-0s-to-a-binary-file. The 104 bases / 208 Bits (removed repetetive "T" and primers) are actually a multiple of 8 (respectively 26 Bytes) so I could be on the right way (even if not 32 Bytes?).
De novo Assembly didn't work and I found no obvious ORF "genes" representing some kind of URL to the video or something which was a neat idea considering the video file would be only ~150MB. (See comments.)

Are the primers all the same and if not are you sure the read direction is not inverted for some? And do any of them align? (I.e. maybe you have to assemble it like a genome).

@MatteoFerla Yes, all primers are the same. And there are not many overlapping regions, so assembly isn't the right way I think.

This fact plus your re-occuring T, suggests it's a problem to be solved with normal encryption cracking techniques, rather than bioinformatics. Say, what is expected at the beginning of an AVI encoded file —you mentioned text-like metadata.

Is your 500 MB file a text file? If so the data present is 2 bits per byte, which would make it a 125 MB binary, which seems low. Had it been a DNA assembly problem, the information would be even less. So I was suspecting it would be an ORF prediction and translation into AA making a URL kind of gimmick.

It's also worth noting that the order in which the sequences are listed in the file might not be the order in which they are to be read once converted to binary.

@MatteoFerla Yes, 500MB text file. That's also the point which surprised me. How would a 40 minute (maybe even HD) Video fit in ~ 150MB? On the other hand their says something like it's the world's first series stored in DNA and "HERE YOU WILL FIND THE COMPLETE DNA CODE OF THE FIRST EPISODE[...]" (then the download link follows). So I think it should be the binary file converted to DNA. But anyways I tried to find the ORF with ORFfinder (min len 100) but nothing suspicious. Also assembly (Velvet and SPAdes) didn't work (not even one contig).

The adaptors 5-ACACGACGCTCTTCCGATCT and 5-AGACGTGTGCTCTTCCGATCT (RC'ed) are parts of Illumina adaptors but are too short. Is there an application that uses such short ones? The T at 104 makes me think there is barcoding going on at the 3' terminus: are bases 105-127 repeated?

I was quite curious when I saw that as well. I've spent more time than I'd care to admit trying the sort of things you have. I've gotten it decoded now, but I can't really claim any sort of victory, as the problem was pre-solved.
After struggling through some of the same experiments you did, I decided to take a closer look at their explainer video here: https://youtu.be/DMYgjOHgHxc
First I tried to decipher the blackboard sketch around 1:30, but it was quite vague. A closer investigation into the speaker, however, yielded some luck. Googling Dr. Reinhard Heckel brought me to his website here: http://www.reinhardheckel.com/ which shows his most recent publication – a paper on encoding digital data in DNA.
The encoding is relatively complex (something nearly impossible to stumble upon accidentally, I'd reckon) but fascinating. The fragments are indexed and have two, layered error-correcting codes. Perhaps making our job harder (but serving the practical purpose of minimizing homopolymers and unwanted annealing) is the fact that the data is XORed with a pseudo-random noise that shuffles the data. The paper obviously has all of the details.
From my skim, however, I stumbled across this Github repo linked in the paper: https://github.com/reinhardh/dna_rs_coding
As it happens, the README was just updated in this last week, describing how to decode the episode from the file provided. If you have Docker, you just need to copy-paste some commands.
The product is indeed the full and final episode – coming in at a dainty 63.1MB. How did they get 40 minutes of video that compressed? Well, in short, 720x360 at 24 fps. It looks a little awful and is in German only, but is certainly a cool little easter egg. If nothing else, I have a cool paper to read.

I haven't thought about going that way, good job! If I get this correctly they inserted error correction codes (ECCs) to the binary video file using some kind of Reed-Solomon-Encoding (RS) approach. So not having this information makes the decoding challenge almost impossible to solve from scratch. Also some specific parameters like size of the ECCs to encode AND decode data have to be known which favors the insolubility. Haven't heard about RS, seems to be kind of complicated but pretty important technique to recover broken/lost information (used in CDs, QR codes and more).

Here in more detail about RS: https://en.wikiversity.org/wiki/Reed%E2%80%93Solomon_codes_for_coders
And a RS implementation in Python:
https://github.com/tomerfiliba/reedsolomon
Works pretty good I played around with it a little.
Very interesting.


text:
Analyzing transcripts for non-coding RNA
From a RNA seq experiment we look for the counts of transcripts of tRNA but we are also interested of other types of RNA we can find. Because of the protocol used for the library prep we can't analyse miRNA but we can analyse ncRNA.
I never analyzed transcripts for ncRNA and I have some questions:

The name of the transcripts are like ENST00000611544.1, ENST00000387069.1, since they are non coding, are they linked to a gene? Should I rename the name of each ncRNA with a gene name or working with these long names is ok?
When analyzing transcripts of RNA I did enrichment analysis to see what pathways are more active. But when analyzing transcripts for ncRNA is possible to do something similar to enrichment analysis? How can analyse the functions of the ncRNA that we get from the study?

Thank you very much in advance

Hi @Mee it might be that a RNAseq expert understands your concerns precisely, however could you provide a bit more background for everyone else please? For example, what package are you performing the analysis with, any description of the samples?

ncRNAs are by itself not special. Sure they might be more unstable and less abundant than other transcripts but for the analysis it is more or less the same. You should take your transcript level counts and aggregate them to the gene level. The tximport package comes to mind. Every transcript (given it is annotated) is linked to a gene. ENST00000611544 is connected to the U2 gene but these snRNAs are also short and I doubt you see them in normal RNA-seq.
In the end you can process data as any RNA-seq. Quantify your data against a transcriptome, e.g. with salmon, then aggregate transcript level abundance estimates (so basically counts) to the gene level and see whether your genes of interest have counts you can work with or whether you mostly have zeros. A transcriptome can be obtaine e.g. from GENCODE. The available annotation files will tell you which genes these transcripts originate from. Check for the GTF files.
Towards enrichment analysis, well it all comes down to the annotations you have. If the ncRNAs are annotated in any pathway then you can perform standard analysis with tools such as goseq or gprofiler2. If they are not annotated, well then there is no information available.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
MM/PB(GB)SA calculations for DESMOND trajectories
I did molecular dynamic simulations with DESMOND.
Now, I want to obtain the interaction free energy using MM/PB(GB)SA.
The free tools I found to perform such calculations are AmberTools and g_mmpbsa.
AmberTools work only on Amber-type files, whereas g_mmbbsa works with GROMACS-type. This imply that I have to convert the DESMOND trajectories to AMBER or GROMACS formats.
I did that using VMD. In case of AMBER, I got the files prmtop and mdcrd whereas for GROMACS I got gro and top files. I also tried to use InterMol, but I am getting an error on conversions.
For both tools, I need to prepare the input files. In case of AmberTools, there is a Python script (ante-MMPBSA) that extract the solvated complex, the complex, the ligand and the receptor from the trajectory files. This scrip is falling to do that from the converted trajectories. In case of g_mmpbsa, it also need the index.ndx file.
My question: Is there a way to accomplish this task? Any suggestion is welcome!

Try here: https://mattermodeling.stackexchange.com/. There's way more VMD questions there, plus there's an entire tag just for GROMACS, and AMBER is on topic there too.

Hi @user1271772 (or I must said @Nick? This is Camps, lol.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
About getting rs id from chromosome and position
Hi I have a dataset with chromosome and position.
VARIANT_ID   chromosome position
17_26797415_147499  17  26797415
17_26797556_147500  17  26797556

How can I assign rs id using R? This data have more than 500 SNPs.

Update: BioMart was used as suggested by @RamRS
'NA NOT FOUND' was a resulting error.
@RamRS recommended inputing na.omit in the input vector
This appear to have resolved the situation

You may want to explore reutils or biomaRt. See this thread on biostars: https://www.biostars.org/p/323319/

Thank you. I tried, but I got a strange error message .
Query ERROR: caught BioMart::Exception::Usage: Filter NA NOT FOUND

That's not strange - one your inputs is NA. Wrap the input vector with na.omit and you should not see the error.

Thank you! I understand it.


text:
PDB residue numbering
The residue numbering in PDB (in my understanding) coincides with UniProtKB align. However, this might not be consistent with that in PDB fasta file. To complicate it further, UniProtKB align can be a few mutations away from the PDB sequence. Disordered/unidentified residues in structures also create gaps in numbering.
Is there a way/suggestion to create a matched residue numbering between sequence and structure? For example
12345678 (residue numbering)
ATPWQMSQ (sequence)
---W-MSQ (PDB residues) "-" for missing residue in structure

Unfortunately, PDB numbering can be one of different schemes. The main three:

Uniprot numbering (canonical isoform)
Construct numbering with or without tag removal
Start from 1 for what has density

However, the are many many corner cases. In Proteopedia there is a lovely article about weird PDB numbering.
The best way to find the details about what part of the sequence is covered is to use the Swissmodel API. It's actually better than either the PDBe and RCSB PDB APIs. Swissmodel is a homology modelling service that maintains close threaded models (>30% ID), but also has PDB data. Here is a Python3 method I use, (provider='pdb') is what is required here).
def get_data(uniprot: str, provider='swissmodel') -> List[dict]:
        """
        Gets data from Swissmodel. Convert each entry with convert_SM_entry
        """
    assert provider in ('swissmodel', 'pdb'), 'Provider has to be swissmodel or pdb'
        url = f'https://swissmodel.expasy.org/repository/uniprot/{uniprot.strip()}.json?provider={provider.strip()}'
        data = requests.get(url).json()
        return data['result']['structures']  # [0] # from to coordinates description

The returned dictionary contains a start and end of the uniprot <--> PDB. It does not however, say where the missing density residues are. But in deposited structures the PDB numbering skips where these residues ought to be... So might not be an important problem.
Another database of interest is SIFTS.
To make a multiple sequence alignment of a non-deposited PDB, you need to extract the sequence and do the alignment here is a snippet using pymol module (from conda):
import pymol2

pdbfile = '' #your pdb file name..
selector = 'chain A' # or whatever chain you want
with pymol2.PyMOL() as pymol:
    pymol.cmd.load(pdbfile)
    fasta = pymol.cmd.get_fastastr(selector)

seq = ''.join(fasta.replace('?','').split('\n')[1:]) # remove header and non AA residues

Biopython has a dozen sequence alignment tools. But the original one still works fine
from Bio import pairwise2
from Bio.Seq import Seq

alignment = pairwise2.align.globalxx(''.join(seq).replace('?',''), "ACQQG")[0]

alignments is a tuple of the MSA sequences.
If you want to do the latter without coding. Extract the sequence from the PDB using PyMOL the normal application and load the pdb and then print cmd.get_fastastr('chain A') in the command line part of the GUI. Copy this. Go to Muscle or other onliner aligner and align!

Thanks @MatteoFerla. The swissmodel API is an interesting one. However, since the rcsb fasta does not have a UniProtKB id, perhaps I have to do an additional alignment to obtain the UniProtKB id? I am surprised rcsb has the sequence data (with the alignment, UnitProtKB id, missing residues, etc), but there is no API/downloadable available.

To go from PDB 4-letter code to Uniprot accession use the SIFTS API or mapping table from their FTP site. Or look manually at the PDB entry.
Do note that a uniprot record has an accession number, say for HRas P01112, and a entry name RASH_HUMAN. The "id" word is an ambiguous and can refer to either, but use to mean the latter. To confuse further, the entry name can change if it is moved from Trembl to Swissprot —say the genome of that organism was sequenced and before the protein was known from an EST read or something.

The RCSB PDB has changed recently, dropping a few annotation from its protein viewer, which I am guessing will be re-added —once that happens I assume a sensible API will appear. How it used to work was that HTML/JS chunks were assembled in the backend, so there was no API for uniprot data and you could not reverse-engineer an ajax call into an API.

Thank you so much. It is indeed a very comprehensive explanation.


text:
workaround for error: Bad CPU type in executable c
I use clustal() from library ape every now and then, but now it doesnt work...
code:
sylvia.clus = clustal(sylvia.seq)

Error:
> sh: /opt/X11/bin/clustalw2: Bad CPU type in executable cannot open
> file
> '/var/folders/46/vt_whvms133190kjmlx5gf4m0000gn/T//Rtmpa82i72/clustald5d11dadcdf.fas':
> No such file or directoryError in file(con, "rb") : cannot open the
> connection

Is there any workaround? (I try reinstalling ape library, but I got the same error)
Screenshot from Rstudio:

Did you recently update the OS on your mac? Try opening running xcode-select --install and try your code again after you run that.

Also you could just open the XCode app on the Mac - it will "install additional components" and then your code might work fine. I do see that there was a recent update.

I get the error in Rstudio (I added screenshot).

I ran:    xcode-select --install   in the terminal and I get 'xcode-select: error: command line tools are already installed, use "Software Update" to install updates'

Did you try opening XCode the app?

re-install ape and clustalw?


text:
Male DNA genotyped with Illumina GSAMD gives two alleles for certain loci within sex chromosomes
I have a tab-separated file (the columns are rsID CHROM POS GENOTYPE) generated by Diganomics with Illumina GSAMD (v2.2) from male DNA, and there are a few hundred loci within X or Y chromosome that contain two alleles, even heterozygous, e.g.
[cut]
rs28463388  X   191998  TT
rs28669107  X   195014  GC
rs6655397   X   201935  GG
rs7890186   X   202458  TG

[skipped]
rs1016793029    Y   8619592 GA

[skipped]
rs527288580 Y   8771750 GA

[skipped]
rs917158084 Y   14871973    TG

[…etc]
There are 1270 such double alleles (7 in Y chromosome and 1163 in the X chromosome) from a total of 734162 number of lines, or 0.17%. What is remarkable is that these 1163 doubles within the X chromosome tend to be clustered, i.e. come consequently one after another mostly, but the 7 in the Y were sparse single lines.

Do you have a link to the manual for the Illumina caller you used? I searched a bit for it but couldn't find anything.

Could it be related to the pseudoautosomal regions of the chromosomes?

@JamesSchinner It seems that you are right. For example, rs28463388 is in the PLCXD1 gene which is a pseudoautosomal gene.

@JamesSchinner, here is a relevant reply: https://www.snpedia.com/index.php/Talk:Gs266

Nice, first time my Biomedical Science degree has been useful!


text:
how to estimate phylogeny without DNA sequences?
I am new to bioinformatics, I am reading ‘Analysis of Phylogenetics Second Edition and Evolution with R’ from Emmanuel Paradis.
I can create phylogeny from DNA sequences, by first calculating the distances and then using (ie nj() function).
My question is: if instead of DNA I have trait data (ie height and weight of species), how can I estimate a phylogeny? (is the same process valid? somehow calculate distances and pass that to nj()?)

You can pass any distance matrix to the NJ algorithm, but if the distances are calculated on height & size you might be building a tree that discerns the entries by age  or health or something other than phylogeny.

Sequence data is basically always better and at this point in history it is de rigueur for any serious study, but originally phylogenies were all based on (ideally discrete) phenotypic characters, so there's no reason why it can't work in principle. As @Pallie mentions in a comment, the clustering methods in e.g. NJ or UPGMA will work on any set of distance measures, and they were in fact originally developed on phenotypic characters rather than sequence data.
I'd recommend taking a look at chapter 24 of Joe Felsenstein's book on phylogenies, which is digitized here (yes the link looks sketchy, sorry). There is a section on p. 425 "Inferring phylogenies and correlations", which discusses how you can in principle use quantitative characters to infer phylogenies. The theory and math for why this is supposed to work is earlier in the chapter.
In many cases people prefer to code the characters as discrete because it makes it easier to work with. But there are also methods that work with the raw quantitative characters.
In R, you should be able to do something like this:
# some simulated data, species are rows and characters are columns
> data
      char1 char2 char3
spec1     0     3     6
spec2     1     4     7
spec3     2     5     8

> library(ape)
> tr = nj(dist(data))
> plot(tr)


text:
Python equivalent to R GRanges and IRanges
Problem:
I am trying to convert some codes written in R to Python and part of that conversion process is find classes equivalent to the GRanges and IRanges from the GenomicRanges R package in Python. https://bioconductor.org/packages/release/bioc/vignettes/GenomicRanges/inst/doc/GenomicRangesIntroduction.html#granges-genomic-ranges
I couldn't find any equivalent library written in Python that does the same operations as the aforementioned R classes. There is however pyranges library in Python but it is not as flexible as GRanges and doesn't have IRanges implemented in it. https://github.com/biocore-ntnu/pyranges
For example in R I could add the following in the GRanges class:
gr = GRanges(
seq = Rle(DF$chr),
ranges = IRanges(DF$Start,DF$End),
ref = DF$ref,
alt = DF$alt,
sampleID=DF$sampleID,
seqlengths = chromSize)

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

@TheUnfunCat, do comment or provide your take, since it's your package :)

I am not sure if your question has anything to do with IRanges. If I get it correct, the limitation you have with pyranges now is adding extra columns to the meta data, like adding values to @values slot in granges.
You can do it using insert or setattr :
import pyranges
import pandas as pd
import numpy as np
gr = pyranges.random(10)

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

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

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

Or directly adding them like this (thanks to @ChrisRands for pointing it out):
gr.alt = np.random.choice(['A','T','G','C'],10)

Thanks @Chris_Rands, I have added it to the answer

thanks alot for your answer @StupidWolf. What about mimicking the IRanges in R? does pyranges support that too?

what function do u need from IRanges? as in.. mimicking IRanges is really broad. Is there something you need to do specifically with IRanges and you want to reproduce that in python?

Mimicking this type of call in R, IRanges(start=1:10, end=11)?

whats the difference between that and creating a pyranges object with pyranges.from_dict({'Chromosome':1,'Start':np.arange(10),'End':11})

you cannot mimick a call... it's really weird... you have a function and it must do something, which is my question in the first place

Yeah, functionality can be "mimic"-ed i.e. with an equivalent implementation, but one cannot expect to Find-Replace function names and expect code to be ported over to a new language. I guess OP has some R code and wishes to implement that in py.

In pyranges you would not need a call like the one you show. You can just convert a dataframe with the columns you mention above to a pyranges like this:
import pyranges as pr
gr = pr.PyRanges(your_df)

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


text:
plink: --update-name vs. editing the BIM
I am working with data from the HGDP public files. These are massive VCFs from whole genome sequencing. To begin with, I converted these to BEDs using plink (1.9). The next problem is that these files use rsid names instead of chrpos names. I would like to change to chrpos. I am wondering whether one can simply edit the BIMs or whether I have to do use plink's --update-name call. The problem with the latter is that it doesn't support duplicates. The HGDP files are full of duplicate variant names because many of the variants don't have rsid's and are just named ".":
$> plink --bfile hgdp_chr21 --update-name updateids --make-bed --out hgdp_chr21_chrpos
PLINK v1.90b6.17 64-bit (28 Apr 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to hgdp_chr21_chrpos.log.
Options in effect:
  --bfile hgdp_chr21
  --make-bed
  --out hgdp_chr21_chrpos
  --update-name updateids

64259 MB RAM detected; reserving 32129 MB for main workspace.
1086522 variants loaded from .bim file.
929 people (0 males, 0 females, 929 ambiguous) loaded from .fam.
Ambiguous sex IDs written to hgdp_chr21_chrpos.nosex .
Error: Duplicate ID '.'.

Thus, one gets an error when trying --update-name and I am wondering if I can simply edit the BIM files to be chrpos instead of rsid. I can't find an answer to this anywhere. It seems that if the BED files don't contain a reference to the variant IDs, but only to the row in the BIM file, then this approach should be possible.
Here's the current BIM file:
user@desktop:/media/luks8tb1/data/genomics/HGDP$ head hgdp_hg38.bim
1   rs1331781659    0   10153   G   A
1   rs1179689498    0   10163   C   T
1   rs201694901 0   10180   C   T
1   rs199706086 0   10250   C   A
1   rs111200574 0   10257   C   A
1   rs145427775 0   10291   T   C
1   rs868413313 0   10297   T   C
1   rs112750067 0   10327   C   T
1   rs1035171912    0   10330   A   C
1   rs868804735 0   10333   T   C

As can be seen, very few lines have the "." variant name. In fact, in the first 100 lines there are 0.

I would reccomend using bcftools on the original vcf files before you convert them to plink, to fill in missing IDs using the command:
bcftools annotate --set-id +'%CHROM\_%POS\_%REF\_%FIRST_ALT' file.vcf

This means you won't have any duplicate ID names and plink should be happy about that. It will preserve the positional information in the .bim file as well.

I like this alternative, but seems inferior to doing it in plink2.0 because here you end up with near-duplicated massive VCFs. In your case, you didn't add compression to avoid this, but one can do this. I ended up using bcftools annotate --set-id +'%CHROM\\_%POS\\_%REF\\_%FIRST_ALT' {f} --output-type z --output hgdp_chr{f_chr}_chrpos.vcf.gz where {f} and {f_chr} is the filename and chromosome in my loop.

yes there is a tradeoff between plink (faster but uses more memory) and bcftools (slower but streams through buffers and has low memory footprint). of course you can pipe together bcftools commands so that you dont have to make as many intermediate files.

I found this 2019 question on biostars in which user zx8754 mentions that plink2 has a command for this purpose --set-all-var-ids, from the plink2 docs:

Whole-exome and whole-genome sequencing results frequently contain
variants which have not been assigned standard IDs. If you don't want
to throw out all of that data, you'll usually want to assign them
chromosome-and-position-based IDs.
--set-missing-var-ids (which just replaces missing IDs) and --set-all-var-ids (which overwrites everything) provide one way to do this. The parameter taken by these flags is a special template string,
with a '@' where the chromosome code should go, and a '#' where the
base-pair position belongs. (Exactly one @ and one # must be present.)
For example, given a .pvar file starting with
#CHROM  POS     ID  REF  ALT chr1    10583   .   G    A chr1    886817  .   T    C chr1    886817  .   C    CATTTT
"--set-missing-var-ids @:#[b37]" would name the first variant
'chr1:10583[b37]', the second variant 'chr1:886817[b37]'... and the
third variant also gets the name 'chr1:886817[b37]'.

Thus, in my case, I used:
user@desktop:/media/luks8tb1/data/genomics/HGDP$ plink2 --bfile hgdp_hg38 --set-all-var-ids @:# --make-bed --out hgdp_hg38_chrpos
PLINK v2.00a3LM AVX2 Intel (27 Jul 2020)       www.cog-genomics.org/plink/2.0/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to hgdp_hg38_chrpos.log.
Options in effect:
  --bfile hgdp_hg38
  --make-bed
  --out hgdp_hg38_chrpos
  --set-all-var-ids @:#

Start time: Tue Aug 25 19:16:10 2020
64259 MiB RAM detected; reserving 32129 MiB for main workspace.
Using up to 8 compute threads.
929 samples (0 females, 0 males, 929 ambiguous; 929 founders) loaded from
hgdp_hg38.fam.
56475980 variants loaded from hgdp_hg38.bim.
Note: No phenotype data present.
Writing hgdp_hg38_chrpos.fam ... done.
Writing hgdp_hg38_chrpos.bim ... done.
Writing hgdp_hg38_chrpos.bed ... done.
End time: Tue Aug 25 19:17:41 2020

Results, BIM:
emil@test-desktop:/media/luks8tb1/data/genomics/HGDP$ head hgdp_hg38_chrpos.bim
1   1:10153 0   10153   G   A
1   1:10163 0   10163   C   T
1   1:10180 0   10180   C   T
1   1:10250 0   10250   C   A
1   1:10257 0   10257   C   A
1   1:10291 0   10291   T   C
1   1:10297 0   10297   T   C
1   1:10327 0   10327   C   T
1   1:10330 0   10330   A   C
1   1:10333 0   10333   T   C


text:
How to know which ortholog database is the best for my work?
I am working on four data sets of yeast genes and one data set of mouse gene. I need the list of human orthologs of those genes. There are a lot of databases available online which is in use for this purpose
https://questfororthologs.org/orthology_databases
I don't know which one to choose and why. Can anyone tell me on what basis should I choose a database over any other database? Thank you.

To recover orthologous sequence for model species
You should be able to find your orthologs using BioMart and your gene IDs.
There are numerous way to do it, directly on the ensembl website or through API. You can start from these links (quite old though) :

http://www.ensembl.info/2009/01/21/how-to-get-all-the-orthologous-genes-between-two-species/
https://www.youtube.com/watch?v=EwfOQsr2O-M #Get orthologs via BioMart
https://www.biostars.org/p/12245/ #Finding Human/Yeast Orthologs Using Embl?

Which database should I use
It is broad and depends on the nature of your work. In your case, working with "major model species", DIOPT (the first line of the table in Quest For Orthology) has the very feature you need on its front page for going from mouse to human. More eukaryotic/viral oriented database will be used for taxonomic identification, annotation inference, etc.
I hope it helps.

Thanks a lot Paul!


text:
How to annotate gene length to a list of gene symbols using UCSC data?
I have a list of HGNC gene symbols, I am looking to get the gene length of each gene. Although I also describe these genes with lots of UCSC datasets as features, so I am wondering if there is a dataset in UCSC I can use to also get gene length from?
I've been looking in the data that is downloadable from UCSC table browser (I've been aiming to find start and ends for each gene to subtract to get gene length) but there are a lot of files and I'm not sure which dataset to take from which will also match to my HGNC gene symbols.

You could do a few command-line operations to answer this question. This assumes the use of hg38 assembly.
First, get a list of genes from GENCODE:
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

Then use grep to filter those genes with your list of HGNC symbols:
$ grep -wFf hgnc_symbols.txt genes.bed > filtered_genes.bed

You can run this through awk to get lengths:
$ awk -vFS="\t" '{ print $3-$2 }' filtered_genes.bed > filtered_gene_lengths.txt

Thank you for this, this looks great. I should've said I am using hg19, my own fault, but I can figure it out from here so thank you :)


text:
Lifting from hg38 to hg37
I am having strange issues lifting this dataset to hg37 from hg38. I am using the HGDP WGS data from ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/. The data set split by chromosome, so I first converted them to BED then merged them to a single file using plink1.9. However, when trying to lift to hg37/19 using either liftOver  or CrossMap v0.5.1 I get strange errors. I tried both the chain from ensemble.org and from UCSC, I tried them both compressed and uncompressed. I am starting to think I need to try the conversion on the VCF files but these are so large I would prefer not to duplicate them unnecessarily.
Data
The data looks fine. If we look at the BIM:
user@desktop:/media/luks8tb1/data/genomics/HGDP$ head hgdp_hg38_chrpos.bim
1   1:10153 0   10153   G   A
1   1:10163 0   10163   C   T
1   1:10180 0   10180   C   T
1   1:10250 0   10250   C   A
1   1:10257 0   10257   C   A
1   1:10291 0   10291   T   C
1   1:10297 0   10297   T   C
1   1:10327 0   10327   C   T
1   1:10330 0   10330   A   C
1   1:10333 0   10333   T   C
user@desktop:/media/luks8tb1/data/genomics/HGDP$ tail hgdp_hg38_chrpos.bim
Y   Y:56887228  0   56887228    C   T
Y   Y:56887285  0   56887285    C   A
Y   Y:56887316  0   56887316    A   G
Y   Y:56887463  0   56887463    A   C
Y   Y:56887491  0   56887491    A   C
Y   Y:56887503  0   56887503    T   G
Y   Y:56887583  0   56887583    G   A
Y   Y:56887631  0   56887631    T   C
Y   Y:56887837  0   56887837    A   G
Y   Y:56887844  0   56887844    C   T

The data are in chrpos variant names here but I had the same result using the original rsid.
Chains
Chains look fine, and I am able to read them into R without problems:
user@desktop:/media/luks8tb1/data/genomics/HGDP$ head GRCh38_to_GRCh37.chain
chain 1 1 248956422 + 10000 297968 1 249250621 + 10000 267719 2
167417 80249 50000
40302

chain 1 1 248956422 + 347968 501617 1 249250621 - 248779253 248932902 3
153649

chain 1 1 248956422 + 585988 12949384 1 249250621 + 521368 13009210 4
1044707 1 0
3709 3 0
user@desktop:/media/luks8tb1/data/genomics/HGDP$ tail GRCh38_to_GRCh37.chain

chain 1 chrY 57227415 + 56821509 57217415 Y 59373566 + 58967656 59363566 11183
85168 0 1
78580 0 1
63496 0 1
3158 1 0
30382 0 1
54147 0 1
80974

CrossMap
user@desktop:/media/luks8tb1/data/genomics/HGDP$ CrossMap.py bed GRCh38_to_GRCh37.chain hgdp_hg38_chrpos.bed hgdp_hg37_chrpos.bed
@ 2020-08-25 20:04:26: Read the chain file:  GRCh38_to_GRCh37.chain
Traceback (most recent call last):
  File "/home/user/.local/bin/CrossMap.py", line 166, in <module>
    crossmap_bed_file(mapTree, in_file, out_file)
  File "/home/user/.local/lib/python3.6/site-packages/cmmodule/mapbed.py", line 32, in crossmap_bed_file
    for line in ireader.reader(inbed):
  File "/home/user/.local/lib/python3.6/site-packages/cmmodule/ireader.py", line 27, in reader
    yield l.decode('utf8').strip().replace("\r", "")
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xdd in position 20: invalid continuation byte

liftOver
user@desktop:/media/luks8tb1/data/genomics/HGDP$ liftOver hgdp_hg38.bed hg38ToHg19.over.chain hgdp_hg37.bed hgdp_hg37_errors.bed
Reading liftover chains
Mapping coordinates
Data format error: expecting at least 3 fields in BED file (hgdp_hg38.bed)

Trying the PED file with liftOver
I converted BED file to PED with plink1.9:
plink --bfile hgdp_hg38 --recode --out hgdp_hg38

Then:
user@desktop:/media/luks8tb1/data/genomics/HGDP$ liftOver hgdp_hg38.ped hg38ToHg19.over.chain hgdp_hg37.ped hgdp_hg37_errors.ped
Reading liftover chains
Mapping coordinates
invalid unsigned integer: "HGDP00001"

There's a few question on Google about this, but nothing that looks useful.
CrossMap does not support PED.

BED format can refer to the ".txt" file format that stores chromosomes and positions of genomic features as well as some other things like name. bed format can also refer to the plink format, which belongs to bed, bim, fam trio. The plink bed format stores all the genetic variants, but the program needs to refer to bim file for their positions.
Genome liftover obviously only needs the coordination of each SNPs, and it does not need individual genotypes at all. In fact, the bed format in PLINK does not provide any useful information. Knowing this means that providing the plink bed file is wrong. And this is pretty clear from the error message.
Moreover, since the three plink file always needs to be used together, most programs refer to them as plink file format. So whenever people say BED format, it is the txt file format.

Not sure about crossMap but liftOver requires the BED format where column 1 is the chromosome, column 2 the 0-based start and column 3 the end coordinate. Additional columns are optional and can be used to store the other columns you have. Simply convert BIM to BED and run liftOver.

BIM files are part of the plink trio of BED, BIM and FAM. The data above are already in BED format.

Data format error: expecting at least 3 fields in BED file (hgdp_hg38.bed)
According to this, they are not.


text:
Populations genetics and dynamics of bacteria on a Graph
Disclaimer:
I'm a physicist and I'm fairly new to Bioinformatics therefore somethings below may not make sense.

My purpose:
I would like to simulate genetic evolution of bacterial populations, implementing interactions among different bacterial strains on a graph.
My starting point is to generate a graph with a certain topology.
Each vertex $V_i$ should contain a population of $N_i$ individuals for which I choose a certain allelic profile from an MLST database. Population numerosity and allelic profile is assigned by random choice.
$$\\\\$$
At starting time each vertex is populated with a "pure" strain of $N_i$ clones, identified with a randomly choosen allelic profile.
Dynamics assumptions
I will assume that:

Topology does not change with time

Vertex dynamic flows as a Moran process:
At time $t$ , one individual is chosen randomly to reproduce and one
individual is chosen to die. The same individual can be chosen to reproduce and
then die. Thus, an individual has either zero, one, or two descendants. Zero and
two with equal probability $p_0 = p_2 = (N_i -1 )/N_i^2$ , and one with probability $p_1 = 1 − 2 p_2$.

Migration from one Vertex to one another is possible for nearest neighborhood with some fixed probability depending on choosen topology. I will assume that when migration occurs $I$ new individuals enters into $V_i$ and $O$ individuals exits from $V_i$, with $I$ and $O$ depending on origin vertex numerosity. Migrants carry allelic frequencies according to population they came from. When migration occurs a new population is established: allelic profile frequencies changes according to that.

Migration and Moran process evolution does not happen simultaneously; This allows me to guarantee that each population performs at least one evolution.

Questions:
(1) Is there something fundamental that I'm missing?
(2) Moran model ensures that the population size remains constant (apart from migration effects obviously). Is this appropriate in order to simulate dynamic in stationary growing phase or should I introduce some demographic effects?
(3) Which could be a proper "timing" for migration and evolution steps?

Just a thought... A key concept in population genetics is fitness of an individual/clonal group under that environment. A change in environment results in a change in the fitness and as a consequence you have a bottleneck, where individuals die when their fitness is lower than a threshold. Specialists have high fitness in one environment, but low in others, while generalists are mediocre at all. So you could choose who dies based on some fitness-threshold effect where fitness is inversely tied to say doubling time, if that is a parameter that is allowed to evolve.

I forgot to mention that my naive implementation of vertex dynamics follow a neutral model because all individuals could be picked as a parent with uniform probability distribution (this should be like an equal fitness).
So do you suggest that I should implement fitness variation as an effect of migration? (e.g. incoming bacteria have lower fitness due to new environment)

I assume you are simulating a null distribution. Are you investigation recombination??
My main advice is to use population genetic terminology rather than geomometry to describe your simulation (e.g. vertex is a migration event between allopatric populations). Migration is investigated between populations investigated via $F_{ST}$. You appear to be assuming clonality, and population movement appear to follow a 'stepping stone model' (I might have misunderstand the question).
Mutation rate ... I'll have a guess at $10^5$ per site per year. Migration rate depends on the bacterial species being modeled, eg. MRSA spread in a very similiar fashion to that described in your model (along a moterway), in contrast soil bacteria don't adhere to allopatry at all (when I looked at them).
Fixing one parameter (population size) whilst investigation a second (migration) would be normal to assess the parameter space initially. Bacterial populations do undergo selective sweeps and using an appropriate growth model I would look at a standard population growth model (e.g. the one Bob May investigated in chaotic interactions within a differential equation). There are several population growth models and fixing migration whilst assessing the model would seem sensible. The parameter space in the combined model of migration, space and population size becomes complex.

Okay I get your model now. Its fine, the thing I'd be cautious about is the population density of the bacterial host because this would skew migration across the nodes of your graph. A quick literature review will demonstrate the importance,  double paeroto log normal distribution and its cousin (the name escapes me .... ) are sort of key words that might flag up. Bacteria are not my active thing so I'm rusty.

Let me explain better:
my graph G(V,E) has nodes (aka. vertex) and edges (aka. links). My dinamics is parametrized by two discrete time-variables, $t_1$ and $t_2$ . $\Delta t_1$ is the moran process typical time and $\Delta t_2$ marks each migration event.
Clearly there's a relation beetwen $\Delta t_1$ and $\Delta t_2$ that I should implement.
Allopatric hypothesis is plausible only at $t=t_0$ when first migration has not occurred yet.
Yes, population movements follows a 'stepping stone model'.
Recombination is not taken into account yet (yes it's my basis null model)

Okay, i see. Your model seems fine. You're looking at bacterial movement between foci under investigation, e.g. hospitals that wort of thing, this is represented by nodes on your graph. I presume you have observed data and are generating a null distribution. The one thing I'd watch ... is there seasonality in the population dynamics of the species? The other thing important is population density of the bacterial host ... of the two I would say the latter is most important.


text:
Sort reads in BAM file based on presence of specific deletion?
I have an indexed BAM file containing long-read sequencing data and I'd like to split the reads contained within into those with a known deletion and those without the deletion (I have the deletion coordinates available to me) when mapped against the hg38 genome. Ideally I'd like to end up with two BAM files (one with reads containing the deletion and one without) for downstream analysis. I've searched a bit online and haven't really found any tools that do exactly that, but this seems like a common enough task that there has to be something.
Thanks in advance for any help/guidance!
UPDATE
In addition to @gringer's great answer, I decided to write a small python script myself that solves my issue using pysam, the fruits of which can be found here. If anyone has any feedback or suggestions on it I'll happily take them, as I'm still very new to this. It's undergone very little testing, so user be warned.

I was asked to write some code sometime near the end of last year to do this, so that viral sequences could be split based on variants at a particular location. The code I wrote doesn't actually do the splitting - it creates SAM read groups based on those variants - but that splitting can be done after processing with this script using samtools split.
As with most of my scripts, it hasn't had much testing / use:
https://gitlab.com/gringer/bioinfscripts/-/blob/master/samVarSplitter.pl
usage:
samtools view -h mapped_reads.bam | ./samVarSplitter.pl [-ref <ref>] [-pos <int>] [options]

example:
samtools view -h mapped_reads.bam | ./samVarSplitter.pl -ref Wuhan-Hu-1 -pos 11083

If the input file contains a lot of reads that don't include the variant location, I'd recommend pre-filtering reads to only include those that cover that location, because expanding the CIGAR string and storing a record of the positions can take some time.
Note: currently the position argument only supports a single location. I could extend this to multiple positions (i.e. a position range) if that were a desired feature.

Thanks so much, this is pretty much what I was looking for! If I understand it correctly, the script only works with a single position though correct? The deletion I'm looking at is 9bp in size, but I suppose I can just use one of the positions as a stand-in for the entire deletion.

The CIGAR string should tell you if a read has a deletion in it.  So filter the bam files based on that.

The OP is looking for a known deletion, so if you look to identify reads in the region with 'D' in the cigar, you might include other deletions


text:
BAM file filteing to remain best isoform
I ran HiSat2, MarkDuplicate, removed reads with the lower quality score than 40 and finally only kept properly paired reads.
After the BAM filtering steps, I used the Scallop results with TransDecoder. Next, I used bamCoverage to create RNA-Seq profile. I noticed that not always the first isoform is the best one compared to the RNA-Seq profile as could be seen below a few examples:
1.

2.

3.

4.

5.

By any chance, is there a way to keep only the best isoform which is in red box?
Thank you in advance,

I'm not familiar with Scallop/TransDecoder, but it would be helpful to know what you mean by "best". I think I agree with you about which one is "best" but formulating a specific quantitative definition of "best" would be a good first step. Is there any definition of "best" in these tools? How would they evaluate what is "best"?

Unfortunately, the additional information was too long to post them on this side, Therefore, I posted them here


text:
tabix indexing a GTEx variant file- ':' delimited
I have a gtex variant file- the head of which looks as follows:
phenotype_id                                    variant_id
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_13550_G_A_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14671_G_C_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14677_G_A_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_16841_G_T_b38  ...

Effectively I would like to be able to lookup the variants inside a particular gene (the file is sorted) and put them in a temporary file:
E.g. If the variants in gene "ENSG00000148481- MINDY3" are on line number 87528225 to 87536766- effectively what I would like is the equivalent to zcat file.gz | sed -n '87528225,87536766p' > MINDY3.txt. However zgrep ENSG00000148481 file.gz is just as fast as the above...
Hence I thought tabix would be the right tool for this-
I would like to tabix index it to make lookups faster. It is gzip compressed and I will firstly have to do:
zcat gtex.txt.gz | bgzip > gtex.txt.bgz

However I am not quite sure how to proceed from there given that the data is not tab-delimited.
As a trial I tried the first 1000 lines:
zcat gtex.txt.gz | head -n 1000 | bgzip > gtex_1000.txt.bgz

./tabix -p bed gtex_1000.gz #index as a bed file
[get_intv] the following line cannot be parsed and skipped: chr1:15947:16607:clu_36198:ENSG00000227232.5        chr1_13550_G_A_b38     ......

./tabix -p vcf gtex_1000.gz #index as a vcf file
Indexing as a bed file results in a warning while indexing as a vcf gives no warning yet either way when I try to retrieve a sequence:

./tabix test.gz chr1:15000:17000

It returns nothing.
I am starting to think that I will just have to write a script that splits on the ':' and writes the data to a new file.... and then index that file- which will take a huge amount of time
... Does anyone know of a trick to index the files with unconventional delimiting?

Suppose that you could successfully creat an index file for your txt file, which program would be used to read it?

@PhoenixMu see the way that I have used sed in the updated question... this was my plan... although this was the same speed as grep

I think tabix only works for bed and vcf format. Your file is not in bed format?

Yes it's not in bed format as the coordinate columns are : delimited and the file is too big to justify producing a copy with properly delimited columns- I am in the process of making an index file myself and will post an answer if it works

I proceeded as follows:
import gzip
import sys

qtls=['sqtls','eqtls']
tissues=['Adipose_Subcutaneous', 'Adipose_Visceral_Omentum', 'Adrenal_Gland', 'Artery_Aorta', 'Artery_Coronary', 'Artery_Tibial', 'Brain_Amygdala', 'Brain_Anterior_cingulate_cortex_BA24', 'Brain_Caudate_basal_ganglia', 'Brain_Cerebellar_Hemisphere', 'Brain_Cerebellum', 'Brain_Cortex', 'Brain_Frontal_Cortex_BA9', 'Brain_Hippocampus', 'Brain_Hypothalamus', 'Brain_Nucleus_accumbens_basal_ganglia', 'Brain_Putamen_basal_ganglia', 'Brain_Spinal_cord_cervical_c-1', 'Brain_Substantia_nigra', 'Breast_Mammary_Tissue', 'Cells_Cultured_fibroblasts', 'Cells_EBV-transformed_lymphocytes', 'Colon_Sigmoid', 'Colon_Transverse', 'Esophagus_Gastroesophageal_Junction', 'Esophagus_Mucosa', 'Esophagus_Muscularis', 'Heart_Atrial_Appendage', 'Heart_Left_Ventricle', 'Kidney_Cortex', 'Liver', 'Lung', 'Minor_Salivary_Gland', 'Muscle_Skeletal', 'Nerve_Tibial', 'Ovary', 'Pancreas', 'Pituitary', 'Prostate', 'Skin_Not_Sun_Exposed_Suprapubic', 'Skin_Sun_Exposed_Lower_leg', 'Small_Intestine_Terminal_Ileum', 'Spleen', 'Stomach', 'Testis', 'Thyroid', 'Uterus', 'Vagina', 'Whole_Blood']

gene_old=''

def index(tissue,qtl,gene_old):
        counter_old,counter=1,1
        if qtl=='sqtls':
                path='/path/to/file/'+tissue+'.v8.sqtl_allpairs.txt.gz'
        else:
                path='/path/to/file/'+tissue+'.allpairs.txt.gz'
        f = gzip.open(path, 'r')
        next(f)
        for line in f:
                try:
                        line=line.decode()
                        line=line[:line.find('\t')]
                        gene=line[line.rfind('ENS'):].split('.')[0]
                        #print(gene,gene_old,counter_old,counter)
                        with open(qtl+'_'+tissue+'.idx','a') as indexed:
                                 indexed.write(str(gene_old)+'\t'+str(counter_old)+'\t'+str(counter)+'\n')
                        gene_old=gene
                        counter +=1
                        counter_old=counter
                except IndexError:
                        print('tissue,qtl,gene_old,counter')
                        counter+=1
        f.close()

for qtl in qtls:
        for tissue in tissues:
               qtl=row.qtls
               tissue=row.tissues
               print(qtl, tissue)
               index(tissue,qtl,gene_old)

The above script writes an index for each gene which should look as follows (took 2 days to finish all files):
ensid    row_start    row_end
......
ENSG00000230337     1526391 1534329
ENSG00000171819 1534330 1542321
ENSG00000198793    1542322 1550178
ENSG00000120942  1550179 1558010
.....

And to get a gene of interest you will have to look up the coordinates for an ensemblid/gene and for example:
zcat Brain_Nucleus_accumbens_basal_ganglia.v8.sqtl_allpairs.txt.gz | sed -n '1550179,1558010p;1558011q' > ENSG00000120942.txt

Tabix files need a name and a genomic coordinate at least. Tabix can be (ab)used by creating a fake genomic coordinate and then it just indexes names but not sure I would recommend this usage. What is your actual expected application of this?

Hi @Colin thanks for getting to me- I just want to be able to take out all variants inside a particular gene- for example if the variants in gene "ENSG00000148481- MINDY3" are on line number 87528225 to 87536766- effectively what I would like is the equivalent to zcat file.gz | sed -n '87528225,87536766p' > MINDY3.txt
Hence I thought tabix would be the right one for this

Here is an example of a tabix file that has been hacked up to allow "free text query". http://target.wustl.edu/ENCFF000ARO_ALL.iteres.loci.gz  So they can query tabix file.txt "A(n)" and because there are many lines starting with A(n) it will return all of them

Can also check out other options for text indexing https://unix.stackexchange.com/questions/3086/command-line-friendly-full-text-indexing


text:
Error When Using biocLite as an installer in rpy2 python library
I am running into an error when I tried to run the following lines of code, the error at line 3:
1 from rpy2.robjects.packages import importr
2 base = importr('base')
3 base.source("http://www.bioconductor.org/biocLite.R")
4 biocinstaller = importr("BiocInstaller")
5 biocinstaller.biocLite("seqLogo")

rpy2.rinterface_lib.embedded.RRuntimeError: Error: With R version 3.5 or greater, install Bioconductor packages using BiocManager; see https://bioconductor.org/install
How can I use BiocManager installer instead of biocLite  and reflect that on my code? I couldn't find any documentation or source that talk about this error from within Python.
Thank you in advance

you need upgrade your R base version.

Definitely recommend installation via Conda as @zorbax points out. This normal stuff for installing rpy2 As a word of warning further errors with rpy2 are also normal and as a general rule of thumb rpy2 is a nightmare. For R's ggplot2 via rpy2 is outdated and will not import the full funcationality of the ggplot2 into Python3. This is irregular given have important ggplot2 is as a graphing package. After my experience I would use subprocessing into R rather than deal with rpy2 again

You're using an outdated version of R. You can use the R 4.0 from conda-forge channel, just create an environment with r-base and then install rpy2, the latest version compatible with R 4.0 is available through pip
conda create -n rpy2 -c conda-forge r-base
conda activate rpy2
pip install rpy2

Now you can use rpy2 with R 4.0 in your python code.
import rpy2.robjects as robjects
robjects.r('install.packages("BiocManager", '
           'repos="http://cran.r-project.org")')
robjects.r('BiocManager::install("seqLogo")')
robjects.r('library(seqLogo)')

Remember deactivate your env with conda deactivate.

Thank you so much this work! but when I do the same thing for "GenomicAlignments" package I get error:
raise embedded.RRuntimeError(_rinterface._geterrmessage())
rpy2.rinterface_lib.embedded.RRuntimeError: Error in (new("nonstandardGenericFunction", .Data = function (expr, envir = parent.frame(),  :
object 'GenomicAlignments' not found

I can't replicate you error. I did robjects.r('BiocManager::install("GenomicAlignments")') and then robjects.r('library(GenomicAlignments)') and it works.

I am unable to do this for pcalg (https://cran.r-project.org/web/packages/pcalg/pcalg.pdf). Can I see a minimal working example?


text:
RCSB API: limited returned result
Is there a limit of query results on RCSB API?
Interestingly, no matter how to adjust my query criteria, only 10 assemblies are returned. Here is an example query json: (truncated from a longer query file)
{
  "query": {
    "type": "group",
    "logical_operator": "and",
    "nodes": [
      {
        "type": "terminal",
        "service": "text",
        "parameters": {
          "operator": "greater",
          "value": 1,
          "attribute": "rcsb_assembly_info.polymer_entity_instance_count"
        }
      }
    ]
  },
  "return_type": "assembly"
}

Example python code requesting query:
url_query = 'https://search.rcsb.org/rcsbsearch/v1/query?json='
query_str = open(query_json, 'r').read()
query_str = query_str.replace('\n', '').replace(' ', '')
url = url_query + query_str

r = requests.get(url)
content = r.content.decode('utf-8')
result = json.loads(content)
print(len(result['result_set']) # always equal to 10

Maybe ask RCSB at https://www.rcsb.org/pages/contactus.

Staff from rcsb kindly reminds me the solution - add the following lines in the json query.
"request_options": {
"return_all_hits": true
}


text:
Microarray Meta-Analysis and Cross-Platform
I have a microarray dataset with Agilent platform, all of the samples are cancerous, I need normal or control samples to compare with cancerous samples, can I combine control samples from another platform (for example Illumina) to Agilent platform? and then remove the batch effect (by ComBat function in the sva package).
I would be grateful if you could guide me to do it.

Hi there isn't sufficient detail really but I'll pop an answer below

In short the answer is no because you have to standardise between the arrays. Thus both arrays would need to run exactly the same control sample and you would then perform statistical standardisation (its a formal term), which is a varience measure where you look at the variation around the 'mean' for all samples in the array and here the 'mean' is the standard sample. If that can be done for both arrays at that point they can be combined.
The problem is intensity per run varies and you have to control against this

In this scenario you  can certainly combine samples by assessing whether the control samples between runs approximates to identity. What you could use is a uniform distribution to assess this, variance must MUCH smaller than the mean for the STANDARD sample(s). I don't really do uniform distribution statistics, so I'd need to research this, because 'uniform' means identity to standard deviations against the uniform distribution is a bit meaningless. However again the uniform distribution is denoted by a very variance in comparison to the mean. The statistic of dividing mean by standard deviation/variance from memory is called the Index of Association (I could be wrong). If this value is at an extreme you're good to go <<<1, if it is remotely close to 1 you're in trouble, if it is >1 you need to use standardisation.
Presumably you have several standards to ensure a mean is meaningful. A T-test and non-parametric equivalents is not stringent enough. If there is a difference assess a correction factor and rerun the uniform distribution check. The same correction will be applied to all samples.
The undisputed method to combine is standardisation, but accept in this scenario the differences between runs is potentially very small, so here what you are attempting is essentially to combine raw data or at least keep one of the runs as pure raw data.

Apologies for the delay in replying. Thank you for your valuable guide and useful information.
Can I combine control samples from the same platform (Agilent) to the desired dataset?

@MortezaHadizadeh I've answered above. If this answers your question please tick as 'accepted'.


text:
Dealing with duplicate keys
I got the following error:

attempting to use to_dict to index a multiple sequence FASTA file with duplicate keys.

Are there methods for working around duplicate keys in biopython? I was thinking to simply remove the duplicate entries, but I would need to delete the entire entry associated with the keys.

Hi @Mike are you looking to rename the duplicate keys?

yes, I was thinking of appending the keys with digits incrementally or something? The keys are uniprot IDs, but the renamed duplicates wouldn't do anything harmful to my output. Might just be best to remove them using a set method or something. Still quite new to this.

Removing keys via a dictionary is easy, it should be automatic (they will overwrite). Pandas will enable duplicate keys are preserved by shifting the index

Yes, Biopython deliberately does not allow duplicate keys with .to_dict(), e.g.
Input file:
$ cat test.fa
>1
ATG
>2
AAA
>1
CG

Code:
from Bio import SeqIO

print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta")))

Traceback:
  File "script.py", line 3, in <module>
    print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta")))
  File "/path/lib/python3.7/site-packages/Bio/SeqIO/__init__.py", line 789, in to_dict
    raise ValueError("Duplicate key '%s'" % key)
ValueError: Duplicate key '1'

However, you can easily write your own function, e.g.:
from Bio import SeqIO

def to_dict_remove_dups(sequences):
    return {record.id: record for record in sequences}

print(to_dict_remove_dups(SeqIO.parse("test.fa", "fasta")))

Output:
{'1': SeqRecord(seq=Seq('CG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[]), '2': SeqRecord(seq=Seq('AAA', SingleLetterAlphabet()), id='2', name='2', description='2', dbxrefs=[])}

Alternatively, you Biopython's .to_dict() accepts a key_function argument, e.g. if you wish to make all the headers unique
from Bio import SeqIO
from itertools import count

c = count()

print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta"),
                    key_function = lambda x: f"{x.id}_seqNum{next(c)}"))

Output:
{'1_seqNum0': SeqRecord(seq=Seq('ATG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[]), '2_seqNum1': SeqRecord(seq=Seq('AAA', SingleLetterAlphabet()), id='2', name='2', description='2', dbxrefs=[]), '1_seqNum2': SeqRecord(seq=Seq('CG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[])}

thank you; return {record.id: record for record in sequences} so the expression garners us only unique IDs it seems, is this the case?

Well, obviously it does, but I guess my question is, what is determining the order (in the first solution)?

@Mike The last key,value pair is nearly* always retained *https://stackoverflow.com/questions/39678672/is-a-python-dict-comprehension-always-last-wins-if-there-are-duplicate-keys

Hi Chris, got it thank you!


text:
COVID-19 GWAS: what is a cross-replicating association?
In a genomic study of patients infected with SARS-CoV-2, the authors  detected cross-replicating associations with rs11385942 at locus 3p21.31 , the association signal spanned the genes SLC6A20, LZTFL1, CCR9, FYCO1, CXCR6 and XCR1
I would like to know what's a cross-replication association.
Source pdf is here

Replication in GWAS is an analytical term, approximating to 'robustness'. It is not biological phenomenon of a genome replicating. What the authors are saying is these genes associate with an allelic variant(s) which correlate with the disease phenotypic of COVID-19.
Three of those genes are very interesting indeed notably CCR9, CXCR6 and XCR1, because these are chemokine receptors. The later for example is associated with regulating between the innate and adaptive (CD8+, i.e. cytotoxic T-cell) immune responses (described in the image below). The innate/adaptive interface is an extremely trendy area in immunology.
The authors are implying the immunological signalling associated with one or more of these chemokines is fundamental to combating COVID-19 infection. 'Replication' is simply the GWA geek for 'its real'.


text:
FPKM, FPKM-UQ, TPM or counts: How do I know which kind of unit should I use?
I'm trying to download data from the TCGA for gene expression analyses in R, but I'm in doubt if I should use FPKM, FPKM-UQ or counts? When the dataset is in counts, I suppose it's raw data, isn't it? So what's the best unit to compare multiple datasets?
I'm planning to use limma or Dseq2 for GE analysis and found that with Dseq2 I need to use count(non-normalized???) data... is that correct? so what's the best package and working strategy?
thank you very much, Fabiano

Type of data you need depends on the downstream applications and since you would like to carry out DEA with DESeq2, you would need raw counts (non-normalized).
There are many ways to import/download TCGA data, one such tool, the TCGAbiolinks package gives a nice interface for not only downloading the read count (or pre-processed) data but also associated clinical data.

yes, I'm using TCGAbiolinks to download the data... what's the difference between limma and DESeq2? many thanks, @haci

I'm planning to compare my gene of interest in different types of cancer with normal tissue.... How many datasets of each kind of cancer and normal tissue should I include?

limma was designed for microarray data, but can be used for RNA-seq by using the voom transformation. However, the gold standard these days is DESeq2. As for how many datasets you need... what's your hypothesis? What are you looking for? Do you just want to blindly do tumor vs normal for multiple cancer types and check what comes up?

Hi @csgroen. yes, I just want to investigate my gene of interest in different types of cancer. Maybe I don't even need an extensive differential expression analysis. can I just compare my gene of interest in cancer cells with normal tissue using t student? it is probably enough to answer my question..

If it is just for gene of interest, you can maybe do it using XenaBrowser... easier than downloading all the data using TCGAbiolinks then doing differential analysis.

@csgroen in XenaBrowser the number of datasets is small, in the TCGA I can access dozens or even hundreds of datasets of the same cancer type...

I'm not sure what you mean. TCGA was a consortium that produced around 30 cohorts of various cancer types. But for each cancer type, there is only one TCGA cohort. All of them are available in XenaBrowser. Now, in GEO, you'll have hundreds of cohorts of other sources, but they are not 'in the TCGA'

@csgroen... no, look at the GDC website (https://portal.gdc.cancer.gov). The TCGA counts with 33 cancer types, but 585.428 datasets are available.. 86.582 datasets are RNA-Seqs..

FYI, it's 86k samples, not datasets. And all the data from the GDC is available in other platforms, like XenaBrowser. But you do as you think is best ;)

You should use "counts" from the list of gene expression metrics that you provided rather than either FPKM or FPKM-UQ. FPKM is a normalized gene expression metric and an acronym for fragments per kilobase per million mapped reads. FPKM normalizes for both sequencing depth and genome size. Those specific packages want non-normalized counts as the expression metric. This means the number of reads mapped to a particular gene or feature. This is a slightly ambiguous metric but it can be the number of single-end reads mapped to a loci, the amount of pair end reads counted towards a loci, or the amount of super-reads (the theoretical sequence that spans the pair of paired-end reads which includes the presumed sequence between the first and second strand of a set of pair end reads) counted toward a loci. These values should be integers based on the assumptions of RNA-Seq and microarray technology which measures RNA transcripts present in a sample and precludes transcripts from spanning multiple genes. It does not matter which of these metrics you use as long as you are consistent. Those Bioconductor packages provide several different normalization methods but rely on read counts as their input.


text:
How are different forcefields developed/evaluated?
The Rosetta software can use a variety of different energy functions:

Ref2015
GenPot
Beta-Nov16

How are these energy functions developed and validated? Is it common to do a single protein design with several of these energy functions?

To misquote Bismarck, "If you like forcefields and sausages, you should never watch either one being made".

Some force-fields rely solely on theoretical physical principles, while others are tweaked specifically for the goal of accurately depicting certain structures like sheets and helices or movement/diffusion aspects at different theoretical temperatures of a set of well-studied proteins. (Top-Down vs. Bottom-Up)
I haven't done MD myself but I am in an MD group and people seem to normally stick to a single force-field for a given MD project but it happens that they do very short simulation runs with force-fields just as a test.
They definitely don't plan on doing the full simulation on different force-fields unless the comparison of those is their main goal.

The Rosetta ref2015 is well described in two papers, which I wish I had read much earlier, see this JChIM paper.
It's a statistical forcefield, where the terms are weighted to best match the empirical total energy in kcal/mol and the paper does a swell job at explaining what they do.
For several Rosetta experiments, you switch between full-atom and centroid (=coarse grain). These have different weights, as you can see from the files in the database folder or via pyrosetta, but majorly they have different "params" (topologies) for the residues as the centroid lacks many atoms as it's a simplification. In some instances, you may switch between dihedral (normal) to cartesian (explained in the FastRelax paper), wherein the former used internal coordinates (e.g. atom A is 1.5 Å away from B and A,B,C is a 107ª angle etc —cf. a params file ICOOR block), while the later is cartesian coordinates ( A is at [x, y, z] —cf. PDB file).
The GenPot uses its own AtomTypes, which means it uses its own params for each residue —like centroid/fa. The topology of a residue (ResidueType) is based on its conformation (ICOOR), bonding and atom definition, which is atom name, gasteiger-marsili charge and AtomType (an amalgamation of element, sp-hybridization, size, H-bond donor/acceptor etc.) These AtomTypes are better for ligands as many groups behave different from peptide atoms —say, a nitrile nitrogen is definitely different from a deprotonated aromatic nitrogen.
Some application, such as hydration, ligand docking etc. have specifically tweaked scorefunctions for use by their movers. FastRelax is a nice example as the cycle weighting changes are spat out by the tracer at vertiginous speed, where the protocol is:
repeat 5
coord_cst_weight 1.0
scale:fa_rep 0.040
repack
scale:fa_rep 0.051
min 0.01
coord_cst_weight 0.5
scale:fa_rep 0.265
repack
scale:fa_rep 0.280
min 0.01
coord_cst_weight 0.0
scale:fa_rep 0.559
repack
scale:fa_rep 0.581
min 0.01
coord_cst_weight 0.0
scale:fa_rep 1
repack
min 0.00001
accept_to_best
endrepeat

As you can see, the repulsion term of the Lenard-Jones is going from 0.04 up. These values were chosen not for thermodynamic reasons, but worked best for the application (energy minimisation of the structure).
However, this switching is all internal. Scores are always reported with the same scorefunction and comparing between different ones is not advised even for the minority of scorefunctions calibrated to be in kcal/mol.


text:
How to measure the total size of a fastq file in base pairs?
Or Kbps/Gbps. It feels like it should be conceptually very simple, but I can't seem to figure out the right combination of keywords to find it via my search engine. Help would be appreciated!
I have BBMAP, SRAtoolkit and MEGAHIT already installed, and also use bash. I'd be very happy if this can be answered with software that I already have, but if not that's perfectly fine.

I've been using this:
cat file.fastq | paste - - - - | cut -f 2 | tr -d '\n' | wc -c
Explanation :

paste - - - - : print four consecutive lines in one row (tab delimited), to merge the info for each read

cut -f2 : print only the second column, to access the sequence after the paste

wc -c : count the characters

tr -d '\n': to remove from count the eventual newline characters

(a tip for your googling: try search for "counting number of bases in fastq file")

The number of bases in a fastq file can be counted in bash with awk and wc
awk 'NR % 4 == 0' ORS="" fastqfile|wc -m

The awk code prints every fourth line (which is actually the quality scores, but that doesn't matter here). wc -m returns the total number of characters.
Edited to add ORS="" to prevent count of newlines (the double quotes are not strictly required)

Using Perl:
Print the sequence lines (line number 2, 6, 10, etc). Remove the newlines with chomp. Count the bytes (here, bases) using wc -c:
perl -ne 'if ( $. % 4 == 2 ) { chomp; print; }' | wc -c

Using a fastq parser in Python (for instance the one provided by mappy), it is rather straightforward to compute the sum of the lengths of the sequences present in a fastq file:
$ python3 -m pip install mappy  # Only if mappy is not installed yet
$ python3 -c 'from mappy import fastx_read; print(sum(len(seq) for (_, seq, *_) in fastx_read("SRR077487_2.filt.fastq.gz")))'
2386161200

python3 -c '<insert some python code here>' uses the -c option of python3 that tells Python to execute the provided code (instead of starting an interactive interpreter, or running code present in a file).
Here the code has two steps (separated by a ;, in a standard script, we could use a new line instead): first making the fastx_read function available (because it is part of an optional module), second displaying the total length of the sequences.
fastx_read("path_to_a_file") generates (name, sequence, quality) triplets when a fastq file is provided, or (name, sequence) pairs when a fasta file is provided.
The (_, seq, *_) syntax is a form of "tuple unpacking" where we store the second element (the sequence) in a seq variable.
The _ in is just a way to ignore the first element in the triplet/pair. The *_ is a way to ignore whatever is after the second element in the triplet/pair.
len(seq) for (_, seq, *_) in fastx_read("SRR077487_2.filt.fastq.gz") is an expression representing the lengths (obtained with the predefined len function) of the sequences present in the file "SRR077487_2.filt.fastq.gz".
So the code uses the predefined sum function to compute the sum of values that are the lengths of the sequences extracted from the fastq file, and this is displayed using the print function.
This should also work with either fasta or fastq format, and with or without gz compression.
I tested more parsers on a slightly different problem here: https://bioinformatics.stackexchange.com/a/380/292 and you can find inspiration with the other approaches proposed in the other answers there.

seqkit stats fastq.gz. Get SeqKit from https://bioinf.shenwei.me/seqkit/usage/

With bbmap:
reformat.sh </path/to/fastq>


text:
which marker identifications method is recommended when single-nuclei RNA seq datasets are integrated across regions
I'm analyzing Sn-RNA seq datasets which include different brain regions and spinal cord in patients and normal controls. I would have different comparison within and between individuals. I've got a question which I could not find the answer. First, I know that in the clustering step FindConservedMarkers is the best choice when you are comparing controls & patients. I was wondering does that work the same when I want to compare two different regions like the brain and spinal cord? if not what do you recommend for this comparison? as there might be different cell population in two different regions and those regions are not the same

Hi @Paria I can't quite see how this is a question, its either a 'yes'  or a 'no;' (I'd steer towards yes). Maybe if you explain how you arrived at the conclusion that FindConservedMarkers was 'best choice' an informed point of view can be forwarded.

Hi @Michael. I am saying based on the Satija lab workflow. I have single-nuclei 10X genomics datasets. FindConservedMarkers code in Seurat helps to identify the shared gene expression across condition and cell-type identification. In my case I'm comparing two different region within individuls (same condition different regions) and I am asking wether this work the same in my case?


text:
Compute the significance of the overlap between 2 or more gene sets
I was able to compute the significance of the overlap between 2 gene sets using the cdf function of scipy hypergeometric distribution.
I wish to be able to perform the same calculation for more than 2 gene sets; should I use the multivariate hypergeometric distribution cdf function for that?
Are there any websites that provides the same calculations over gene sets so I can validate my results?

Not so trivial, if you want to use a distribution. https://stats.stackexchange.com/questions/52004/overlapping-sets-3-significance-test

Easiest way is to simulate the data, and see how many times you end up with an intersection set larger than what you see

This is a shot at it, first an example dataset:
import matplotlib.pyplot as plt
import numpy as np
import functools
from matplotlib_venn import venn3

# define universe
uni = ["gene"+str(i) for i in range(1000)]
# some overlap
gs1 = uni[250:300] + uni[900:950]
gs2 = uni[:300]
gs3 = uni[250:500]

The ever amazing venn diagram:
venn3([set(gs1),set(gs2),set(gs3)],set_labels=["gs1","gs2","gs3"])

Then a function to draw a set with length equivalent of each set, randomly from the universe and find length of intersection (all 3):
def sim_intersect(uni,set_lengths):
    randomsets = [np.random.choice(uni,n) for n in set_lengths]
    return len(functools.reduce(np.intersect1d,randomsets))

We run this 1000 times:
permuted_values = [sim_intersect(uni,[len(gs1),len(gs2),len(gs3)]) for i in range(1000)]

plt.hist(permuted_values,bins=range(50))

The probability of observing the starting result, using (B+1)/(M+1) as estimator, see this post:
(sum(np.array(permuted_values)>obs_n)+1)/(1000+1)
0.000999000999000999

How would a population size parameter fit in this approach? Doesn't the p value should be affected by the total 'available genes' ?
For example,. if I have two gene sets A, B.
|A| = 250, |B| = 220, |A&B| = 65 As far as I understand,
the p value should change dramatically if the population size is 300 or 30,000.
Could you please explain where/how does it fit in this approach?

Yes you are right, the total available genes would matter. In the code above, you would define the set of available genes as uni . It is not like a hypergeometric where you input a number of genes parameter. You simply simulate the results and see how likely you end up with a result as extreme as yours

easiest thing to do is to re-run the above simulation with a small "universe"


text:
Core genome too short
I'm trying to create a snp distance matrix from M tuberculosis isolates to try and infer transmission networks. My plan was to use snippy to make a core genome alignment and then snp-dists to create the matrix but have hit a snag with creating a core genome alignment. When I try it with all 972 of my samples, snippy halts when trying to run snp-sites, saying Warning: No SNPs were detected so there is nothing to output. The issue is that with all the samples included, the core genome size is 0.
When looking at the alignment statistics provided by snippy, it looks like some isolates have very low numbers of aligned base pairs (??due to low coverage or contamination). I have progressively filtered out samples with low numbers of aligned base pairs and can then get snippy/snp-sites to work. As I increase my filter stringency, the resulting core genome is still short. Eg Filtering out those with aligned base pairs of <90% of the reference had a core genome of only 5747bp (from a reference of 4.4 million bp), while excluding 88 (9%) of my samples.
snippy's author recommends using it's output file core.txt to figure out which samples are the "bad" outliers. That file provides 1) the length of the reference, 2) the number of aligned base pairs, 3) the number of unaligned base pairs, 3) the number of variant sites, 4) the number of heterogenous sites, 5) the number of masked sites, and 6) the number of low coverage sites.
#>   ID                        LENGTH ALIGNED UNALIGNED VARIANT   HET MASKED LOWCOV
#>   <chr>                      <dbl>   <dbl>     <dbl>   <dbl> <dbl>  <dbl>  <dbl>
#> 1 R15795_CATCAAGT_S34_L006 4411532 4192432      6818     755   428 209178   2676
#> 2 R15842_GTCTGTCA_S49_L006 4411532 4187344     10304     784   441 209178   4265
#> 3 R15876_CGCTGATC_S36_L006 4411532 4176662     18992    1292   660 209178   6040
#> 4 R15951_ATTGGCTC_S7_L002  4411532 4170649     14732    1281   980 209178  15993
#> 5 R16019_TGGAACAA_S78_L001 4411532 4190733      8132     712   715 209178   2774
#> 6 R16046_GACTAGTA_S6_L002  4411532 4186069      4309    1257  1084 209178  10892

My question is: what heuristic would you use to filter out samples prior to creating a core genome?
And secondarily, what would a reasonable core genome size be for Mtb?

Targeting Illumina reads with >70X depth-of-coverage is a good place to start.  Once a robust database is established lower coverage isolates can be used effectively too.
More specific to individual SNPs... filtering on VCF QUAL value, Map Quality and/or AC value have shown to be useful.
Core genome will likely be around 1,000 SNPs, but it will depend on the lineage the reference belongs to and the aligned sample lineage.  If both are from the same lineage you're likely to be less than 700 SNPs.  When from different lineages SNP counts will obviously be higher in addition to the alignment error seen.
Pertaining to filtering individual SNPs...
vSNP is a tool developed for SNP analysis.  It can be applied to different specie types but was originally developed for use with TB complex isolates.  Code detail is here.
vSNP is a 2-step process.  First, generating VCF files.  Second, using those VCF files to output SNP tables and trees.  The goal of the script is to achieve a high resolution genotyping analysis that can be validated with minimal effort.  Some steps are included to help achieve this such as reference selection, automatic grouping and position filtering.
Pseudocode:
For set of VCF files created from the same reference...
def find_positions()
    AC=2
    QUAL > 300

For each found position...
def get_snps()
    AC=2 (solid call)
    QUAL > 50
        THEN --> SNP
    AC=1 (mix call)
    QUAL > 150
        THEN --> AMBIGIOUS
    AC=2
    QUAL >= 50
       THEN --> REF
    QUAL < 50
       THEN --> N
    ALT == "None"
       THEN --> -

First find all possible positions that could be informative in all samples/VCF files: find_positions(), then define each position to be included in the alignment: get_snps().
Then it can be convenient to filter consistently poor alignment regions such as PPE/PGRS.


text:
How to get the genomic sequences from a blat result?
Assume we have a query.fa file that contains sequences and we run:
blat -stepSize=5 -repMatch=2253 -minScore=20 -minIdentity=0 -out=pslx /genomes/mm10.fa.qz query.fa output.pslx
the output output.pslx file looks like this:
match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block       blockSizes      qStarts  tStarts
        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
20      0       0       0       0       0       0       0       +       seq     20      0       20      chr9    124595110       44046930        44046950   20,      0,      44046930,       aaaagtatcagtgtgtatag,   aaaagtatcagtgtgtatag,
20      0       0       0       0       0       0       0       +       seq     20      0       20      chr9    124595110       44046930        44046950   20,      0,      44046930,       aaaagtatcagtgtgtatag,   aaaagtatcagtgtgtatag,

What would be a reasonable way to get the genomic contexts (5bp upsteam and 5bp downstream) for each aligned sequence.
For example, assume that blat found that the seq: AAATTGGGGAAAA aligns to chr2:100-113, so the question is how to get chr2:95-118 easily without reinventing the wheel.

I couldn't make it work with bedtools, because my genome's index file is corrupted, but this should work for others who have successfully used bwa or samtools to index their reference genome:
blat -stepSize=5 -repMatch=2253 -minScore=20 -minIdentity=0 -out=pslx /genomes/mm10.fa.qz query.fa output.pslx
awk 'NR>5 {print $14 "\t" $16-10"\t" $17+10}' output.pslx > regions.bed
bedtools getfasta -fi /genomes/mm10.fa.gz -bed regions.bed

https://biopython.org/DIST/docs/api/Bio.SearchIO.BlatIO-module.html

Creat a bed file with chr2 95 118 inside, and then use it with bedtools getfasta to extract your region

@user3479780 I tried that. Add some issue with indexing my reference genome. But this should work too. You are right.

Via BEDOPS convert2bed (psl2bed) and bedops operations:
$ psl2bed < hits.psl | bedops --range 5 --everything - > answer.bed

The file answer.bed will contain target intervals from the PSL (BLAT) input, padded up- and downstream by five bases.
This BED file can be run through samtools faidx or similar to get sequence data.
References:

https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/psl2bed.html

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#range-range

I couldn't get it to work with varied: pcl2bed, samtools, bedtools commands. I think this biopython should work:
import gzip
from Bio import SeqIO
n=5

exact_matches = {}
blat_qresult = SearchIO.read('blat_output.pslx', 'blat-psl', pslx=True)
for hit in blat_qresult:
    for f in hit.fragments:
        if f.hit_span == f.query_span:
            chrom = hit.id
            if chrom not in exact_matches:
                exact_matches[chrom] = []
            exact_matches[chrom].append(f.hit_range)

with gzip.open('/genomes/mm10.fa.gz', "rt") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        if record.id in exact_matches:
            for subseq in exact_matches[record.id]:
                print(record.seq[subseq[0]-n:subseq[1]+n])


text:
DESeq2 - Error in lfcShrink: 'coef' should specify same coefficient as in results 'res'
I have created a DESeq object (dds) and now I want to access the results and apply the lfc shrinkage on them.
res <- results(dds, alpha=0.05)
resultsNames(dds)
[1] "Intercept"                                                  "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages"
[3] "Cell_type_T_cells_vs_Alveolar_macrophages"

Here we see all the coefficients I got from the results function.
Now I try to use the lfcShrink() function on the second coefficient:
res_im_vs_am <- lfcShrink(dds, coef=2, res=res, type="apeglm")
Error in lfcShrink(dds, coef = 2, res = res, type = "apeglm") :
  'coef' should specify same coefficient as in results 'res'

it doesn't work, so I try with the name directly copy-pasted from the above code:
res_im_vs_am <- lfcShrink(dds, coef = "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages", res=res, type="apeglm")
Error in lfcShrink(dds, coef = "Cell_type_Interstitial_macrophages_vs_Alveolar_macrophages",  :
  'coef' should specify same coefficient as in results 'res'

Same error.
According to the ?lfcShrink for coef: the name or number of the coefficient (LFC) to shrink, consult resultsNames(dds) after running DESeq(dds).
So I do not understand why it is not working. Anyone can help?

We use an example dataset:
library(DESeq2)
mat = matrix(rnbinom(9000,mu=100,size=1),ncol=9)
rownames(mat) = paste0("gene",1:1000)
grps = rep(c("A","B","C"),each=3)

dds = DESeqDataSetFromMatrix(mat,data.frame(grps),~grps)
dds = DESeq(dds)

resultsNames(dds)
[1] "Intercept"   "grps_B_vs_A" "grps_C_vs_A"

When you call results, by default it takes the last coefficient, in this example is C vs A:
res <- results(dds, alpha=0.05)

log2 fold change (MLE): grps C vs A
Wald test p-value: grps C vs A
DataFrame with 6 rows and 6 columns
              baseMean     log2FoldChange             lfcSE               stat
             <numeric>          <numeric>         <numeric>          <numeric>
gene1 74.3974631643997 -0.439258650876538  1.22842645044656 -0.357578307367818
gene2 99.4576039995999    1.2903547180366  1.12500005531808   1.14698191518912
gene3 87.1593850381565  -3.11763876685787  1.19848628038255  -2.60131368868297
gene4   51.47821785473  -2.16622171024749 0.953571440023498  -2.27169315200349
gene5 62.4630151830988   1.05879388186169  1.31999477669387  0.802119751195993
gene6 111.967714581788  -1.23916751032406  1.17612862134365  -1.05359863524824

If you want to look at the second comparison:
res <- results(dds, alpha=0.05,name="grps_B_vs_A")
res_lfc <- lfcShrink(dds, coef=2, res=res, type="apeglm")

log2 fold change (MAP): grps B vs A
Wald test p-value: grps B vs A
DataFrame with 6 rows and 5 columns
              baseMean      log2FoldChange             lfcSE            pvalue
             <numeric>           <numeric>         <numeric>         <numeric>
gene1 74.3974631643997 0.00143892076074213 0.275057978370776 0.987946388824352
gene2 99.4576039995999  -0.031292774512063 0.276439771853046 0.632995722892249
gene3 87.1593850381565 0.00806769808245902 0.274672487627629                NA
gene4   51.47821785473 -0.0394804572502463 0.274244572145665 0.609999235800452
gene5 62.4630151830988   0.065875708916621 0.288411208595865 0.234448669110134
gene6 111.967714581788 -0.0871847461207727 0.295733383995144 0.153081913500716

or simply:
res_lfc <- lfcShrink(dds, coef=2, type="apeglm")


text:
How to refine this distribution comparison script
My data look like this

They are compound similarity estimations(the whole file is around 10GB). What I am trying to achieve is to compare the similarity distributions of each compound using the Kolmogorov-Smirnoff test so as to use the data for AIC calculation.
This is my code:
import pandas as pd
from scipy.stats import ks_2samp
import re

with open('ecfp4_file.csv', 'r') as f, open('Metrics.tsv', 'a') as f_out:
    f_out.write('compound_1' + '\t' + 'compound_2' + '\t' + 'Similarity' + '\t' + 'KS Distance' + '\n')
    df = pd.read_csv(f, delimiter = ',', lineterminator = '\n', header = 0)

    for m in range(10, 99):
        dc1 = []
        dc2 = []
        m1 = m/100
        m2 = m/100+0.009
        df2 = df.loc[(df['ECFP4'] >= m1) & (df['ECFP4'] < m2)]          #Work on different thresholds each time

        for i in range(0, df2.shape[0]):
            if len(df2.index[df2['id1'] == df2.iloc[i, 0]].tolist()) >= 5:
                dc1.append(df2.iloc[i, 0])
                dc2.append(df2.iloc[i, 1])
        for i in range(0, len(dc1)):
            dc1_l = []
            dc2_l = []
            df3 = df.loc[df['id1'] == dc1[i]]
            dc1_l = df3.iloc[:, 2].tolist()
            df4 = df.loc[df['id1'] == dc2[i]]
            dc2_l = df4.iloc[:, 2].tolist()
            x1 = re.findall(r"statistic=(.*)\,.*$", str(ks_2samp(dc1_l, dc2_l)))
            f_out.write(str(dc1[i]) + '\t' + str(dc2[i]) + '\t' + str(m1) + '\t' + str(x1[0]) + '\n')

which so far has taken 10 hours to run and it has not yet passed the 0.1 threshold, so it is highly unpractical. Does anyone know of a better way to do this? Or can help me refine my code for faster and better results?
I also tried to refine it and run it in R, but I get a vector memory exhausted (limit reached) error when I run it
install.packages("ggplot2")
library(ggplot2)
PS22_ECFP4 <- read.csv("./PS22_ECFP4.csv")

df <- data.frame(matrix(ncol = 4, nrow = 0))
cols <- c("id1", "id2", "Similarity", "KS Distance")
colnames(df) <- cols

for (i in PS22_ECFP4) {
  sim_df <- PS22_ECFP4[PS22_ECFP4$id1 == i,]
  dist1 <- sim_df[ , 3]
  for (k in 1:nrow(sim_df)) {
    sim <- sim_df[k , 3]
    sim_r <- floor(sim * 100) / 100
    comp2 <- sim_df[k, 2]
    comp2_df <- PS22_ECFP4[PS22_ECFP4$id1==comp2,]
    dist2 <- comp2_df[ , 3]
    if (length(dist1) > 1 & length(dist2) > 1) {
      ks <- ks.test(dist1, dist2)
      df[nrow(df) + 1,] = c(i, comp2, sim_r, ks$statistic)
    }
  }
}

write.csv(df, "./Metrics.csv", row.names = FALSE)

I am not sure how useful this answer is. But may be you could try this.
For your python script:
You use f_out.write within a loop. I think usually this I/O operation is resource intensive. I guess the speed would improve significantly if you store the string in a variable and write the output only once outside the for loop like this:
out_str = ''
for i in range(0, len(dc1)):
            dc1_l = []
            dc2_l = []
            df3 = df.loc[df['id1'] == dc1[i]]
            dc1_l = df3.iloc[:, 2].tolist()
            df4 = df.loc[df['id1'] == dc2[i]]
            dc2_l = df4.iloc[:, 2].tolist()
            x1 = re.findall(r"statistic=(.*)\,.*$", str(ks_2samp(dc1_l, dc2_l)))
            out_str +=  (str(dc1[i]) + '\t' + str(dc2[i]) + '\t' + str(m1) + '\t' + str(x1[0]) + '\n')

f_out.write(out_str)

Regarding the memory error in your R script, are you using 32bit build of R?
Please see this:
http://web.mit.edu/~r/current/arch/i386_linux26/lib/R/library/base/html/Memory-limits.html#
They say "Under most 64-bit versions of Windows the limit for a 32-bit build of R is 4Gb: for the oldest ones it is 2Gb. The limit for a 64-bit build of R (imposed by the OS) is 8Tb."
Conclusion: Please try using your python script, but use the write operation outside the loop.

If I was doing something like this in R (10G file), I'd use data.table
I can't quite work out what you're R code is doing (does it do what you want it to?)
I think the below might do what you want
library(data.table)
dat <- fread("./PS22_ECFP4.csv")
dat[,Similarity:=floor(ECFP4*100)/100]
dat[,'KS Distance':=lapply(1:nrow(dat),function(i)
  ks.test(dat[id1==dat$id1[i],ECFP4],dat[id2==dat$id2[i],ECFP4])$statistic
)]
dat[,ECFP4:=NULL]
fwrite(dat,"./Metrics.csv")

Edit:
The script loops through the data, and for each row creates two subsets of the data matching id1 and id2. Then does the ks.test on the similarity columns of the two subsets.


text:
How do I write tests for a snakemake pipeline?
I've been playing with snakemake, which is really powerful, however, I'd like to create a 'unit test' for the whole pipeline. e.g. a way to call the pipeline that tests that all the rules work together, all the steps complete, etc.
I assumed I could find some guides about how to do this, but I don't find any clear information about how to tests pipelines.
Any advice or suggestions would be welcome.

I made a feature request here: https://github.com/snakemake/snakemake/issues/599

You can run CI on a snakemake workflow when provoding a minimal example dataset. Maybe an example from the snakemake-workflows repo will be helpful: https://github.com/snakemake-workflows/dna-seq-varlociraptor

This is what I do when I want to check whether a program runs from start to finish and produces the expected output on test data. Here I used it to test a C program but the idea can be applied also to a snakemake pipeline as in your case.
The idea is pretty simple really: use python's subprocess module to run the pipeline then check the results using the unittest module (or your favourite unit test package).
For example, save this in test.py:
#!/usr/bin/env python3

import unittest
import os
import subprocess as sp

class Tester(unittest.TestCase):

    def testPipelineMakesStuff(self):
        # Run the pipeline
        p= sp.Popen('snakemake [some parameters]', shell=True, stdout= sp.PIPE, stderr= sp.PIPE)

        # Optionally, get stdout and stderr
        stdout, stderr= p.communicate()

        # Check exits code and other expected output
        self.assertEqual(0, p.returncode)
        self.assertTrue(os.path.isfile('some.pipeline.output.txt'))

        ...Maybe read some.pipeline.output.txt and check it looks ok

    def testPipelineReturnsWithError(self):
        p= sp.Popen('snakemake [some parameters]', shell=True, stdout= sp.PIPE, stderr= sp.PIPE)
        stdout, stderr= p.communicate()
        self.assertTrue(p.returncode != 0)
        self.assertTrue('Something went wrong' in stderr.decode())

    def ...more test cases...

if __name__ == '__main__':
    unittest.main()

Run this script as:
 python test.py

and check all tests pass ok. Of course, you can embed python test.py in travis or other continuous integration system as I've done here.