Why do ten rows (Figure_1) correspond to 2 bits (Figure_2) in a sequence logo? Following this question, I'm confused with the computation of sequence logo

Following data comes from the book "Machine Learning - A Probabilistic Perspective (Figure_1)"

here is the corresponding sequence logo (Figure_2).

ten rows represent sequences of DNA (e.g. row 1 could be a human sequence, row 2 could be the equivalent mouse sequence etc.)

Each column corresponds to a particular position.

Why do ten rows (Figure_1) correspond to 2 bits (Figure_2)?

Bits are not frequencies. If a position only contains an A (position 3 for example) then you would need 2 questions (bits) to derive that information without a priori knowledge.

Is it a G or C? if no > then it is a A or T
Is it T? If no then it is an A

Position 1 can be derived by 1 question only:

is it a G or C?  If no then it is an A or T

In this case the number of bits can be calculated using the frequenci)es.

For 4 bases, A only (100%):

4 * -0.25 * log2(0.25) - (-1 * log2(1)) = 2 (bits)

For 8 letters, A only (100%):

8 * -1/8 * log(1/8)- (-1 * log2(1)) = 3 (bits)

Mixed base possibility (A/T 50% each) with 4 bases:

4 * -0.25 * log2(0.25) - (-0.5 * log2(0.5)+-0.5 * log2(0.5)) = 1 (bit)

As mentioned in the comments, you can follow the wiki link for an example.

Or, an example with a little more explanation from the biology SE site.

sequence logo calcuation

Hopes this helps


Can I model technical replicates in DESeq2? I’d normally use collapseReplicates (or do the collapsing upstream) to handle technical replicates.

However, in my current RNA-seq experimental design, samples were sequenced twice using different library preparation protocols, leading to marked differences in the resulting count estimates (in fact, the choice of different protocol explains most of the variance in the bulk data according to PCA). I would therefore like to model these differences by adding them as a covariate to the DESeq2 model.

I thought I could just add another factor to the design formula, equivalent to the “pasilla” example in the DESeq2 vignette, which adds the factor type for the library type (single-end vs paired-end). However, the pasilla vignette states that these different libraries are actually independent biological replicates, not technical replicates as I had previously assumed.

I’m now concerned that having technical replicates in the design might artificially inflate the power. Can the design be salvaged?

It’s worth noting that we also have biological replicates, and a full rank design; i.e. every combination of treatment and library prep, with biological replicates.

Practically speaking, there's no way to include the technical replicates in that design (in DESeq2 at least). Your concern regarding inflating the power is exactly correct and the only way to combat that would be to add a pairing factor like one might do with case-control or tumor-normal studies. That is, something like:

  group libraryPrep sample
1    WT           A      1
2    WT           B      1
3   MUT           A      2
4   MUT           B      2
5    WT           A      3
6    WT           B      3
7   MUT           A      4
8   MUT           B      4

Here sample pairs the technical replicates together. However this ends up either being rank deficient or in the end not more informative.


RunUMAP in Seurat not working: module 'umap' has no attribute 'UMAP' I am trying to run UMAP in the following way:

RunUMAP(seurat_object, dims.use = 1:15, reduction.use = "pca", reduction.name = "umap",
    reduction.key = "UMAP", n_neighbors = 30L, min_dist = 0.3, metric = "correlation")

Or just this:

RunUMAP(seurat_object, dims.use = 1:15, reduction.use = "pca")

It is giving me the error:

  Error in py_get_attr_impl(x, name, silent) :
    AttributeError: module 'umap' has no attribute 'UMAP'

I already installed umap-learn successfully in this way: conda install -c conda-forge umap-learn and checked that Rstudio is using the same python version as is provided by default in the terminal: 3.6.5

This is resolved by installing the more recent version of the code from GitHub.

git clone https://github.com/lmcinnes/umap
cd umap
pip install --user -r requirements.txt
python setup.py install --user

As Peter said, this kind of question is best directed to the Issues page on their GitHub repository.


Error in colSums(assay(sce)) : 'x' must be an array of at least two dimensions I am having a SingleCellExperiment (sce) object:

  sce <- SingleCellExperiment(list(counts=UMI_count))

which has two dimensions:

dim(sce)
[1] 46761 68027

When I am trying to get column sums for all rows the following does not work:

colSums(assay(sce), na.rm = FALSE, dims = 1)

Giving me an error:

  Error in colSums(assay(sce)) : 'x' must be an array of at least two dimensions

Apparently the issue was already described:

https://support.bioconductor.org/p/100954/

However, I can not find the solution. Is there another way of getting sums of columns of either SingleCellExperiment object or just dgcMatrix?

Ultimately I just implemented it manually with the precursor to SingleCellExperiment matrix (UMI_count which is dgcMatrix type) object in the following way:

# Removing genes with 0 - counts in all of the columns
num_removed = 0
for(i in 1:nrow(UMI_count)) {
  idx = i - num_removed
  if(sum(UMI_count[idx,]) < 1) {
    UMI_count <- UMI_count[-idx,]
    num_removed = num_removed + 1
    print(i)
  }
}

sce <- SingleCellExperiment(list(counts=UMI_count))

I am getting sums of the columns for each row and then filtering out rows with the sum < 1. So, sum function for the slices - rows (numeric type) is working.


Transcribing DNA exons and append transcripts into a list I am trying to find a way to read and transcribe a list of DNA sequences (list of lists) only when the for loop finds a start codon (triplet of the list items) and until it finds a stop codon, over and over, but the if function will only transcribe the first codon and won't go through the whole sequence.
I'm trying to do this with a mix of normal python and biopython functions cause I think it's simpler.

This is an example:
Start codon: ATG
Stop codon: TAG
gene_1= 'ACGGACTATTC'
gene_2= 'GGCCATGAGTAACGCATAGGGCCC
gene_3=GGGCCCATGACGTACTAGGGGCCCATGCATTCATAG
So, the first gene does not contain any transcript, whereas the second contains 1 and the third contains 2. I'm trying to get rid of everything outside these reading frames and append these transcripts into a list that should look something like this.
mRNAs = ['AGUAACGCA', 'ACGUAC', 'CAUUCA']
This is what I have so far:
genelist=[]

gene_1= 'A','C','G','G','A','C','T','A','T','T','C'
gene_2= 'G','G','C','C','A','T','G','A','G','T','A','A','C','G','C','A','T','A','G','G','G','C','C','C'
gene_3='G','G','G','C','C','C','A','T','G','A','C','G','T','A','C','T','A','G','G','G','G','C','C','C','A','T','G','C','A','T','T','C','A','T','A','G'

genelist.append(gene_1)
genelist.append(gene_2)
genelist.append(gene_3)

def transcription(ORF):
    mRNA= ''
    for i in range(0, len(ORF), 3):
        codon= ORF[i:i+3]
        if codon != 'ATG':
            next(codon)
            if codon == 'ATG':
                mRNA=codon.transcribe()
                if codon == 'TAG':
                    break
    return(mRNA)

mRNAs=[]
for gene in genelist:
    for codon in gene:
        mRNA= transcription(codon)
        mRNAs.append(mRNA)
print(mRNAs)

But it is not really giving anything back, I wonder if the code it's too redundant and I really don't need to define a function here, do you know any better way to do this?
Thaaanks!!

I'm afraid none of your code makes much sense. You have written a function that expects an array, and then you are passing it a single letter
for gene in genelist:
    for codon in gene:
        mRNA= transcription(codon)

Then, the function also doesn't make sense. What does next(codon) mean? How would your script know what a codon is or how to move to the next one? And your if codon == 'ATG': is within the if codon != 'ATG': block, so it will never be executed by definition.
I think what you are looking for is something like this, using the same basic approach as your script:
#!/bin/python3
from Bio.Seq import Seq
from Bio.Seq import transcribe

genelist=[]

gene_1= 'A','C','G','G','A','C','T','A','T','T','C'
gene_2= 'G','G','C','C','A','T','G','A','G','T','A','A','C','G','C','A','T','A','G','G','G','C','C','C'
gene_3='G','G','G','C','C','C','A','T','G','A','C','G','T','A','C','T','A','G','G','G','G','C','C','C','A','T','G','C','A','T','T','C','A','T','A','G'

genelist.append(gene_1)
genelist.append(gene_2)
genelist.append(gene_3)

def transcription(ORF):
    mRNA= ''
    foundStart = False
    foundEnd = False
    for i in range(0, len(ORF), 3):
        codon= "".join(ORF[i:i+3])
        if codon == 'ATG' and not foundStart:
            foundStart = True
        if foundStart and not foundEnd:
            cc=transcribe(codon)
            mRNA = mRNA + transcribe(codon)
        if codon == 'TAG':
            foundEnd = True

    return(mRNA)

mRNAs=[]
for gene in genelist:
    mRNA = transcription(gene)
    mRNAs.append(mRNA)
print(mRNAs)

Which returns:
$ foo.py
['', '', 'AUGACGUACUAG']

Of course, it only returns one mRNA since you only have one gene with an ORF in the 1st frame and your script was only looking at the forward strand and the first frame. In reality, ORFs can be in any of the 6 possible frames. You are also assuming that all ORFs will end with a TAG which is perhaps the case in the species you are working in? I don't remember if there is any genetic code that has TAG as the only stop.
If you do actually need to find ORFs in any frame, then I hope the script above should be enough to get you started.


Unable to install bedtools on windows 10 ubuntu I am trying to install bedtools on windows 10, but I get an error I don't understand. How can I fix it?

Building BEDTools:
=========================================================
DETECTED_VERSION = v2.27.1
CURRENT_VERSION  =
Updating version file.
 * Creating BamTools API
- Building in src/utils/FileRecordTools
  * compiling FileRecordMgr.cpp
make[1]: g++: Command not found
Makefile:36: recipe for target '../../../obj//FileRecordMgr.o' failed
make[1]: *** [../../../obj//FileRecordMgr.o] Error 127
Makefile:186: recipe for target 'src/utils/FileRecordTools' failed
make: *** [src/utils/FileRecordTools] Error 2

My personal advice would be to avoid compiling tools like this if you can avoid it, especially where you are using a "barebones" setup like the one you get with WSL Ubunutu as you'll almost always run into dependency issues.

The linux version of conda installs just fine in WSL and bioconda hae pretty much any bioinformatics package you might need, already compiled.


Statistical test design question I have condition as such mouse liver and spleen data.
So in liver I have WT (Wild Type) and KO (knock out) as well as in spleen WT (Wild Type) and KO (knock out).

I have to compare the parasite load in WT and KO condition in both liver and spleen.

In terms of data points what I got is this:

10 mice for WT, 6 mice for KO in case of liver
10 mice for WT, 5 mice for KO in case of spleen.

WT 20 8 21 9 22 10 23 11 24
KO 12 25 26 27 28

So far what I know is I have to perform a parametric test.

But I'm not sure if I have to go for a paired t-test or an unpaired t-test.

Because here I'm comparing with between two groups, only so I don't have to go for anova if I'm not wrong. Secondly there is a difference in number of observation that is available for WT as compared to KO. So would that affect my choice of test as well.

Suggestion would be highly appreciated

So far what i know is i have to perform a parametric test .

How did you determine this?

  But im not sure if I have to go for a

  Paired t-test or Unpaired t-test

If you can't tell the difference between when to do a paired t-test and an unpaired one, you need to stop and talk to someone who can teach you basic statistics.  Throwing statistical tests at numbers when you don't understand them is a recipe for disaster.


RNAseq biological replicates not clustering in PCA plots I have RNAseq data from 4 samples with 3 biological replicates per sample. I am currently trying to do the differential expression analysis with DESeq2 but the biological replicates will not cluster together when I make the PCA plot or correlation heatmap. This is my first time with RNASeq analysis and so am not sure what the best route forward is? I would like to avoid repeating the experiment with new samples if possible!
My pipeline prior to DESeq2 was the following:
FastQC quality check -> Trimmomatic -> Kallisto
I used tximport to convert kallisto files into suitable format for DESeq2

PC1 is 81% of the variance?
This PCA plot confirms that different tissues are different.  You probably already knew that.  I'd make more PCA plots that are tissue specific.  That will be more informative than these.
Personally, I'd also not do DESeq on all of these samples together, unless the goal of your experiment really is to learn the differences between these two tissues.


seqret: Warning: bad /protein_id value I am using the seqret tool from the Emboss suite to transform produce an EMBL annotation file from a GTF annotation file and the corresponding FASTA sequence:

seqret -sequence ecoli.fa -feature -fformat gff -fopenfile ecoli.gtf -osformat embl -auto -outseq ecoli.embl

It produces my desired EMBL file but it warns against a 'bad protein_id':

...
Warning: bad /protein_id value '"CCE57732"'
Warning: bad /protein_id value '"CCE57733"'
Warning: bad /protein_id value '"CCE57734"'
Warning: bad /protein_id value '"CCE57735"'
Warning: bad /protein_id value '"CCE57754"'

What's wrong with these protein IDs? The organism I am working with is Escherichia coli.

Here are a few lines from my GTF file. The first line has a protein_id.

Chromosome      ena     CDS     349     2808    .       +       0       gene_id "HUS2011_0001"; transcript_id "CCQ26880"; exon_number "1"; gene_name "thrA"; gene_source "ena"; gene_biotype "protein_coding"; transcript_name "thrA-1"; transcript_source "ena"; transcript_biotype "protein_coding"; protein_id "CCQ26880";
Chromosome      ena     start_codon     349     351     .       +       0       gene_id "HUS2011_0001"; transcript_id "CCQ26880"; exon_number "1"; gene_name "thrA"; gene_source "ena"; gene_biotype "protein_coding"; transcript_name "thrA-1"; transcript_source "ena"; transcript_biotype "protein_coding";
Chromosome      ena     stop_codon      2809    2811    .       +       0       gene_id "HUS2011_0001"; transcript_id "CCQ26880"; exon_number "1"; gene_name "thrA"; gene_source "ena"; gene_biotype "protein_coding"; transcript_name "thrA-1"; transcript_source "ena"; transcript_biotype "protein_coding";

You're supplying a GTF feature file, but telling it that the features are in GFF3 format. I think you need -fformat gff2. See https://www.ebi.ac.uk/Tools/sfc/emboss_seqret/help/


Clarification on Gene Enrichment When I run a GSEA analysis on two conditions from the same RNaseq (negative control PBS injection VS positive control CpG injection) from the same dataset/same gene list, I get results that look something like this:

Notice in my example that many gene sets are significantly enriched in PBS VS CpG, but none are very significantly enriched in the inverse, CpG VS PBS.

I have a basic question:

If we're comparing two items, say A and B, shouldn't a up-regulated gene set that's enriched in A have a corresponding down-regulated gene set that's enriched in B?

I'm consistently confused about how to interpret the fact that one of my two conditions has many more enriched gene sets than the other. What am I missing/misunderstanding?

Thanks!

If you compare A vs B the genes's fold change will have the opposite sign to B vs A. So will be the gene set up or down-regulated depending on the comparison to take.

The gene set test analyze if a given group of elements is sorted in a certain way in a list (I am talking about the GSE like the one you performed or the one in the Broad Institute). Usually it is used to say that gene set X is up-regulated in the comparison A vs B, that is, gene set X (in general) is more expressed in A than in B. The enrichment is measured in a enrichment score (or a normalized enrichment score, the higher, the clearer trend/distribution group X has. To asses X is really differentially expressed a p-value is also calculated.

Despite that the number of the up-regulated gene sets sum up all the gene sets tested as pointed by juod, the number of gene sets with a a p-value below 5% and 1% is not the same in both comparisons which indicates that you are not really using the same samples in both gene sets.

Check how you performed your comparison. Are the genes in each comparison with the same value and opposite sign?
Check also the software used to run the GSEA analysis (if I recongize the type of report is from a package in Bioconductor [I don't recognize which one now]) as you (or me) might have misunderstood something or there is an error if the comparisons are done properly. For instance it could be that the report of p-values below the threshold would be from only the up-regulated gene sets.


Holistic enzyme activity determination with computation I'm currently working on a program which will determine the activity difference between an original enzyme and a variant with just 1 or 2 variants. Of course, I'm not talking about the "kinetic activity", but some indications of a probable activity difference should be given. The system has to be sensitive enough to determine the probable activity difference between unthought-of variants.
From the information that I've got from my research, there is two option. I can use one of 3 quantitive enzyme reaction analysis tools (such as EVB, FEB and LRA) or their hybrids.
EVB => Empirical Valence Bond
FEB => Free Energy Perturbation
LRA => Linear Response Approximation
Source:
https://www.researchgate.net/post/Can-the-activity-of-proteins-be-studied-using-molecular-dynamic-simulations
Now, my questions are: - Does these application be used for such purpose? - For an accurate prediction of  enzyme activity difference, are these programs good enough or require developement? - Which one of them can give the best result? - And lastly, are QM/MM methods a must for an accurate answer?
I'm a high school student so sorry if my questions sound dumb.
If you can answer my inquiry's I would appreciate it:)
Thanks

Your title says holistic. This is a tad problematic as there's layers upon layers. Say, post-translation regulation, inhibiting metabolites, interacting protein etc. That is why when talking of an enzyme inhibitor, at the biochemical level one speaks in terms of k_i (inhibition constant), while at the cellular level ("holistic") one talks of IC50.
This could be a biocatalysis/synthetic biology problem (make an enzyme better) or translational medicine problem (why is this mutation pathogenic).
In terms of in silico enzyme design, one can have all sorts of issues that need to be addressed:

catalysis
side reactions
stability
aggregation
solubility

Also there's the expression problems, which may be due to subpar promoter, insufficient cofactors and poor codon usage.
For clinical variant assessment, a variant may:

destabilise leaving to degradation
remove a post-translational modification
disrupt the active site, affecting catalysis
lose a binding partner
result in incorrect localisation
result in solubility issues

So force-field scores are not the sole solution. Furthermore, a forcefield method gives one a ∆∆G (difference in relative Gibbs free energy), not a K_M (Michaelis constant) and a k_cat (catalytic rate), although it is somewhat possible to calculate these, or a relative variant of these.

a variant with just 1 or 2 variants

A variant protein can have one or two mutations (e.g. G10R) in its sequence. A mutant and a variant are near synonyms, but variant is generally intended to be a mutant in a population.

Does these application be used for such purpose?

For enzyme design MD simulations are not generally the first step. Rosetta pmut_scan for example is a nice starting step as the mutational landscape of a protein is N^19, while a MD simulation takes weeks and 6+ core per single variant.
This calculates only an increase in stability (relative difference in Gibbs folding energy), not reactivity. For that you can use a transition state, whose results are inaccurate, but majorly affected by the set-up, i.e. junk-in junk-out plays a big part. There's a saying 'a week in the lab saves you an hour at the library': this is backwards for sarcasm as it happens too often, namely that is background research is always a lifesaver. An QM-MD simulation on the wild type reaction definitely can help in elucidating things if it is not clear.
Online tools
Generally, however, there are resources out there to do many stability calculations already.
SynBio If you are interested in picking a residue or two to make an enzyme more stable, without resorting to a thermophilic enzyme, PROSS is a very nice tool.
Clinical (Caveat: If your variant is a variant that affects you, a friend or a relative, you/they should speak to a genetic councillor and not do research on it yourself) There are several tools that can describe a variant at the protein level. Do remember that ∆∆G is a potential so more negative is more stable —except in SDM where the score is inverted.

Missense3D
Venus
Miscast
SDM


Searching tool to calculate phase/switch error rate I'm looking for a tool which, given a truth vcf file and a test vcf file, calculates the phase/switch error rate. I performed phasing of a vcf using WhatsHap and want to compare the outcome to some ground truth phased vcf I have. I can't find a tool, and I don't want to write it myself. It seems like a rather common thing... so I'd expect it to exist, but my google-fu is weak today.

Example of test lines:

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample

chr10   100588  .       G       A       .       PASS    AF=0.18;AC=1;NS=150;AN=2;EAS_AF=0.14;EUR_AF=0.39;AFR_AF=0.08;AMR_AF=0.31;SAS_AF=0.2;VT=SNP;DP=11907     GT:PS   1|0:48005
chr10   102385  .       A       G       .       PASS    AF=0.07;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.1;AMR_AF=0.19;SAS_AF=0.01;VT=SNP;DP=22665      GT:PS   0|1:48005
chr10   105170  .       G       A       .       PASS    AF=0.03;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.0;AMR_AF=0.19;SAS_AF=0.0;VT=SNP;DP=17140       GT:PS   0|1:48005
chr10   105365  .       G       T       .       PASS    AF=0.19;AC=1;NS=150;AN=2;EAS_AF=0.14;EUR_AF=0.39;AFR_AF=0.09;AMR_AF=0.31;SAS_AF=0.2;VT=SNP;DP=21845     GT:PS   1|0:48005
chr10   106057  .       T       C       .       PASS    AF=0.13;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.32;AMR_AF=0.19;SAS_AF=0.01;VT=SNP;DP=22559     GT:PS   0|1:48005

Example of truth lines

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample
chr10   100554  .       C       CT,CTT  30      .       .       GT:AD   1|2:0,1,1
chr10   100588  .       G       A       30      .       .       GT:AD   0|1:1,1
chr10   102385  .       A       G       30      .       .       GT:AD   1|0:1,1
chr10   102636  .       T       C       30      .       .       GT:AD   1|1:0,2
chr10   102757  .       T       C       30      .       .       GT:AD   0|1:1,1
chr10   104815  .       GCC     G,GC    30      .       .       GT:AD   2|1:0,1,1
chr10   105170  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   105365  .       G       T       30      .       .       GT:AD   0|1:1,1
chr10   106057  .       T       C       30      .       .       GT:AD   1|0:1,1
chr10   106110  .       C       G       30      .       .       GT:AD   1|1:0,2
chr10   106261  .       A       G       30      .       .       GT:AD   1|0:1,1
chr10   108612  .       T       G       30      .       .       GT:AD   1|0:1,1
chr10   108646  .       A       AT      30      .       .       GT:AD   0|1:1,1
chr10   108834  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   110840  .       C       T       30      .       .       GT:AD   1|1:0,2
chr10   111743  .       A       AT      30      .       .       GT:AD   1|0:1,1
chr10   112016  .       T       A       30      .       .       GT:AD   1|0:1,1
chr10   112262  .       C       T       30      .       .       GT:AD   1|1:0,2
chr10   113006  .       C       T       30      .       .       GT:AD   0|1:1,1
chr10   113031  .       G       A       30      .       .       GT:AD   0|1:1,1
chr10   113136  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   113359  .       A       T       30      .       .       GT:AD   1|1:0,2
chr10   113583  .       C       G       30      .       .       GT:AD   0|1:1,1
chr10   113997  .       G       C       30      .       .       GT:AD   1|1:0,2

Assuming your files are as you show, and the genotype field is always the 1st entry of the sample information field (the 10th), then you can just do this directly with a simple awk script:

$ awk -F"\t" -vOFS="\t" 'BEGIN{print "#CHROM","POS","ID","REF","ALT","Truth","Test"}{ if(/^[^#]/){var=$1$2$4$5; gt=substr($10,1,3); if(NR==FNR){a[var]=gt; next;} if(var in a && a[var]!=gt){print $1,$2,$3,$4,$5,gt,a[var]}}}' truth.vcf test.vcf
#CHROM  POS     ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0

Here's the script split into separate lines for clarity (you can still copy paste directly into a terminal):

awk -F"\t" -vOFS="\t" '
    BEGIN{
        print "#CHROM","POS","ID","REF","ALT","Truth","Test"
    }
    {
        if(/^[^#]/){
            var=$1$2$4$5;
            gt=substr($10,1,3);
            if(NR==FNR){
                a[var]=gt;
                next;
            }
            if(var in a && a[var]!=gt){
                print $1,$2,$3,$4,$5,gt,a[var]
            }
        }
    }' truth.vcf test.vcf

The basic idea is:

awk -F"\t" -vOFS="\t": set the input and output field separators to a tab.
BEGIN{print "#CHROM ... }: print a header.
if(/^[^#]/){: skip the header lines of the input vcf files
var=$1$2$4$5: get a "signature" for this variant, by concatenating the chromosome, position, ref and alt fields.
gt=substr($10,1,3); : get the genotype value.
if(NR==FNR){: if this is the 1st file/
a[var]=gt; next: save the genotype for this variant and skip to the next line.
if(var in a && a[var]!=gt){: this part will only be run for the 2nd file because of the 'next' in the previous block. If this variant was in the truth set and if its genotype there isn't the same as in the test set, print the variant and the two genotypes.

If you need the error rate, you can extend the script thusly:

awk -F"\t" -vOFS="\t" '
BEGIN{
    print "#CHROM","POS","ID","REF","ALT","Truth","Test"
}
{
    if(/^[^#]/){
        var=$1$2$4$5;
        gt=substr($10,1,3);
        if(NR==FNR){
            a[var]=gt;
            next;
        }
        tot++;
        if(var in a && a[var]!=gt){
            bad++;
            print $1,$2,$3,$4,$5,gt,a[var]
        }
    }
}
END{
    printf "Correct: %s of %s (%s%)\n",tot-bad,tot,(tot-bad)*100/tot;
    printf "Wrong: %s of %s (%s%)\n",bad,tot,bad*100/tot;
}' truth.vcf test.vcf

On your example dataset, that prints:

#CHROM  POS ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0
Correct: 0 of 5 (0%)
Wrong: 5 of 5 (100%)


How to uninstall hmmpress from my conda environment? I have installed hmmpress on my Mac OS X. It is here:

which hmmpress
/Users/user/anaconda/bin/hmmpress

It is available in my root environment:

hmmpress -h
# hmmpress :: prepare an HMM database for faster hmmscan searches
# HMMER 3.1b2 (February 2015); http://hmmer.org/
# Copyright (C) 2015 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Usage: hmmpress [-options] <hmmfile>

Options:
  -h : show brief help on version and usage
  -f : force: overwrite any previous pressed files

I don't remember exactly how I installed it but I thought that I used conda because it is in the anaconda folder. However, conda search hmmpress returns no packages. I thought it could be part of the hmmer package. Indeed, the hmmer package exists in the bioconda channel. However, if I try conda uninstall hmmer, it says:

conda uninstall hmmer
Solving environment: failed

PackagesNotFoundError: The following packages are missing from the target environment:
  - hmmer

Does anyone know how I can remove hmmpress? I thought of rm -rf the /Users/user/anaconda/bin/hmmpress but it does not seem like a clean solution. There might be dependencies elsewhere that I could also get rid of.

This time conda uninstall hmmer in the root environment did the job. The sequence of commands that led me there is shown below. Thanks Devon Ryan.

conda list | grep hmmer
hmmer                     3.1b2                         3    bioconda

conda uninstall hmmer
Solving environment: done

## Package Plan ##

  environment location: /Users/user/anaconda

  removed specs:
    - hmmer

The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    conda-4.5.9                |           py36_0         626 KB  conda-forge

The following packages will be REMOVED:

    barrnap: 0.9-0        bioconda
    busco:   3.0.2-py36_6 bioconda
    hmmer:   3.1b2-3      bioconda
    prokka:  1.12-3       bioconda

The following packages will be UPDATED:

    conda:   4.5.8-py36_1 conda-forge --> 4.5.9-py36_0      conda-forge
    openssl: 1.0.2o-0     conda-forge --> 1.0.2o-h470a237_1 conda-forge

Proceed ([y]/n)? y

Downloading and Extracting Packages
conda-4.5.9          |  626 KB | ################################################################################################################# | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

Checking that hmmpress is not available anymore:

hmmpress
zsh: command not found: hmmpress


Annovar index out of date I'm trying to use Annovar to annotate some variants with their CADD and FATHMM scores. I've downloaded the latest versions of the software and the databases but when I run it I get an error saying the index is out of date.

  WARNING: Your index file annovar_latest/humandb/hg19_cadd.txt.idx is out of date and will not be used. ANNOVAR can still generate correct results without index file.

It still runs without the index file but it is extremely slow. How can I update the index file?

Your safest bet is to just redownload the annotation. I've not seen official documentation of the index used by annovar, but apparently it's just a text file with a single line header and chromosome-bin coordinates. Over on SEQanswers, there's a thread discussing this issue with a single perl script that allegedly recreates the index.


Parsing SAM/BAM files for additional information I used BWA-MEM to alignment and I would like to gather the some informations like total % of match, mismatch, insert/delete etc. I am wondering if there is any existing tools that produces this information?

Other options... samtools stats maybe gives you some of the stats you are after. Chances are you already have samtools so nothing additional to install. Also useful may be pysamstats, in particular the variation sub-command.


Does the DNA or RNA of some ethnic groups map better than others' to the human reference genome sequence? I believe that the human genome reference took DNA samples from different people and that some natural variation is included in extra contigs.

However, the human reference genome comes from a limited number of people.

Can the limited number of donors and the linear representation of the genome (as opposed to a graph representation that could potentially include variation) bias an analysis in favour of one ethnic group over another? Are there studies that have compared the alignment rates/mappabilities of reads of people coming from different geographic regions? Do some ethnic groups' reads align better to the reference?

For simple variants like SNPs it would not really be a problem to use the current genome assembly for other ethnic groups. But for more complex variants this could be indeed problematic, however not only for ethnic groups but also for individuals within the same population. Think of very complex regions such as HLA or KIR.

In studies where they compare different ethnic groups (e.g., 1000 G) they rather not discuss such matters, but act like the analysis is 100% perfect.

There is a nice comment in Genome Biology about this issue.


Pathifier bioconductor package recommends to use technical noise for an argument. How to determine it? I'm running the pathifier approach against C2 pathway curated
database for a specific microarray dataset.
As I was reading the documentation of the pathifier, in order to configure it properly for my dataset, I saw that for the min_std argument they suggest to use the technical noise.

min_std:   The minimal allowed standard deviation of each gene. Genes
with lower standard deviation are divided by min_std instead of their
actual standard deviation. (Recommended: set min_std to be the
technical noise).

But how am I going to calculate it ? At first run I just used the value 0.2254005 from the example they provide with that package, but this might not be the right for my dataset.
Update:
So from here it says: The minimal standard error (min_std) was set as the first quartile of the standard deviation in the data. How am I supposed to calculate this first quartile of the std? Any idea? Should I calculate the sd() for each gene and then find the first quartile?
What do you think about this function?

Answer from @llrs, converted from comment:
Probably is a question for the maintainer but I would guess that is a variance of the gene in the control dataset. So, if in normal samples the gene has a standard deviation of 0.5 this is the expected too in the tumoral cells. calculate the sd only on the control samples. So if the healthy/control samples have such variation the altered ones should have these or more according to the article (if I understand it correctly).
std <- apply(data[, controls], 1, sd, na.rm = TRUE); quantile(std, 0.25) should be much faster and efficient
But  I would check it in support.bioconductor with the authors/maintainers.


What is the exact definition of "breakend"? I'm reading the Manta User Guide
→ it is a method to discover structural variants and indels from next generation sequencing data.
The word 'breakend' is used in several contexts
I have searched for a definition in duckduckgo, qwant and google but have not found a concise definition.
Thanks for you help

A structural variant is an inversion, deletion, duplication, or translocation. A "breakend" is one endpoint of a structural variant. This paper (Alkan et al.  Nature Review Genetics 2011) may help you understand: https://www.nature.com/articles/nrg2958
The Manta paper and documentation seem sparse on diagrams. DELLY, another structural variant caller, has diagrams in their publication. Their paper may be more of use: https://academic.oup.com/bioinformatics/article/28/18/i333/245403
Specific example: Say you're sequencing a chronic myeloid leukemia (CML) tumor. You will likely find reads that half map to the ABL1 gene on chr9 and half map to the BCR gene on chr22 (https://en.wikipedia.org/wiki/Philadelphia_chromosome). The genomic position of the location where chr9 turns to chr22 is called the breakpoint. The positions on chr9 and chr22 are the breakends.


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'.


Definition of "seed" in sequence alignment I would like to know what is meant by "seed" for various sequence aligners. How is it important?

The seed is the subset of a read used in the first step of an alignment. Many aligners work by a seed-and-extend model, wherein they first find all regions matching the "seed" and then extend the alignment around that allowing mistmatches and indels until it either gives up (and therefore uses a different seed) or finds a sufficiently good alignment.


How to fix RepeatMasker.lib.nsq missing error when running RepeatModelor? I have run RepeatModelor and I don't see the classified output file. I got this error.
Missing /home/mudith/.conda/envs/braker/share/RepeatMasker/Libraries/RepeatMasker.lib.nsq!
Please rerun the configure program in the RepeatModeler directory before running this script.

I think I have to run configure file. I found it. But it doesn't run when I just try
perl ./configure
Can someone please explain me how to run this? I'm new to this.
Also I found this in github. But I don't understand.
cd $(dirname $(which RepeatMasker))/../share/RepeatMasker
# ./configure downloads required databases
echo -e "\n2\n$(dirname $(which rmblastn))\n\n5\n" > tmp && ./configure < tmp

I don't understand the third line.

What that code snippet is doing:
That snippet you found on github (where is that from exactly?) looks like it tries to run RepeatMasker's configure script, feeding it some options that are needed by RepeatModeler.  The cryptic third line generates and dumps a block of text to a file named tmp:

2
path/to/bin/rmblastn

5

And then passes that in to RepeatMasker's own configure script as though it were entered on the keyboard.  It's the equivalent of hitting enter to accept a default, entering "2" to select the RMBlast option, typing in the path to rmblastn, hitting enter to set RMBlast as the default, and then typing "5" to save and quit.
That all looks fine, but I think this might be a problem with your (and mine, see below) RepeatMasker install rather than RepeatModeler itself.  These tools are new to me but I'm interesting in seeing this work too, so I gave a shot at troubleshooting it here.
The way I was able to get it working (note I assume you're working with conda packages here) was:

pushd $CONDA_PREFIX/share/RepeatMasker
./util/buildRMLibFromEMBL.pl Libraries/RepeatMaskerLib.embl > Libraries/RepeatMasker.lib to "automatically create the old RepeatMasker.lib FASTA style database" as the script's comments put it.  NOTE: you might want to first delete or move the old empty file rather than alter it because of how conda shares data between environments.  If you update it in place (which > will still do), you'll alter the file for all conda environments using the package, and conda may complaint about unexpected file sizes when checking packages.
perl configure to do RepeatMasker's setup including some makeblastdb commands using RepeatMasker.lib.  (Note this is RepeatMasker's configure script and not RepeatModeler's.) You can type in the options manually or you could probably use that code snippet you found.
popd and run RepeatModeler/RepeatClassifier and hopefully it finds the BLAST library files now

All the gritty detail:
I see a conda path in your logs.  Are you using the bioconda packages named repeatmodeler and repeatmasker?  It looks to me like the repeatmasker package in particular includes an empty RepeatMasker.lib file and that's the cause of the trouble.  I don't see any way that RepeatMasker.lib is set up automatically via the configure script, and the BLAST database files like RepeatMasker.lib.nsq that RepeatModeler expects are built from that (empty) file.
The Please rerun the configure program in the RepeatModeler directory message is a red herring, I think, in this case because the $CONDA_PREFIX/share/RepeatMasker/Libraries/RepeatMasker.lib file does exist, just not the other files for a database it assumes should exist.  If I instead run configure for RepeatMasker and update its script to not hide error messages, I see this coming from makeblastdb:
Building a new DB, current time: 10/26/2020 16:10:07
New DB name:   /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib
New DB title:  /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib
Sequence type: Nucleotide
Keep MBits: T
Maximum file size: 1000000000B
BLAST options error: File /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib is empty

And it's true, that's an empty file provided by the conda package (the metadata in $CONDA_PREFIX/conda-meta/repeatmasker-4.1.0-pl526_0.json shows size_in_bytes: 0 for that one so that's apparently "correct" for the package.)
I saw a reference to that filename in a buildRMLibFromEMBL.pl script in the RepeatMasker share directory, and running that produces a FASTA file:
./util/buildRMLibFromEMBL.pl Libraries/RepeatMaskerLib.embl > Libraries/RepeatMasker.lib

Then ./configure sees a non-empty .lib file and builds the database, and RepeatClassifier finally sees the resulting database files.
But why is this needed?
Shouldn't RepeatMasker come with a non-empty FASTA file if it's needed by its own configure script?  I can see bioconda's build.sh doesn't do anything in particular to handle that file, but, maybe it should?  This is probably where somebody who knows these tools better can chime in.)


Gold standard benchmark This page is claimed to contain a gold standard benchmark for viral genome assembly.
https://github.com/cbg-ethz/5-virus-mix
The claim is here:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5411778/
"This data set was recently presented as a gold standard benchmark (Di Giallonardo et al. 2014)"
But at (Di Giallonardo et al. 2014) I can only find that mention:
"We previously estimated the frequencies of the virus strains in the 5-virus-mix by amplifying the protease gene using single-genome amplification (SGA), the current gold standard for studying diversity of virus populations".
That doesn't seem to mean precisely that the dataset is a gold standard benchmark for genome assembly but, as I understand, that SGA is a gold standard for something else.
My question is: is it? And who decided it is? Is still valid? Is there an official list of results somewhere for this benchmark set?
I tried to contact the owner of the page but still no answer.

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


How to generate conformations for a ligand for in silico screening? I have multiple ligands (in pdb format) which I want to try and dock to various proteins.

But for each ligand I have only a single conformation (by which I meant the spatial relationship between it's atoms, not its position relative to the target protein).

How do I generate other possible conformations?
Is there a generally used tool for that?

One option I found was Frog2

https://mobyle.rpbs.univ-paris-diderot.fr/cgi-bin/portal.py#forms::Frog2

https://github.com/tuffery/Frog2

From an input of a mol2 or sdf file it creates a plethora of conformers, which can be outputted as a pdb file


How to get UCSC Genome Browser to show the actual sequence? Is it possible to make UCSC Genome Browser to show the sequence (nucleotide chain) of the current region?

I fiddled with the various viewing options and couldn't figure it out.

If you zoom in at one specific region, it shows the nucleotide chain on the third line (below the chromosome number). It can't be there by default as the size of default region is more than 100,000bp. You can use the scale on the top of the page to see the size of the region.


Adding entries to bigwig file I generate bigwig files using a shell script based on bedtools genomecov (to generate a bedgraph from a bam file) bedmap (to compute means across 10 bp bins) and bedGraphToBigWig to convert the binned bedgraph into bigwig.

Sometimes, the bam file has no data on some chromosomes, this results in a bedgraph file with no entries for the empty chromosomes.

I would like to add the missing entries. How could I do this?

I tried to see what I could do with pyBigWig, opening the bedgraph file in "a" mode, but don't really know how to proceed:

# I have the chromosome sizes for my genome of interest already in a list:
chrom_sizes
#[('I', 15072434),
# ('II', 15279421),
# ('III', 13783801),
# ('IV', 17493829),
# ('MtDNA', 13794),
# ('V', 20924180),
# ('X', 17718942)]
bw_test = pyBigWig.open("path/to/bigwig", "a")
bw_test.chroms()
# {'IV': 17493829, 'V': 20924180, 'X': 17718942}
bw_test.addHeader(chrom_sizes)
# RuntimeError: Received an error in bwCreateHdr

Presuming you still have your BAM file sitting around somewhere:

#!/usr/bin/env python
import pysam
import pyBigWig

bam = pysam.AlignmentFile("alignments.bam")
ifile = pyBigWig.open("your current file.bw")
ofile = pyBigWig.open("your new file.bw", "w")

# Create/Add the header
hdr = []
for chrom in bam.header['SQ']:
    hdr.append((chrom['SN'], chrom['LN']))
bam.close()
ofile.addHeader(hdr)

# Write the bigWig
for _ in hdr:
    try:
        ints = ifile.intervals(_[0])
        starts = [x[0] for x in ints]
        ends = [x[1] for x in ints]
        vals = [x[2] for x in ints]
        ofile.addEntries([_[0]] * len(starts), starts, ends=ends, values=vals)
    except:
        # intervals() will throw an error if the chromosome isn't present
        pass
ofile.close()
ifile.close()

The reason you can't append to a bigWig file is that it has a various indices and zoom levels typically built in. You'd have to add to them without overwriting neighboring blocks of data. In theory this might sometimes be theoretically possible, but not typically and I'm not aware of tools that support this.


What is the best way to start with Structural bioinformatics with zero prior experience? My undergraduate thesis is on the topic "Ligand design from protein-protein interfaces" and my PI has asked me to learn the basics of Structural bioinformatics over the summer so that I'll be able to start with my project as I begin my next semester(last week of July). Can someone please suggest me top resources to learn the basics of structural bioinformatics. Though I have done projects in Genomics, and Systems biology, this will be my first time with structural bioinformatics, so, I'll really appreciate some good leads. Thanks!

As others have said, there is no short cut and you will need to spend some time with a good book.

I would say there are 3 broad areas of understanding you have to have as a structural bioinformatician:

Solid understanding of the basic science. This includes molecular biology but also chemistry, which is often overlooked.
Knowledge of the databases and data formats available in the field. You need to know where to go to get protein structures, protein sequences etc.
Programming and algorithmic skills. You don't have to be a computer whiz but being effective in a high-level language (e.g. Python, R, Julia) and being able to parse data and carry out basic data analyses are essential.

If you want to take it further and have a bit more time then I would consider implementing something yourself such as a protein secondary structure predictor. Getting this working is challenging and would show you have competency in the above three areas.


What is the best way to account for GC-content shift while constructing nucleotide-based phylogenetic tree? Let's say I want to construct a phylogenetic tree based on orthologous nucleotide sequences; I do not want to use protein sequences to have a better resolution. These species have different GC-content.

If we use a straightforward approach like maximum likelihood with JC69 or any other classical nucleotide model, conserved protein coding sequences of distant species with similar GC-content will artificially cluster together. This will happen because GC-content will mainly affect wobbling codon positions, and they will look similar on the nucleotide level.

What are possible ways to overcome this? I considered the following options so far:

Using protein sequence. This is possible of course, but we lose a lot of information on the short distance. Not applicable to non-coding sequences.
Recoding. In this approach C and T can be combined into a single pyrimidine state Y (G and A could be also combined in some implementations). This sounds interesting, but, first, we also lose some information here. Mathematical properties of the resulting process are not clear. As a result, this approach is not widely used.
Excluding third codon position from the analysis. Losing some short-distance information again. Also, not all synonymous substitution are specific to the third codon positions, so we still expect to have some bias. Not applicable to non-coding sequence.

It should be possible in theory to have a model which allows shifts in GC-content. This will be a non time-reversible Markov process. As far as I understand there are some computational difficulties estimating likelihood for such models.

There are models that take into account compositional heterogeneity both under the maximum likelihood and Bayesian frameworks. Although the substitution process is not time-reversible, the computations are simplified by assuming that the instantaneous rate matrix can be decomposed into an "equilibrium frequency vector" (non-homogeneous) and a symmetric, constant exchange rate matrix.

I guess all your suggestions are also valid, and I remember recoding being used successfully to reduce the GC-content bias (examples in the references above and here).


Calculate average parental genotype I am interested in creating an averaged value for genotype across two parents:

E.g. Input (three unrelated couples, 1, 2 and 3):

Individual    SNP1    SNP2    SNP3    SNP4 ...
Father_1      1       0       1       2    ...
Mother_1      0       0       2       1    ...
Father_2      1       0       1       1    ...
Mother_2      2       1       0       1    ...
Father_3      1       2       0       1    ...
Mother_3      1       .       0       2    ...

E.g. Output (the midparental genotypes for couple 1, 2 and 3):

Individual    SNP1    SNP2    SNP3    SNP4 ...
Midparent_1   0.5     0       1.5     1.5  ...
Midparent_2   1.5     0.5     0.5     1    ...
Midparent_3   1       .       0       1.5  ...

I can think of a way to do this using e.g. a PLINK .raw file and manipulating that in R, but am concerned that would become memory- and storage-heavy when assessing 100000s of SNPs and 1000s of pairs.

Does anyone know of an existing tool that can do this?

Thanks

Beside the correct question by @Kamil, this sounds like a job for datamash.

First let's get the number of columns in the file

$ head -n1 genotypes.txt| awk -v FS="\t" '{print NF}'

Now lets calculate the average per column with datamash

$ datamash -H mean 2-<n> < genotypes.txt > mean.txt

Replace <n> with the number of columns we received with command above.

For your example the output will look like this:

mean(SNP1)  mean(SNP2)  mean(SNP3)  mean(SNP4)
0,5 0   1,5 1,5

One can remove the mean with sed

$ sed -i -e 's/mean(//g' -e 's/)//g' mean.txt


The state, limitations and comparisons of large variant stores Background: We're increasingly needing some way of storing lots of variant data associated with lots of subjects: think clinical trials and hospital patients, looking for disease-causing or relevant genes. A thousand subjects is where we'd start, there's talk of millions on the horizon. With various genomic medicine initiatives, this is likely a wider need.

The problem: While there's plenty of platforms out there, it's a rapidly evolving field. It's difficult to get a feel for how (and if) they perform and how they line up against each other:

What's scalable and can handle a lot of data? What sort of limits?
What's robust and not a teetering pile of hacked-together components?
What has a large community behind it and is actually used widely?
What makes for easy access and search from another service? (Commandline, REST or software APIs)
What sort of variants they handle?
What sort of parameters can be used in searching?

Solutions I've seen so far:

BigQ: used with i2b2, but its wider use is unclear
OpenCGA: looks the most developed, but I've heard complaints about the size of data it spits out
Using BigQuery over a Google Genomics db: doesn't seem to be a general solution
Gemini: recommended but is it really scalable and accessible from other services?
SciDb: a commercial general db
Quince
LOVD
Adam
Whatever platform DIVAS & RVD run on: which may not be freely available
Several graphical / graph genome solutions: We (and most other people) are probably not dealing with graph genome data at the moment, but is this a possible solution?
Roll your own: Frequently recommended but I'm sceptical this is a plausible solution for a large dataset.

Anyone with experience give a review or high-level guide to this platform space?

An epic question. Unfortunately, the short answer is: no, there are no widely used solutions.

For several thousand samples, BCF2, the binary representation of VCF, should work well. I don't see the need of new tools at this scale. For a larger sample size, ExAC people are using spark-based hail. It keeps all per-sample annotations (like GL, GQ and DP) in addition to genotypes. Hail is at least something heavily used in practice, although mostly by a few groups so far.

A simpler problem is to store genotypes only. This is sufficient to the majority of end users. There are better approaches to store and query genotypes. GQT, developed by the Gemini team, enables fast query of samples. It allows you to quickly pull samples under certain genotype configurations. As I remember, GQT is orders of magnitude faster than google genomics API to do PCA. Another tool is BGT. It produces a much smaller file and provides fast and convenient queries over sites. Its paper talks about ~32k whole-genome samples. I am in the camp who believe specialized binary formats like GQT and BGT are faster than solutions built on top of generic databases. I would encourage you to have a look if you only want to query genotypes.

Intel's GenomicDB approaches the problem in a different angle. It does not actually keep a "squared" multi-sample VCF internally. It instead keeps per-sample genotypes/annotations and generates merged VCF on the fly (this is my understanding, which could be wrong). I don't have first-hand experience with GenomicDB, but I think something in this line should be the ultimate solution in the era of 1M samples. I know GATK4 is using it at some step.

As to others in your list, Gemini might not scale that well, I guess. It is partly the reason why they work on GQT. Last time I checked, BigQuery did not query individual genotypes. It only queries over site statistics. Google genomics APIs access individual genotypes, but I doubt it can be performant. Adam is worth trying. I have not tried, though.


Are there free cloud computing platforms for biology projects? I want to implement the analysis found in the paper RNA-Seq of Tumor-Educated Platelets Enables Blood-Based Pan-Cancer, Multiclass, and Molecular Pathway Cancer Diagnostics The project id is 281708. It has 285 samples comprising about 300 Gb of read data. The reads are available on sra.

I want to apply a machine-learning pipeline to find important features in these data but I do not have access to computing resources to finish the project.

Question

Are there free cloud computing platforms that I can use to perform this research?

If you're getting started with a project and need some free computational resources you should look at https://galaxyproject.org/. There are some free servers that you may be able to use. However, Galaxy is designed to have pre-fabricated pipelines and a GUI to avoid having to deal with command-line tools. Therefore any sort of custom analysis will have to go through the galaxy API.

Otherwise, most of the free computational resources that were available several years ago like iPlantCollaborative are no longer available. Running servers costs money, after all, and the funding for such projects often runs out after a few years.

You may consider asking researchers at your local universities to collaborate and provide you with access to their servers.


Removing data from OrthoFinder output prefix I am using OrthoFinder to find orthogroups between species. I am using the latest version unpacked through tar as per their instructions.
Orthofinder creates an output folder according to the current date, ie, Results_Dec25. However,
I want to integrate orthofinder in snakemake, so having a stable output directory is necessary.
The -o option allows specification of an output prefix, but the dated folder is then created within.
Is it possible to force a given name to the output folder?
Or is there a way in snakemake to allow ambiguous output directory names?
I know this would be possible via editing the source code of OrthoFinder, but this workflow needs to be reproducible without faff.
Many thanks
This question has already been asked here but remains unanswered.
Is there an option in Orthofinder to choose the output directory name?

The latest version solves the issue.
Source : https://github.com/davidemms/OrthoFinder/issues/403
Latest Version: https://github.com/davidemms/OrthoFinder/releases/tag/2.4.0


Very low probability densities I'm generating a density plot using the ggplot2 library:

ggplot(df_chr20, aes(x=start1)) + geom_density() + theme_bw()
ggplot(df_chr21, aes(x=start1)) + geom_density() + theme_bw()

The data frame for df_chr20 has 1573 rows, and I wanted to plot the probability density of genomic regions as x = start1. The same for df_chr21, with 1271 rows. What I see in the plot is what I expect, but is there a reason why the y-axis probablilty densities are so low? The area under the curve is surely not equal to 1.

I expected the high localised high density in chr 21, and both the plots have similar scales in terms of y-axis values. But since they are so low, how can I interpret this plot?

I found some information here, but my values are much lower than the examples.

Any help would be very kind!

It is not like your conventional frequency which adds up to one. The density is low because the width of your bins is huge, and the number of observations you have is low. From this by Wickham, the basic kernel is:

where K is the kernel and h is the bandwidth (Scott, 1992b)

If you don't specify it, by default h will revert to nrd0:

library(rtracklayer)
gr_narrowPeak <- import("https://www.encodeproject.org/files/ENCFF072FBH/@@download/ENCFF072FBH.bed.gz",format = "BED",extraCols = extraCols_narrowPeak)
gr_chr21 = gr_narrowPeak[seqnames(gr_narrowPeak)=="chr21"]

da = as.data.frame(gr_chr21[rep(1:length(gr_chr21),gr_chr21$signalValue)])
dens = density(da$start)
dens$bw
[1] 1133889

plot(dens)

ggplot(da,aes(x=start))+geom_density()

You can decrease h such that you get something smaller (but at the cost of looking coarse). Important thing is it doesn't sum up to one.

If you would like to explain it easier like a histogram, consider converting it into counts, like explained in this post:

binW = 50000

ggplot(da,aes(x=start))+
geom_density(aes(y = ..density..*(nrow(da)*binW)),bw=50000) +
geom_histogram(binwidth=binW,fill="steelblue",alpha=0.5) +
theme_bw()

ggplot(da,aes(x=start))+
geom_histogram(aes(y=..count../nrow(da)),binwidth=50000) +
geom_density(aes(y = ..density..*binW),bw=binW,fill="steelblue",alpha=0.5) +
theme_bw()


Oxford Nanopore mapping Quality and sequencing error I'm analyzing some Oxford Nanopore sequencing reads, with several articles stating the sequencing error is between 5% and 15%. Since its 1D and not 1D^2 sequencing i would assume the sequencing error is closer to 15%.

I'm using Minimap2 aligner to align the reads to the reference genome, and then im choosing those reads with a mapQ of 60, to analyze.

Where i know with the phred score mapQ of 60 means 99.9999% of the bases are called correctly.

So my question is, if the mapQ is 60, then does it really matter if my sequencing error is 15%? Or how is the sequencing error correlated with mapQ? Becuase i guess you could still obtain a good mapping quality with a few mismatches?

Or am i completely off target.

Thanks for your time and help.

Where i know with the phred score mapQ of 60 means 99.9999% of the bases are called correctly.

That's not at all what the MAPQ means, it means that the sequence is very likely to have originated from that place in the genome, whether the underlying sequence is correct is another question all together.

For long nanopore reads, sequence quality is not well correlated with MAPQ. MAPQ will tend to be either very low (0 or 1) or very high (60) since you tend to have kb+ reads and those will largely be uniquely assignable in the genome. The sequencing quality of nanopore reads is only an issue when you care about something at the base level, such as SNP calling or finding exact splice junctions. Even then, you mostly need sufficient coverage to get around this.


making a customer reference by using cell ranger Anyone here knows how to add a new gene (GFP) to customer reference by using cellranger? We use the Cre/loxp system by using R26RmTmG mice and want to use GFP to isolate specific cells. I have read the cellranger workflow but still, have no idea how to do it. Like it is said that I should update the GTF file as shown in the screenshot.

Step 1:
Add the the GFP sequence to your fasta file

>GFP_name
ATG.....TAA

Step 2:
Annotate the gtf so that it contains all the necessary information to annotate your information

Add this line to your GTF file but modify where necessary:
"Make sure that all columns are tab separated"

GFP_name    AddedGene   exon    1   1966    .   +   1   gene_id "GFP"; transcript_id "GFP"; gene_name "GFP"; transcript_name "GFP"

The first column has to correspond to the header of the fasta sequence
start and end indicates the length of your sequence. I added a 1966 bp long sequence to my gtf

Code example:
step 1:

cat reference.fa gfp_seq.fa > new_reference.fa

step 2:

cat genes.gtf gfp.gtf > new_genes.gtf

-  gfp_seq.fa is your custom gene in fasta file
-  gfp.gtf is your custom gtf line or lines if you want to add more

You can also open your current reference fasta and gtf file in a text editor and copy paste the sequence and gtf line in there.

Then run cellranger mkref on your new_reference and new_genes file


Trouble using biomaRt to retrieve hgnc symbols from Ensembl transcript ids I have a matrix of gene counts which I'm going to use as input for DESeq.  Right now, each gene is labeled by its Ensemble transcript ID, but I'd like to convert these to their HGNC symbols before I input them into DESeq for analysis. I'm attempting to do this conversion using biomaRt, however, this is my first time using the program, and I'm running into issues.  Here is my current code:

library("biomaRt")
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "mmusculus_gene_ensembl")
transcript.ids <- rownames(txi.kallisto$counts)
hgnc_symbols <- getBM(attributes = "hgnc_symbol", filters = "ensembl_transcript_id", values = transcript.ids, mart = mart)

Printing hgnc_symbols shows that it isn't what I want:

>ghnc_symbols
[1] hgnc_symbol
<0 rows> (or 0-length row.names)

Lastly, here is what the input transcript IDs look like:

> head(rownames(txi.kallisto$counts))
[1] "ENSMUST00000178862.1" "ENSMUST00000178537.1" "ENSMUST00000196221.1"
[4] "ENSMUST00000179664.1" "ENSMUST00000177564.1" "ENSMUST00000179520.1"

You need to specify the number without the version. Instead of "ENSMUST00000178862.1" just "ENSMUST00000178862":

You can do this with one more line:

g <- gsub("\\..*", "", rownames(txi.kallisto$counts))
(hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = g, mart = mart))
##  hgnc_symbol chromosome_name ensembl_transcript_id
##1          NA              14    ENSMUST00000177564
##2          NA               6    ENSMUST00000178537
##3          NA               6    ENSMUST00000178862
##4          NA              12    ENSMUST00000179520
##5          NA              14    ENSMUST00000179664
##6          NA              14    ENSMUST00000196221

With these transcripts it doesn't seem that there is an equivalent hgnc_symbol. But is recognized, otherwise they are not recognized:

(hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = rownames(txi.kallisto$counts), mart = mart))
## [1] hgnc_symbol           chromosome_name       ensembl_transcript_id
##<0 rows> (or 0-length row.names)

However you might be interested in doing something like this:

human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
getLDS(attributes = "ensembl_transcript_id",
       filters = "ensembl_transcript_id", values = g,mart = mouse,
      attributesL = "hgnc_symbol", martL = human)
## The query to the BioMart webservice returned an invalid result: the number of
## columns in the result table does not equal the number of
## attributes in the query. Please report this to the mailing list.


How to find/build the evolutionary history of a protein from its sequence? I'd like to build the evolutionary history of a protein, given its sequence. Namely, given a FASTA entry how can I build an evolutionary tree? Here is the 5wxy protein as an example:
>5WXY:A|PDBID|CHAIN|SEQUENCE
MGHHHHHHMMKTKLPILGVLGGMGPVVTAEFLKSIYEYNPFIDKEQESPNVIVFSFPSAPDRTGSIDSGKEREFIDFIQV
NLEHLNKLADCIVIGSCTAHYALPQIPENLKDKLISLIKIADQELQEYNEPTLLLASTGTYQKKLFQEGCTTADLIISLS
ESDQKLIHEMIYKVLKRGHDPLSILRDIEALLEKYNTRSYISGSTEFHLLTKSLKLKGIDSIKAIDPLSTIAQNFSQLII
KQAQVDLVTDCHQPSNPKSP

Do I need to do some queries to Entrez in order to build the evolution tree?
I saw this paper, but I am not sure how I can use it with a single mammal's protein primary sequence.

Question given a protein name, is it possible to look for its evolutionary history or build evolution tree using other databases beside PDB on top of bioinformatics algorithms.

The general procedure is:

Find the sequence for the same (or as similar as you can find (see mentions of blast in the comments above)) gene in other species. You can use uniprot or any other large sequence database for this.
You now have a large set of sequences related to yours, which you'll need to compare. Search the methods section of your favorite paper to see how they reconstructed likely evolutionary history from this, there are a number of packages out there.


Compare two networks of the same genes between two species I had a set of genes for honey bees. I found the homologues of those genes in the house fly and created a network of the honey bee genes and a network of the homologous fly genes using STRING.

I know that I could do a topological comparison between the networks and check if the same set of genes give the same hubs when in two different species. However, I can't find a single study where something similar has been done.

What biological inference can I draw from the topological parameters? For example, if the diameter of the fly network is smaller than the honey bee network, what inference could I draw?

UPDATE: It would be great if atleast some publications could be given where network comparisons have been done between different species.

So, you can give a try to this Cytoscape app that will allow to perform comparison between networks.

This app is not specifically designed for comparing networks across different species but it will enable you comparing any kind of networks.

These are the papers they published, which can give you some hints.

About your specific problem...I can't help you about what to look for when you're comparing networks from different species (not an expert in this specific topic) but as far I know you should look for something related to patterns occurring across different networks (i.e. species), network motifs, distribution of network centralities, conserved structures. In general, to compare networks, I would look for global properties (e.g. diameter, clustering coefficient...) while I would use local properties (e.g. degree, closeness, and network centralities in general) to compare the role of a node that you retrieved in all your networks.

Also, since you're looking for some literature, these links may help:

About conserved patterns of protein interaction in multiple species
About comparing genomic expression patterns across species
This one is slightly out of topic but can give you
some other ideas about comparing networks.


size of the pathways for analysis and filtration I have recently started working on a substance's effect on a cell line in different dosages. for this, there is a tool called bmdexpress2 that I am using. Its input is the normalized counts from RNASeq for each dosage as a big matrix. When it comes to the pathway analysis step, unlike hallmarks of GSEA, this tool uses some databases for defining pathways which involve pathways of even 2-3 genes.

So what I want to discuss here is the strength of small pathways and filtering thresholds. How can we bioinformatically decide on the importance of a pathway made of 3 genes? Should we just filter them out as there are much more comprehensive pathways? Or do their GO Levels matter? Also there are many cases that out of those 3 genes, 1 of them is differentially expressed; is having 1 gene diff. expressed out of X genes in a pathway enough to say that the pathway is enriched?

Here is an example GO term for a small pathway.

And here is an example of the output I get from the tool. all genes platform column is the number of genes involved in that pathway, input genes column is how many of the genes in that pathway are significantly differentiated. I believe the rest is self explanatory.

Looking at the table, can we really say, ESCRT III Complex, for example is a good pathway to focus while having 1 gene out of 13 in a pathway but that change is significant?

Thanks

The problem here is that any one GO term with a small number of members is unlikely to include any differentially expressed genes by chance, but if you test a large number of GO terms, some of them definitley will. The solution is to correct the p-values for multiple comparisons. If you do this, then none of those GO terms with a single DE member will be significant. The larger GO terms with multiple DE members in the chart you show probably won't survive correction either. Hopefully you have some that are more significant.


What does instantaneous rate matrix mean? PhyML is a tool used for maximum likelihood estimation analysis.
I was running a PhyML analysis to check the rate of nucleotide change in DNA sequences and in the PhyML output, there are rate parameters and an instantaneous rate matrix. Could you please clarify the difference between the two?

What I've understood so far is that the rate parameters indicate, for ex: that the A -C transition rate is 1.95 and so on. If so what does the instantaneous rate matrix tell us?

Reversible At its simplist the GTR rate, is the General Time Reversible model and most importantly infers the matrix is "reversible", thus as many A mutations will go to C as C mutations go to A and the number of mutations is relative to the G-T mutation rate. The G-T mutation rate is one of the least common mutations (transversion) so all other values are above 1. If you plotted this as a matrix the top right-hand matrix is identical to the bottom left hand matrix because its reversible, i.e. A <-> C . This is cool because you don't need to tell the algorithm which sequence is the root.
Directional The "instantaneous" matrix is a directional matrix, thus A -> C is not the same as C -> A. In order to obtain an accurate directional matrix you need to assign a root, its like lets start from the deepest node in the tree (the one between the root/outgroup and all ingroup sequences) and now - if that is correct - we can remove the restrictions of reversibility and impose directionality on the mutation process. This is needed if you are performing molecular clock type calculations. However, the key is assigning the root/outgroup - if it assigns it e.g. by mid-point rooting that is not necessarily a good idea.
Advice. My advice is if you are building a basic tree and haven't assigned a root/outgroup, obviously using GTR which is a good move, then you don't to worry about the intananeous matrix because it might not be meaningful. If you are running a molecular clock/dating/rate estimate type model and a very clear what the root/outgroup is and that is specified in the input, then the instananeous matrix becomes important.

The mutation rates of GTR and instananeous are not directly comparable because the GTR rate is with respect to G <-> T mutation, whilst instananeous is primarily concerned with the directionality of mutation and is a 12-parameter model, as opposed to GTR which 6-perameters, where one of those parameters is the reference.
Thus G-T is reversible in GTR (the reference) and so e.g. A<->C is reversible with respect to G<->T.
In instananeous matrix G -> T and T -> G are two different rates as are A-> C and C -> A, so one approximation based on an assumption with respect to a second approximation based on the same assumption can't directly equate to a mutation rate which removes that assumption.


Use Htslib to extract specific region of vcf file I am working in c++ using the htslib library to parse different vcf files as well as fasta files.
My goal is to read in a specific region from a fastafile and then extract all the VCF sites within that region from a seperate vcf file, in this case only containing chr14. I can extract the fasta region in several ways, one of which is to use the function faidx_fetch_seq() and define a region of chr14:20000000-22000000.
Then for the bcf files i am able to iterate through the entire file and extracting the different information regarding the reference, alternative allele and so on. With a simple example shown below.
void faVCF(const char* fastafile,FILE *fp){
  htsFile *test_bcf = NULL;
  bcf_hdr_t *test_header = NULL;

  test_bcf = bcf_open("chr14.vcf", "r");
  test_header = bcf_hdr_read(test_bcf);
  bcf_hdr_set_samples(test_header, "20000000-22000000", 0);
  bcf1_t *test_record = bcf_init();
  while(bcf_read(test_bcf, test_header, test_record) == 0) {
              std::cout << bcf_hdr_id2name(test_header, test_record->rid) << "\t" <<
                          test_record->pos << "\t" <<
                          test_record->n_allele << "\t" <<
                          std::endl;
          }
}

However, i doesn't work to extract all the variant sites in a specific region. I tried bcf_hdr_set_samples with the argument "20000000-22000000", to extract that given region from my chromosome. As far as I understand this should be able to subsample both samples and region queries using the header information, and then upon initiation extract that region, however this doesnt work.
So my question is: Can anyone help me with what function I can use to extract all the vcf sites from the region chr14:20000000-22000000.?

For iterating over reads in a region in a single BCF file, you want to use the iterator interface: bcf_itr_queryi() or bcf_itr_querys(), together with bcf_itr_next(). For a VCF file, it appears you need to use the similar TBX iterator interface: tbx_itr_queryi() or tbx_itr_querys(), together with tbx_itr_next(). See HTSlib's tabix.c for reasonably straightforward examples of both.
Alternatively you could use the routines in htslib/synced_bcf_reader.h, which encapsulates both and is what bcftools mostly uses.


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

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

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

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

import pysam

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

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

This outputs lists of tuples of the parse output:

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

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

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


How do I get the adjacency matrix from gene similarity network? I am exploring application of Graph convolution network on gene similarity network. I am trying to explore
genemania.org and https://string-db.org/

Please let me know where can I find the data for my research in gene interaction network, which allows me to download the data in the form of an adjacency matrix.

The STRING download page is here. I am not sure that they provide any data in the form of an adjacency matrix, as the number of nodes (rows/columns) would be quite large, and for this reason they are likely to represent the data as lists of links, for compactness.

Their data will likely have to be cleaned for any purpose you are interested in. You will have to make some decisions on how to filter the data. For example, as suggested in the comments, you may want to filter on the quality of the evidence supporting interactions, or based on species of the protein, or other things.

It's possible that you may have better luck going with another database that has done some of these steps already. For example, there are relatively compact human interactome datasets here, or yeast interactome datasets here. In this case you will still need to convert from a list of links to an adjacency matrix, however this is not very difficult if you can use a scripting language.

Notably, you may have to figure out how to distinguish missing data from tested non-interactions in these cases (or decide that missing data is ok). In most cases people will not have exhaustively tested every cell of the all-by-all interactome adjacency matrix.

Hope that helps.


How can atomic RMS differences / atomic RMS distributions help us understand the structure of a protein? For example, I have a paper that says: "the atomic RMS distribution about the mean coordinate positions is 0.41 +/- 0.04 A for the backbone atoms and 0.82 +/- 0.04 A for all atoms"
And I have no idea what any of this means, I found several papers with this kind of information but none of them explain why this is relevant in order to understand the structure of a protein.

The "RMS distribution about the mean coordinate positions" could be restated as "certainty of the structure". In other words, the main backbone that defines the large-scale shape in the structure is expected to be accurate to within about 0.41 angstrom, while the atoms more on the periphery (or otherwise not forming the backbone of the structure) are less accurately solved. How important this is will depend on what you want to do with the structure. If you want to get an idea for what it looks like then the accuracy is less important, if you want to do docking then the more accurate the atoms around any sort of pocket/interaction surface are the better.


How to do BLASTP for short sequences (<20 aa) effectively? I am trying to find sequence homology between viral sequences and my protein of interest. I have the sequences of their epitopes which varies from 5 to 500 amino acids long. For shorter sequences, it is recommended to use blastp-short but I am confused until which length it is okay to run blastp?

For shorter sequences, what are the parameters needs to be added for effective results?

blastp -task blastp-short -query POI.fasta -db virus.fasta -out output.txt -outfmt 6 -evalue 20000

You should increase the -evalue 10 default option to a higher value because this value it's adjusted to the sequence database size. Check this table, maybe this can help you to decide which parameters are the best option for you.

    OPTION      TASK          DEFAUL    DESCRIPTION
    word_size   blastp        3         Word size of initial match.
    word_size   blastp-short  2         Word size of initial match.
    gapopen     blastp        1         Cost to open a gap.
    gapopen     blastp-short  9         Cost to open a gap.
    gapextend   blastp        1         Cost to extend a gap.
    gapextend   blastp-short  1         Cost to extend a gap.
    matrix      blastp        BLOSUM62  Scoring matrix name.
    matrix      blastp-short  PAM30     Scoring matrix name.
    threshold   blastp        11        Minimum score to add a word to the BLAST lookup table.
    threshold   blastp-short  16        Minimum score to add a word to the BLAST lookup table.
    window_size blastp        40        Multiple hits window size, use 0 to specify 1-hit algorithm.
    window_size blastp-short  15        Multiple hits window size, use 0 to specify 1-hit algorithm.


How do you visualize minimap2's PAF output format? As described by the developer of minimap2:

PAF is a text format describing the approximate mapping positions between two set of sequences.

He also suggest to use either Ribbon to visualize PAF files, or to use minidot from miniasm:

minimap2 -c -x??? asm1.fa asm2.fa > out.paf
minidot out.paf > out.eps && epstopdf out.eps

What does -x??? stand for? And how do you visualize PAF files? Could you think of a way to use IGV to visualize PAF files?

As described in the minimap2 manual, -x STR is used to specify a preset option like -x map-pb for PacBio/Oxford Nanopore read to reference mapping (-Hk19).
In general, I use the -a flag to output alignments in sam format that can then be  sorted and converted to bam format, and viewed in any genome browser. The support for PAF format when it comes to visualization appears to be limited. A few options that I have found, but not tested are:

minidot
pafr
pafpy -- this one may not have the plotting capabilities though


Accuracy of the original human DNA datasets sequenced by Human Genome Project? The Human Genome Project was the project of 'determining the sequence of nucleotide base pairs that make up human DNA, and of identifying and mapping all of the genes of the human genome'. It was declared complete in 2003, i.e. 99% of the euchromatic human genome completed with 99.99% accuracy.

Are the datasets provided by HGP still accurate, or as accurate as was claimed in 2003?

Given the technology in the past (such as using old techniques), or any other reason (newer research studies), is it possible that the datasets are not as accurate as originally expected?

The HGP developed the first "reference" human genome - a genome that other genomes could be compared to, and was actually a composite of multiple human genome sequences.

The standard human reference genome is actually continually updated with major and minor revisions, a bit like software. The latest major version is called GRCh38, was released in 2013, and has since had a number of minor updates.

  Are the datasets provided by HGP still accurate?

Yes, in a sense, but we certainly have better information now. One way to measure the quality of the assembly is that the initial release from the HGP had hundreds of thousands of gaps - sequences that could not be resolved (this often occurs because of repetitive sequences). The newest reference genome has less than 500 gaps.


processing multiple fastq files with cutadapt I have DNA sample from 5 pools, having 25 fastq files each. I am running cutadapt to remove the primers using this command
cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG A7_S7_R1_001.fastq -o sb_A7_S7_R1_001.fastq --discard-untrimmed
However, this is too time consuming.  I tried my hand at writing shell script

for file in /dir/*
do
  cmd cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG "$file" >> results.out
done

But it is not running .
Can someone help.
Thanks

Try this one

for file in /dir/*
do
  cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG "$file" >> results.out
done

Also you have omitted --discard-untrimmed from the original command. If you want that as well along with everything to run in parallel:

for file in /dir/*
do
  cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG --discard-untrimmed "$file" -o "`basename -s .fastq $file`_out.fastq" &
done


How do I validate a single sample ArrayCGH result? We have arrayCGH (aCGH) results for one sample. There is a 0.5 Mb terminal duplication on chromosome 19 (62995490-63407936, according to NCBI36/hg18). The duplication is rare: a literature review suggests there are only 3-4 samples with clinical information.

What are the steps to validate the results? How do we ascertain that this duplication is the cause of the clinical symptoms?

I have some ideas:

aCGH the parents. Not sure how this would help.
whole genome exome sequencing. Worried this might make it more difficult to pinpoint genetic cause.
whole genome sequencing?
other ideas?

Note: I am new to aCGH and high-throughput sequencing, any advice is welcome.

Why not just good old qPCR? That's (A) quick, (B) cheap and (C) easy to analyze. If you care about the exact location of the break point (I'm guessing from the context that you don't), then targeted sequencing with a custom capture kit would work.

Regarding validating the biological relevance of this, there are multiple (parallel) routes one can take. Firstly, screen unaffected family members for this same alteration. If you find this in unaffected individuals then it's obviously not the causative alteration. Ideally, one would also do either a cell-line experiment or a mouse (or other model organism) experiment to see if you can reconstitute at least some component of the clinical phenotype. This may not always be possible, of course.


What datasets are available out there for prediction based on DNA sequences? I am looking for publicly available data for a genomics deep learning project. My goal is to compare different architectures to predict biological insights from DNA sequences.
I have heard about Janggu before, but I am looking for a few datasets which are like the MNIST equivalents in genomics. Kaggle doesn't seem to have many such datasets, and I would be happy if you share your insights about potential datasets that one could use.
For example, it could be datasets that map 100-Mers to their ChIP-seq values or datasets that annotate a DNA sequence with the probability that a transcription factor would bind to it.
Choong and Lee lists a few datasets on their IEEE paper from 2017:

Good question and the Janggu link is cool. Thanks!
Kaggle is/was heavily image recognition. There isn't an MNIST type data set in genomics, in case anyone is wondering this is a classic natural language processing (NLP) type analysis to identify handwritten single digit numbers.
There has been several FDA precision challenges, one being described here:
FDA precision challenges
This could provide something towards a 'standard'.
A key issue remains the accurate vectorisation of sequence and whether you think that is/was satisfactory in my personal opinion.
There is an old review article here


Get gene exons from UCSC browser by HGNC name I have a set of gene names in HGNC standard and I want to get all their exons (in BED format) from the UCSC browser but it seems like the browser doesn't index genes by their HGNC name. How can I do that?

You could grab exons via Gencode v27 (hg38) GFF annotations and convert them to BED via BEDOPS convert2bed/gff2bed:

$ wget -qO- ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/gencode.v27.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "exon"' \
    | convert2bed -i gff - \
    > exons.bed

Then use grep to filter exons against your list of gene names:

$ grep -w -f hgnc.txt exons.bed > exons.filtered.bed


Lisp in bioinformatics After many years of programming in Python, C/C++ and other mainstream languages, I am looking for a new language of pleasure, mainly as a part of mental hygiene. Since I really admire the ideas behind Lisp, including the original paper from the 1950's, I would like to learn some Lisp dialect.

Which Lisp dialects are currently suitable for bioinformatics (e.g., for developing simple command-line tools)? I have the following requirements:

Portability across main Unix platforms and reproducibility (this seems to be a huge issue in many Lisp-like languages)
Resulting programs being installable using Bioconda (it should be technically possible at least in future if not now)
Negligible starting time (this disqualifies Clojure, which requires approximately 10 seconds)

I like the idea but I’m not aware of a strong Lisp presence anywhere in bioinformatics. I think a good starting point would be the paper How the strengths of Lisp-family languages facilitate building complex and flexible bioinformatics applications by Khomtchouk &al (and the Hacker News discussion).

Excluding Clojure from consideration will unfortunately close off a great many modern resources (although I fundamentally agree with your technical reason for doing so). Apart from this, Scheme (in particular R7RS) is probably the most-used Lisp dialect, ahead of Common Lisp. Or at least the most actively developed one. As Pierre mentioned, there’s also GNU Guile but, to be honest, despite it being the default for GNU software tools I don’t think it’s frequently used outside GDB.

That said, Clojure’s startup time is not anywhere near 10 seconds on my system. It’s not exactly instantaneous but not much worse than other dynamic languages, despite being backed by the rather sluggish JVM.


What is a signal 'chunk' in the context of Nanopore sequencing? Does one chunk correspond to the raw signal of a single read read? Why does messing with this parameter affect basecalling results?

If you have access to the nanopore community website you can read more about this here: https://community.nanoporetech.com/posts/not-chunky-enough (I'm not sure to which extent I can copy paste from there so I won't).

A TLDR version: basecalling doesn't happen simultaneously for an entire read, but it is broken in segments which are processed separately and combined afterward. Those segments are chunks. Basecalling results are affected mostly by low-complexity regions, in which a heavily biased region (in terms of nucleotide content) leads to artifacts in the per-chunk basecalling. Those issues are solved by taking a chunking-size larger than most low-complexity regions.


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

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


Given a .gb file and a locus - how to get relevant annotations in Python? Given a .gb file and a specific locus in the genome - how can I retrieve the relevant annotations in Python (i.e., annotations that include that locus)?

I could retrieve the features using:

SeqIO.read('my_gb_file.gb', 'gb').features

and then scan them to find the relevant ones, but it feels like reinventing the wheel.

Is there a function in Biopython that does that?
Or in any other well-maintained package?

Biopython is the main package for this. It's only a few lines, so it is not reinventing the wheel. Unfortunately, this feels like a homework question, so nobody can write the code for you as per no homework policy.

But pointers are okay...

So you want to iterate across the features table (a list in Biopython) of a record and find the case where the qualifier['locus'][0] matches your query.

The things to watch out for are:

filter by type (CDS?)
entries where there is not key among the qualifiers, so add error catching (try: ... except Exception: pass
the values of the qualifiers of a feature are a list. so add a [0]
multithreading.Pool might help if speed is a worry (which it really shouldn't). Asyncio and thread are not the way (both single core).


Defining methylation level of CpG island Differentially methylated site is technically measured by each probe. However, CpG island within gene promoter is comprised of numerous probes. Do we have a consensus of how to determine the methylation level of CpG island using methylation level of probes within corresponding CpG island?

In my experience, the most common approach to find DMR is to use the DSS R package to call and merge DML into regions.
DMC/DML (Differential Methylation Cytosine/Loci) is often calculated by comparing the counts of methylated and unmethylated cytosine at one position between samples.
Steps are to first call DMLs with the callDML() function. Then call DMRs from these first results with callDMR(). Finally annotate the DMLs or DMRs using a GFF file. From there you can look for DMR falling into CpG island. Some of the parameters that can be tuned are the minimum region length, minimum number of differential CpG sites, percentage of CpG sites required to be significant. However I have never heard of what would be the standard values for these parameters. It may not be clearly established yet.
So because the process of calling DMLs and DMRs need robust statistical testing, it is suggested to first call the differential methylated states and then annotate it.


Is there a standard way to clean a PDB file and re-number its residues? Is there a pre-existing tool which will tidy up the numbering of a PDB file?

Firstly, I would like to re-number the residues on inserts to make the icode an actual residue in the chain (by that I mean it's own number, shifting everything after it, so 45A becomes 46).

Secondly, I would also like to re-number so that each residue is one after another, as there are skips in the numbering, in spite of the two residues being physically next to each other.

I could do this in-house with a python script, but I thought I would ask specifically if there is a library or pre-existing tool?

Worst case scenario I have found this Perl script that does similar.

Thanks.

Several years ago we have written a tool for this: http://kodomo.fbb.msu.ru/~avkitex/projects/pdbparser.html
https://github.com/avkitex/pdbparserTools

The description is only in Russian but you can contact the author.


select unique names while removing the duplicates from a column I have a file with Ensembl gene and Transcript IDs :

I want to retain only the first entires from both the column so that the output looks like this:

ENSG00000003137 ENST00000001146
ENSG00000003402 ENST00000309955

and ignore the rest.
Have tried awk awk '!a[$1$2]++' and versions of it. Nothing works.
Kindly suggest some way to get the desired output.

regards.

That sound like a job for a perl one-liner:

perl -ane 'unless(exists $seen{$F[0]}){print $_; $seen{$F[0]}++}'

Explanation:
The input is read line-wise and printed back to the command line if the ID in the first column has not been seen before.

perl -ane ' reads input line-wise and splits fields to array @F
unless(exists $seen{$F[0]}){ checks if ID exists in hash %seen and continues if that is not the case
print $_; prints the original input to STDOUT
$seen{$F[0]}++}' adds the ID to the hash %seen