How to split multifasta based on partial fasta header I have download a batch of refseq fasta files and want to split them based on strain. this is complicated by accession numbers, project IDs etc. I downloaded the genomes via eDirect from NCBI. I know I can split the output but I want them together for some preliminary parsing (remove ambiguous and plasmid sequences etc.)

See example below.

>NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
>NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
>NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
>NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC

I want two files:
Clostridium botulinum BKT028387.fas and Peptoanaerobacter stomatis OBRC8.fas

My actual dataset is much larger (~70,000 genomes) and not too sure about the headers in each case.

Im thinking the following could work but Im not sure how to execute it:

Using a grep command to extract the lines but split the output.

If I download the summary list from refseq and use the names as the query example:

grep -Fwf refseq.list sequences.fa

  Summary file

  Peptoanaerobacter stomatis strain OBRC8, whole genome shotgun sequencing project 7,810,226 rc other DNA  NZ_AJXR00000000.1 GI:938474759

  Clostridium botulinum BKT028387, whole genome shotgun sequencing project 4,134,593 rc other DNA  NZ_NOJY00000000.1 GI:1231817869

  But I could edit it down to

  Peptoanaerobacter stomatis strain OBRC8

  Clostridium botulinum BKT028387

Am I one the right track. How do I finish this?

The SeqBuddy --pull_records command can grab sequences based on a regular expression against record IDs and/or metadata. Given your example:

$: seqbuddy.py sequences.fa --pull_records "Clostridium botulinum.*BKT028387" full > Clostridium_botulinum_BKT028387.fas
$: seqbuddy.py sequences.fa --pull_records "Peptoanaerobacter stomatis.*OBRC8" full > Peptoanaerobacter_stomatis_OBRC8.fas


Accessing .bam/.cram files from AWS S3? What's the best/easiest way for accessing .bam/.cram files from S3?

I have .bam/.bai .cram/.crai and .bed + .gff files that I want to make available easily to others so they can browse in IGV and other solutions from their computers. There seem to be different solutions after cursory Googling, what would you recommend?

Here is what I've found so far:

IGV:

https://twitter.com/igvteam/status/1187766079457644544
https://umccr.org/blog/igv-amazon-backend-setup/

Samtools: samtools directly supports AWS

http://www.htslib.org/doc/htslib-s3-plugin.html

The full URL format is:

    s3[+SCHEME]://[ID[:SECRET[:TOKEN]]@]BUCKET/PATH

You could set up a public read-only bucket in an S3 account and put your files there. URLs would look like:

https://${username}.s3.amazonaws.com/public/${filename}

For samtools, the request for a BAM file's index should be built off that URL. For IGV, I think you can specify separate, custom URLs for the data and the index resources.

Public buckets are available to all, without need for credentials. So that might be an issue to think about, with respect to billing for web requests and data transfer.

I'm still not 100% certain that samtools/tabix support https requests reliably (I have a prior question on here about that), but that might be one approach.


How to get Uniprot and swissprot cross references using Biopython I am struggling to use biopython to gather cross references such as GO annotation from Uniprot/swissprot text files ("DR" in the .txt files) using biopython.

From the source code:

def _read_dr(record, value):
    cols = value.rstrip(".").split('; ')
    record.cross_references.append(tuple(cols))

This leads me to think that the following code should work, as this is how other features like annotation and features work.

for record in SeqIO.parse(filename, "swiss"):
    print(record.cross_references) #does not work, should return a list of DR rows
    print(record.features) #works, returns a list of parsed FT rows

I believe the attribute you should be looking for is dbxrefs:

[x for x in record.dbxrefs if x.startswith('GO:')]

Apparently SeqIO handles uniprot text files differently than the SwissProt module:

from Bio import SwissProt
record = SwissProt.read(open('single_record_from_uniprot.txt'))
record.cross_references

So you should change you code to use the SwissProt module to get the behavior you expect


variant calling on ChIP-seq style data: samtools mpileup with minimal filters I am running samtools mpileup (v1.4) on a bam file with very choppy coverage (ChIP-seq style data). I want to get a first-pass list of positions with SNVs and their frequency as reported by the read counts, but no matter what I do, I keep getting all SNVs filtered out as not passing QC.

What's the magic parameter set for an initial list of SNVs and frequencies?

EDIT: this is a question I posted on "the other" website, but didn't get a reply there.

I used this in the past for ChIP-seq data and it generated SNVs:

samtools mpileup \
--uncompressed --max-depth 10000 --min-MQ 20 --ignore-RG --skip-indels \
--fasta-ref ref.fa file.bam \
| bcftools call --consensus-caller \
> out.vcf

This was samtools 1.3 in case that makes a difference.


How do you find differences when your metagenome sample size is small? I received a matrix of genomics features (KEGG annotated metagenomes) abundance from 6 samples belonging to 2 groups. My collaborator is interested in finding differences between the two groups.
The matrix, after scaling and removing zero-variance, is roughly 4000 features. Via shapiro test I already know some 30% of the features are not normally distributed across the samples, so I assume PCA is out of option and any test with the normal distribution assumption.
I already made some preliminary analysis using NMDS to cluster my samples and was thinking of trying also t-SNE and recursive feature elimination. However I am not sure how to find "significant" differences given the small sample size. Any suggestion? Would something like LEfSe work in this instance?

I think there are several ways to figure this out. I wouldn't worry too much about normality as most of the methods are either nonparametric or pretty robust to violations.

Since you are using KEGG, I would suggest working at the pathway or module level of annotation (if your features are e.g. KEGG Orthology groups). In principle, you can then combine information across features as well as across samples. Simplest I can think of, without even looking at the literature, would be to use permutations (shuffle row and column labels) to look for modules different in abundance between groups. (Many existing methods are doing some more sophisticated version of this.)
I think that LEfSe in theory should work quite well here, it's been a while since I thought about it but I believe this is the problem that it is designed to solve, assuming that I understand your question.
There are a lot of methods that are previously described for doing this kind of work. I would suggest reading some reviews of methods, for example the section "Multivariate differential abundance testing" in this article, or here, looking for some inspiration. To some extent it depends on which questions you are interested in answering and what exactly your features represent.


Is there a computational tool or possibility to identify mRNA isoforms from the count matrix of a bulk RNA sequencing dataset? I have the counts matrix of an RNA sequencing dataset of fibroblasts and I wish to identify isoforms of a particular gene of interest in it. Can anyone please hint me on a bioinformatics method to achieve this?

In literature, 9 different isoforms of my gene of interest have been reported. My colleague who does the wet lab experiments has identified 3 of those isoforms through western blotting in the fibroblasts that later got sequenced. There might be more isoforms in these fibroblasts.

I will be grateful for any computational hint to identify gene isoforms from my dataset, irrespective if they are differentially expressed or not.

Thank you in advance for your kind help.

You can back to your bam file and use cuffdiff tools it is very able to do your job follow this to liks
http://cole-trapnell-lab.github.io/cufflinks/manual/

https://github.com/Jeanielmj/bioinformatics-workshop/wiki/Differential-Analysis-with-Cuffdiff

My best


Simulate fast5 MinION run I'm planning on using Oxford Nanopore MinION sequencer to test a script that I'm working on. However, I would not like to run multiple MinION runs to validate the tool. I'd like to simulate various minION runs and get raw reads and all pore events in a fast5 file.  Do you know of any Oxford Nanopore (MinION) sequencing simulator that allows to get fast5 as output?
I looked into NanoSim but it seems like there is no option to get fast5 as an output.

Using a simulator for example MinION output doesn't make sense when there are an abundance of freely-available FASTQ and FAST5 sequences for general purpose use (and especially for the purposes of testing out algorithms). Simulation of sequences can only ever express the components of DNA that are explicitly modelled, which is never going to fully represent reality.
As long as you're doing your own private research, you can use the FAST5 sequences from Clive Brown's own genome:
https://github.com/nanoporetech/ONT-HG1
Bear in mind the license attached to that data:
https://github.com/nanoporetech/ONT-HG1/blob/master/LICENSE.md

License to use for research use only (i.e. non-commercial) governed by the law of England and Wales.
You agree to this License by downloading the data on this site whether in whole or in part.
You cannot re-distribute modified copies of these data, or derivative works, without my express written permission. Where such permission is granted, you must cite the original work and its Copyrights and you must include this license with any work that includes any of my genome.

For a more open license, you can use the nanopore-wgs reads here (derived from human cell line NA12878):
https://github.com/nanopore-wgs-consortium/NA12878

We encourage the reuse of this data in your own analysis and publications which is released under the Creative Commons CC-BY license. Therefore we would be grateful if you would cite the reference below if you do.

If you don't want human data, there's mouse parasite data here:
http://www.ebi.ac.uk/ena/data/view/PRJEB20824
And E. coli data here:
https://www.ebi.ac.uk/ena/data/view/PRJEB22772
... and if you just want 8h of sensor noise as a FAST5 file, I've got that here in one of my DropBox folders (as a ~70MB file):
https://www.dropbox.com/s/20j1anlz0motrs4/MINICOL233_20170118_FN__MN15625_sequencing_throughput_ONLL3127_84599_ch320_read10302_strand.fast5?dl=0
Reminder: if you do ever present the results of analyses based on a data source, it's a good idea to acknowledge the source (even if not expected). If nothing else, it will let the creator know that someone cares about their data, and perhaps encourage them to release more in the future.


What software can I use to develop a population haplotype network? I looking for software to conduct a haplotype network using cytochrome oxidase (CO1) mitochondrial sequence data to study population structure in Anopheles mosquitoes.

The software you are looking for is Network 5.0.0.3

Phylogenetic Network Software . Network © Copyright Fluxus Technology Ltd 1999–2018. Retrieved from https://www.fluxus-engineering.com/sharenet.htm

Whilst its old code it is still regularly used in pop gen for example there guys here recently used it. It will do the job, although personally I think Google Fusion tables are very cool, but that is complicated to use properly. R has some cool stuff in this context BTW, but in context I'd keep it simple (i.e. Network 5)


How to extract column from a matrix matching the another file (sample file) I have two files master file.txt and sample.txt
master.txt
Name    GTEX.1117F.3226.SM.5N9CT    GTEX.111FC.3126.SM.5GZZ2    GTEX.111FC.3326.SM.5GZYV    GTEX.1128S.2726.SM.5H12C    GTEX.1128S.2826.SM.5N9DI    GTEX.117XS.3026.SM.5N9CA
ENSG00000223972 1   0   0   1   0   0
ENSG00000227232 298 168 197 106 221 184
ENSG00000278267 0   1   1   0   0   0
ENSG00000243485 0   2   0   1   0   1
ENSG00000237613 0   0   0   1   0   1
ENSG00000268020 1   0   0   1   0   1
ENSG00000240361 1   3   0   4   2   1
ENSG00000186092 3   1   1   1   1   2

Sample file.txt
GTEX.1117F
GTEX.111FC
GTEX.111XS

Desired output
Name    GTEX.1117F.3226.SM.5N9CT    GTEX.111FC.3126.SM.5GZZ2    GTEX.111FC.3326.SM.5GZYV    GTEX.117XS.3026.SM.5N9CA
ENSG00000223972 1   0   0   0
ENSG00000227232 298 168 197 184
ENSG00000278267 0   1   1   0
ENSG00000243485 0   2   0   1
ENSG00000237613 0   0   0   1
ENSG00000268020 1   0   0   1
ENSG00000240361 1   3   0   1
ENSG00000186092 3   1   1   2
ENSG00000238009 2   2   1   1
ENSG00000233750 0   3   3   0
ENSG00000268903 103 44  76  24

I tried using the R-code form this link
Print specific columns in a matrix on the basis of sample id's in the header
However, I am not able to figure out how to extract columns matching only half part of ID (like for GTEX.111FC I want both GTEX.111FC.3126.SM.5GZZ2 and   GTEX.111FC.3326.SM.5GZYV)
Thank you

Assuming your "master" files looks like this:
df <- data.frame(Name = LETTERS[1:5],
                 Col1 = sample(1:10,5),
                 Col2 = sample(1:10,5),
                 Col11 = sample(1:10,5),
                 Col3 = sample(1:10,5))

  Name Col1 Col2 Col11 Col3
1    A    5    4     3   10
2    B    1    6     5    6
3    C    3    1     9    4
4    D    9    9     8    9
5    E    8   10     1    1

For a sequence of pattern you want to extract, you can use base r function such as grep:
match_pat = c("Col1","Col2")

m <- unlist(sapply(match_pat,function(x) grep(x,colnames(df))))

Col11 Col12  Col2
    2     4     3

And then, extract your columns of interest by doing:
> df[,m]
  Col1 Col11 Col2
1    5     3    4
2    1     5    6
3    3     9    1
4    9     8    9
5    8     1   10

In your example, it will be something like:
match_pat <- c("GTEX.1117F","GTEX.111FC","GTEX.111XS")
m <- unlist(sapply(match_pat,function(x) grep(x,colnames(df))))

master[,m]
```


How can I convert codon coordinates to genomic position? I am looking for a given mutation in IGV, which accepts coordinates in the form of

chr<X>:<Y>

Where X is the chromosome number and Y is the base location in that chromosome.

Eg. chr6:26091179 which has the following attributes:

Gene: HFE|LOC108783645
Marker: rs1799945
Position: chr6:26091179

Additionally, an interval can be written like

chr<X>:<Y>-<Z>

So basically from base located at Y to Z.

My questions are:

1) How do I convert HFE|LOC108783645 to chr6:26091179 (and vice-versa)?

2) How do I convert rs1799945 to chr6:26091179 (and vice-versa)?

3) On a different example, how do I convert something like codon 109 of exon 3 of beta globin gene (GTG->TTG) to chr<X>:<Y> (or chr<X>:<Y>-<Z>)?

Why are there so many ways to specify positions in genes? For me the most clear way is chr<X>:<Y>, all the others leave me guessing...

HFE|LOC108783645 is https://www.ncbi.nlm.nih.gov/gene/108783645. You can use biomaRt to go from the LOCxxxx ID to location. Searching by HFE will give you a different location though, so you might have to account for such eventualities. It is a long and often manually intensive process.
You should be able to use entrez-direct/eutils to go from rsID to position
Unless there's some sort of natural language processing API, codon 109 of exon 3 of beta globin gene cannot be parsed by existing APIs (Unless it's always codon X of exon Y of Z, in which case you could use grep too get the components. Even so, going from exon + codon + codon base change to exact base location will require a lot of wrangling.


Transcript Coordinate Ranges to Genomic Coordinates I have 2 GFF3 files:

Features using transcript IDs as the landmarks. i.e. "CDS" feature types using coordinates from transcript space.
Features using chromosome IDs as the landmarks. i.e. "exon" feature types using coordinates from chromosome space.

I would like to transform the coordinates of the features in file 1 into the coordinate space of file 2. i.e. transcript-based coordinates into genomic-based coodinates.

Here is an example for file1:

cat transcript_orfs.gff3
##gff-version 3
##sequence-region Tx.1 1 4000
Tx.1    ORF_finder  gene    1   4000    .   +   .   ID=1
Tx.1    ORF_finder  mRNA    1   4000    .   +   .   ID=2;Parent=1
Tx.1    ORF_finder  exon    1501    2500    .   +   0   ID=3;Parent=2
Tx.1    ORF_finder  CDS 1501    2500    .   +   0   ID=4;Parent=2

gt sketch -addintrons transcript_orfs.png transcript_orfs.gff3

Here is an example for file1:

cat genome.gff3
##gff-version 3
##sequence-region chr3A_part1 1 454103970
chr3A_part1 genome_assembly gene    1001    6000    .   +   .   ID=Tx
chr3A_part1 genome_assembly mRNA    1001    6000    .   +   .   ID=Tx.1;Parent=Tx
chr3A_part1 genome_assembly exon    1001    3000    .   +   .   ID=Tx.1.exon1;Parent=Tx.1
chr3A_part1 genome_assembly five_prime_UTR  1001    2000    .   +   .   ID=Tx.1.utr5;Parent=Tx.1
chr3A_part1 genome_assembly CDS 2001    3000    .   +   0   ID=Tx.1.cds1;Parent=Tx.1
chr3A_part1 genome_assembly exon    4001    6000    .   +   .   ID=Tx.1.exon2;Parent=Tx.1
chr3A_part1 genome_assembly CDS 4001    5000    .   +   2   ID=Tx.1.cds2;Parent=Tx.1
chr3A_part1 genome_assembly three_prime_UTR 5001    6000    .   +   .   ID=Tx.1.utr3;Parent=Tx.1

gt sketch -addintrons genome.png genome.gff3

I want to convert the coordinates of the features from file1 into chromosome genomic coordinates. I expect to get something like the following output:

cat output.gff3
##gff-version 3
##sequence-region chr3A_part1 1 454103970
chr3A_part1 ORF_finder  gene    1001    6000    .   +   .   ID=1
chr3A_part1 ORF_finder  mRNA    1001    6000    .   +   .   ID=2;Parent=1
chr3A_part1 ORF_finder  exon    2501    3000    .   +   0   ID=3.1;Parent=2
chr3A_part1 ORF_finder  CDS 2501    3000    .   +   0   ID=4.1;Parent=2
chr3A_part1 ORF_finder  exon    4001    4500    .   +   0   ID=3.2;Parent=2
chr3A_part1 ORF_finder  CDS 4001    4500    .   +   1   ID=4.2;Parent=2

gt sketch -addintrons output.png output.gff3

I have looked into using mapFromTranscripts() from Bioconductor's GenomicRanges library but I have made little progress trying to decipher the manual.

This task has the same "flavor" as many I've done before, but each case is so subtly different that it's impossible to write a generalized tool that will work correctly for all circumstances.

R isn't my wheelhouse, but I was able to throw something together pretty quickly with standard Python. This should get you at least 95% of the way to where you're going. :-)

#!/usr/bin/env python3

from argparse import ArgumentParser, FileType
from collections import defaultdict
from re import search, sub

cli = ArgumentParser()
cli.add_argument('txspace', type=FileType('r'))
cli.add_argument('chrspace', type=FileType('r'))
args = cli.parse_args()

exons = defaultdict(list)
for line in args.txspace:
    if line.startswith('#') or line.strip() == '':
        continue
    values = line.split('\t')
    featuretype = values[2]
    if featuretype not in ('CDS', 'exon'):
        continue
    txid = values[0]
    exons[txid].append(values)

for line in args.chrspace:
    values = line.split('\t')
    if len(values) != 9:
        print(line, end='')
        continue
    featuretype = values[2]
    if featuretype not in ('gene', 'mRNA'):
        continue
    print(line, end='')
    if featuretype == 'mRNA':
        txid = search('ID=([^;\n]+)', values[8]).group(1)
        start = int(values[3])
        txexons = exons[txid]
        for exon in txexons:
            exon[0] = values[0]
            exon[3] = str(int(exon[3]) + start)
            exon[4] = str(int(exon[4]) + start)
            exon[8] = sub(r'Parent=[^;\n]+', 'Parent=' + txid, exon[8])
            print('\t'.join(exon), end='')


Can we describe synthetic biological sequence data as belonging to a "species" or "taxon"? For the purpose of testing phylogenetics software, we're using synthetic sequence data (i.e., it's computer generated, and there's no actual biological organism with that e.g. DNA sequence).  For real-world biological sequence data, we would say the sequence belongs to a "species" or "taxon" (or "taxonomic group"), but this feels inaccurate since there is no actual species or taxon.

Question: Can we describe synthetic biological sequence data as belonging to a "species" or "taxon"?

I'm wondering if this is standard practice, despite the synthetic nature of the data.

Phylogenies can be made with/of non-species groups. I think "taxon" is perfectly fine. You might also consider the use of the "leaf / leaves" or "terminal" if you want to emphasise the tip nature of the organisms. Even "sequence" would be fine.

Ninja edit to my own suggestion: since you're making a phylogeny from sequences of unknown context or relationship (i.e. you can't say that they're sequences from a "species" or "OTU" or from anywhere), agnostic terms like "taxon" or "sequence" are more correct.


Why use "robust" FPKMs? Both DESeq2 and edgeR have an FPKM/RPKM function that by default uses normalized library sizes ("robust" option in DESeq2). FPKMs have their own issues, but I thought the main benefit was to have comparable units across independent experiments. If you are adding an experiment-specific normalization factor, then why even use FPKMs?

FPKM are inherently experiment specific and can not be used to compare across samples. Let's consider the following two sequencing runs. Let $E1$ and $E2$ be the true, underlying expression in two samples of genes 1-6. Let $S1$ and $S2$ be the observed expression in our sequencing.
$$
\begin{matrix}
Gene & E1 & S1 & E2 & S2 \\
  G1 & 100 & 10 & 100 & 20 \\
  G2 & 100 & 10 & 100 & 20 \\
  G3 & 100 & 10 & 100 & 20 \\
  G4 & 100 & 10 & 0 & 0 \\
  G5 & 100 & 10 & 0 & 0 \\
  G6 & 100 & 10 & 0 & 0 \\
\end{matrix}
$$
Our totals are 60 counts for $S1$ and 60 for $S2$. We sequenced both libraries to the same depth, but as there are fewer genes expressed in $E2$, we just naturally capture the genes that are expressed more often. For simplicity's sake, let's assume the genes all have the same length of 1, so we ignore the "K" in FPKM. Let's also forget the scaling by a million, since that's just to make nice numbers:
$$
\begin{matrix}
Gene & E1 & S1 & E2 & S2 \\
  G1 & 100 & 0.167 & 100 & 0.333 \\
  G2 & 100 & 0.167 & 100 & 0.333 \\
  G3 & 100 & 0.167 & 100 & 0.333 \\
  G4 & 100 & 0.167 & 0 & 0 \\
  G5 & 100 & 0.167 & 0 & 0 \\
  G6 & 100 & 0.167 & 0 & 0 \\
\end{matrix}
$$
Now let's calculate size factors. In DESeq2, we create a pseudo-sample which is the geometric mean of all counts of each gene, so $\sqrt{10 \cdot 20} = 14.14$ for genes 1-3 and genes 3-6 are ignored as they have a member with 0-valued counts. The size factor for each library is the median ratio for each library compared to that pseudo-sample:
$$
med([\frac{10}{14.14}, \frac{10}{14.14}, \frac{10}{14.14}) = 0.71
$$
$$
med([\frac{20}{14.14}, \frac{20}{14.14}, \frac{20}{14.14}]) = 1.41
$$
So where do we end up if we divide by the size factors?
$$
\begin{matrix}
Gene & E1 & S1 & E2 & S2 \\
  G1 & 100 & 141 & 100 & 142 \\
  G2 & 100 & 141 & 100 & 142 \\
  G3 & 100 & 141 & 100 & 142 \\
  G4 & 100 & 141 & 0 & 0 \\
  G5 & 100 & 141 & 0 & 0 \\
  G6 & 100 & 141 & 0 & 0 \\
\end{matrix}
$$
So in this simple example, size factors produce a better estimate than FPKM or TPM to compare values across experiments, hence the "robustness". Please note that this is still not great for any statistical testing, and in general should be used with caution anywhere.
Finally, I'll leave some more reading material concerning FPKMs and cross-sample comparisons here:

[1] Dillies M.A., Briefings in Bioinf. 2012
[2] Soneson & Delorenzi, BMC Bioinf. 2013
[3] Lior Pachter’s lecture at CSHL (Cold Spring Harbour Laboratory) from minutes 32 on: https://www.youtube.com/watch?v=5NiFibnbE8o&t=32m
[4] From Answer: "Over-correction" in the size-factors of the DESeq2 package: https://support.bioconductor.org/p/80733/#80736
[5] http://www.biomedcentral.com/1471-2105/14/370
[6] https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/
[7] http://rnajournal.cshlp.org/content/22/6/839


Removing PCR duplicates in RNA-seq Analysis After reading some of the forum posts in Biostar and SeqAnswers I find it very confusing whether to filter out the duplicate reads from aligned files or not. As far I understand it's very difficult to distinguish between highly expressed genes and duplicate reads and we may lose important information during the filtration process.

So, is it really necessary to remove the duplicates in differential expression analysis using RNA-seq data?

For normal RNA-seq PCR duplicates are normally kept in, but the duplication rate can be used as a quality control: The higher the duplication rate, the lower the quality. For expression analysis, it is probably best to discard high duplication rate samples, rather than deduplicate them.

In general, the smaller the amount of RNA input into the library prepartion the worst the duplication. Many protocols for very low input quantities (such as single cell) include random barcodes called UMIs (Unique Molecular Identifiers). These allow PCR duplicates to be distinguished from genuinely independent molecules that just happen to have the sample mapping position.


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

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

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

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


Inheritance annotation for nsSNPs I am relatively new to SNP analysis. Is there a database to find annotation about the Mendelian inheritance of SNPs?
I have a small list of nsSNPs (non-synonymous SNPs) and I need to find how they are they passed to from parents to offspring.
Thank you in advance.

If you only have a small number of genes, OMIM is a very easy way to do this.
For example:

If you have too many to do manually, I suspect there is an easy way to query lists, although I'm not sure exactly the way myself. Perhaps someone can add a new answer if so.


Assembly reads having a copy of their beginning in their tail I am analyzing the reads for the SARS-CoV-2 assembly with id SRR11140748. Apparently these reads were obtained with parallel sequencing by Illumina and Oxford Nanopore Technologies.
I have found these reads
GGTAAGTATGTACAAATACCTACAACTTGTGCTAATGACCCTGTGGGTTTTACACTTAAAAACACAGTCTGTACCGTCTGCGGTATGTGGAAAGGTTATGGCTGTAGTTGTGATCAACTCCGCGAACCCATGCTTCAGTCAGCTGATGCACAATCGTTTTTAAACGGGTTTGCGGTGTAAGTGCAGCCCGTCTTACACCGTGCGGCACAGGCACTAGTACTGATGTCGTATACAGGGCTTTTGACATCAGTATGTACAAATACCTACAACTTGTGCTAATGACCCTGTGGGTTTTACACTTAAAAACACAGTCTGTACCGTCTGCGGTATGTGGAAAGGTTATGGCTGTAGTTGTGATCAACTCCGCGAACCCATGCTTCAGTCAGCTGATGCACAATCGTTTTTAAACGAGTATGTACAAATACCTACAACTTGTGCT

As you can see, the second one has an exact copy at the tail of its beginning:
AGTATGTACAAATACCTACAACTTGTGCT

I have found other issues apparently caused by rearrangements of the virus genome, but this case seems to be different. It seems unlikely that the virus mutated in such a way and then a read was cut exactly this way. What is the most likely explanation for this read?
It's not the only one displaying the same or similar behavior. for example this other one has a copy at the tail but from bps not exactly at the beginning (CACCGGAAGCCAATATGGATC):
CAGTTACACCGGAAGCCAATATGGATCAAGAATCCTTTGGTGGTGCATCGTGTTGTCTGTACTGCCGTTGCCACATAGATCATCCAAATCCTAAAGGATTTTGTGACTTAAAAGGTAAGTATGTACAAATACCTACAACTTGTGCTAATGACCCTGTGGGTTTTACACTTAAAAACACAGTCTGTACCGTCTGCGGTATGTGGAAAGGTTATGGCTGTAGTTGTTACACCGGAAGCCAATATGGATC

If you google AGTATGTACAAATACCTACAACTTGTGCT you'll find it's a primer sequence: https://artic.network/resources/ncov/ncov-amplicon-v3.pdf


How to simulate biological pathways? I would like to simulate a biological pathway in Matlab, Mathematica or Python. To study the response of the pathway using Control Theory. Is there a standard way to simulate pathways by using differential equations or any other formalism.

For example I would like to try to simulate this KEGG pathway map. Is there something more recent than this work from 2005?

For example to simulate Petri networks?

In matlab this can be done using the SimBiology toolbox.
In mathematica this can be done using the SystemModeler.
In python there are multiple packages, e.g. (in no particular order) PySB, Tellurium, PySD


Cluster is split in 2-3 locations on tsne plot - Suerat I am running a single cell dataset (count data - exon) through Seurat. After running tsne I see a cluster (13) split in 3 different locations on the plot. Here are the commands I am running:

library('Seurat')
## rows (genes) that have 99% zeros were filtered out,
## and then columns (cells) that have 99% zeros were filtered out
isp.data <- read.table(file = "genes.cells99.percent.filtered", sep =
            ",", row.names = 1, header=TRUE)
isp <- CreateSeuratObject(raw.data = isp.data, project =
       "brain_regions")
isp <- NormalizeData(object = isp, normalization.method =
       "LogNormalize",scale.factor = 10000)
isp <- FindVariableGenes(object = isp, mean.function = ExpMean,
       dispersion.function = LogVMR,do.plot = FALSE)
hv.genes <- head(rownames(isp@hvg.info), 1000)
isp <- ScaleData(object = isp, genes.use = hv.genes, num.cores = 3)
isp <- RunPCA(object = isp, pc.genes = hv.genes, pcs.compute = 100,
       do.print = TRUE, pcs.print = 1:5, genes.print = 5)
PCElbowPlot(object = isp, num.pc = 100)
isp <- JackStraw(object = isp, num.replicate = 100, display.progress =
                 FALSE, num.pc = 30)
JackStrawPlot(object = isp, PCs = 1:30)
isp <- FindClusters(object = isp, reduction.type = "pca", dims.use =
        1:26, resolution= 0.75, print.output = 0, save.SNN = TRUE,
        force.recalc=TRUE)
isp <- RunTSNE(isp, dims.use = 1:26, do.fast = T)
TSNEPlot(isp, do.label = T, pt.size = 0.5)

Here is the tsne plot:

I have tried varying parameters for resolution in FindClusters (from 0.5-3.5), but that just increases the clusters found. Further I tried to play around with parameters min.dist (0.75,1), n_neighbors (10 to 30) and perplexity(5 to 50) but I still get at least 1-2 clusters split in the tsne plot.

Another thing I have tried is to reduce the dimension (dims.use=10) on which tsne is ran but the clusters seem to be too close together, like so:

Note: Dimensionality of dataset (after filtering) - 23824 genes across 8967 samples.

My question being why does such behavior happen? Is my filtering strategy not correct? Have I not normalized/scaled data correctly?

Your cluster labels come from graph clustering implemented in the FindClusters() function. The resulting clusters are then visualised with a 2D tSNE plot (via RunTSNE() and TSNEPlot()). So, your cluster 13 is not split into three sub-clusters but cells within cluster 13 look somewhat distant from each other on the tSNE plot. Having mentioned distances on a tSNE plot, please see this SE Data Science entry Can closer points be considered more similar in T-SNE visualization? and especially the comment "What I mean is that you can't just use distance in the lower space as a similarity criterion". In your case "the lower space" would be your tSNE plot.

As far as I can see from your code, you are following the widely-adopted "Seurat workflow" and all seem fine including your TSNE plot output. I could not spot the QC step though: apart from the one at object initiation, one usually applies filters for number of genes and % mitochondrial gene contribution (using FilterCells() in Seurat v2 and subset() in Seurat v3).

It would be a good idea to try UMAP instead of tSNE, please see Evaluation of UMAP as an alternative to t-SNE for single-cell data (Becht et al, 2019). Seurat v3 vignettes are now using UMAP and not tSNE any more.

Lastly, please see How to Use t-SNE Effectively to learn about points to take into account when using/interpreting tSNE plots. You will see that you would get surprisingly dissimilar tSNE plots from the same underlying data with different parameters (some are hardcoded in the RunTSNE() function).


How can I obtain a list of all NCBI gene ID's along with their full name, symbol, and also known as symbols? I would like to do two things:
(1) Obtain a list of all NCBI (Entrez) gene ID's.
(2) Obtain the "Official Symbol", "Official Full Name", and "Also known as" fields from the NCBI gene pages for each of those genes.
This information, since it is displayed on the NCBI website, must be located in some databases, but I have not had luck finding a tool that will allow me to extract this information.

All of these data are in the gene_info.gz file on the NCBI FTP site: https://ftp.ncbi.nlm.nih.gov/gene/DATA/
Specifically, you may want to look at the gene_info.gz file. It has the following fields (shown for an example gene below):
                              #tax_id [  1]: 9606
                               GeneID [  2]: 5768
                               Symbol [  3]: QSOX1
                             LocusTag [  4]: -
                             Synonyms [  5]: Q6|QSCN6
                              dbXrefs [  6]: MIM:603120|HGNC:HGNC:9756|Ensembl:ENSG00000116260
                           chromosome [  7]: 1
                         map_location [  8]: 1q25.2
                          description [  9]: quiescin sulfhydryl oxidase 1
                         type_of_gene [ 10]: protein-coding
   Symbol_from_nomenclature_authority [ 11]: QSOX1
Full_name_from_nomenclature_authority [ 12]: quiescin sulfhydryl oxidase 1
                  Nomenclature_status [ 13]: O
                   Other_designations [ 14]: sulfhydryl oxidase 1|quiescin Q6 sulfhydryl oxidase 1|testis tissue sperm-binding protein Li 62n|thiol oxidase 1
                    Modification_date [ 15]: 20200604
                         Feature_type [ 16]: -

If you are starting  with a set of NCBI GeneIDs and want this sort of information for just those, you can use Entrez Direct for that. If you are starting with RefSeq transcript or protein accessions, be sure to check out the other files in the FTP path mentioned above.


pvalue calculation of protein-protein network with permutation test Note: this question has also been asked on Biostars

I have performed a network recreation analysis based on the interactions of proteins from String db.

I wanted to know if the interactions have significance by reshuffling the proteins with permutation against the background network. Since, its first time for me to perform permutation test on PPI, Can someone guide me to a tool or a tutorial to perform permutation test on newly created PPI network.

I tried looking at a youtube tutorial and read some articles like this one, but they were not very useful for my purpose.

It is probably necessary to correct for the degree distribution of the network, I am not sure that a tool such as the Mantel test will quite do this (as suggested by user Michael G.). In naturally occurring networks degree distribution is a strong confounder.

Then again I am not totally sure what shape your data is in or what the question is. You mention a background network, but I am not sure what this is or how you want to use it. Possibly your network reconstruction is a subset of your background network? In which case you just define the rule for how your subsample is drawn and then you compare it to permuted networks.

However, here is some fairly crude but serviceable R code that implements a general purpose network rewiring test, using igraph. I am not using a background network as again I'm not sure what it is.

(Python also has igraph so should be similar rewire() function)

compute_measure = function(network) {
  # i don't know what you are measuring, so you have to write this function.
  # this function takes a graph object as input
  # this function returns a single numeric value (whatever you're measuring) as output
}

# "net" is a network of the graph type from igraph

measure = compute_measure(net)

num_perms = 1000

perm_measures = c()

for (i in 1:num_perms) {
  # permute the network by rewiring, rewire 10x as many times as edges
  # to ensure saturation
  perm_net = rewire(net, with = keeping_degseq(niter = gsize(net) * 10))
  perm_measures[i] = compute_measure(perm_net)
}

# one-sided monte carlo p value is proportion of rewirings (plus one) greater than your measure
p_value = (1 + length(which(perm_measures >= measure))) / num_perms

If on the other hand you are just trying to come up with a single permutation of the edges of your network to compare against, you can run rewire() only once to get that permuted network. It depends on what question you are trying to answer, and what null hypothesis you're trying to get a p value for.


Finding a single open reading frame with ribosomal binding site, using Biopython I'm given a Fasta file, containing a large DNA(over 115,000 long) sequence, and I am tasked with finding a single large open reading frame contained within the DNA sequence using Biopython.

I'm aware this has been asked before consistently however i can't find anything for actually finding the source of my ORF, but apologies for any overlap.

Now from other sources and the Biopython cookbook I've translated my sequence and found six open reading frames (three for each strand) and their positions within the sequence;

def find_orfs_with_trans(seq, trans_table, min_protein_length):
    answer = []
    seq_len = len(seq)
    for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
        for frame in range(3):
            trans = str(nuc[frame:].translate(trans_table))
            trans_len = len(trans)
            aa_start = 0
            aa_end = 0
            while aa_start < trans_len:
                aa_end = trans.find("*", aa_start)
                if aa_end == -1:
                    aa_end = trans_len
                if aa_end-aa_start >= min_protein_length:
                    if strand == 1:
                        start = frame+aa_start*3
                        end = min(seq_len,frame+aa_end*3+3)
                    else:
                        start = seq_len-frame-aa_end*3-3
                        end = seq_len-frame-aa_start*3
                    answer.append((start, end, strand,
                                   trans[aa_start:aa_end]))
                aa_start = aa_end+1
    answer.sort()
    return answer

orf_list = find_orfs_with_trans(record.seq, table, min_pro_len)

for start, end, strand, pro in orf_list:
print("%s...%s - length %i, strand %i, %i:%i" \
      % (pro[:30], pro[-3:], len(pro), strand, start, end))

OUTPUT

 TNRQVYGGTLQSLRTGTGIYSRLASSPTNR...CEA - length 160, strand 1, 7950:8433
 GIPPGRTEGLGRYVHGESIFLEATLVPEPQ...VQS - length 171, strand -1, 19275:19791
 ISVHLRRYFSVLLRAPVALNADSRVAGPLD...ITP - length 190, strand 1, 34079:34652
 GCVNGNFPDYRVAVDDPGALVVGGELGQAL...VNT - length 771, strand -1, 39335:41651
 LLKLLASRLTWLLVLSAALGFLTSVSYRLG...SPG - length 235, strand -1, 39358:40066
 PGDVIVSKVPADNLRPRMSEINDFTNSIII...VLV - length 826, strand 1, 39362:41843
 GPGRQSSSAYERDQRLHQQHHHQRERHQVR...QRS - length 764, strand 1, 39385:41680
 LLLSSTSSASCWLKKFFFRLLLLSRSVALL...LAL - length 158, strand -1, 40459:40936
 TRLIIRPPSHTTLYYDPYPQPLKVSSSVPC...RIT - length 159, strand -1, 54663:55143
 IGAFSRGYNQGWVLIEAQEVRTTSPLFPRS...RDQ - length 201, strand 1, 56104:56710
 KARGSSLGAIYCRRLKRRSLPSWQIESFNC...TNS - length 150, strand -1, 62497:62950
 RRPARAGAELNKDDIRDTCLLSCSEVDNKV...VTL - length 168, strand 1, 72339:72846
 SGPTTVRTSAPQRCWRTSLKPKLSLSDSTE...PVA - length 177, strand -1, 114494:115028
 ERELWLETGPPTPLWGTGPDCSRAALNRQR...PRR - length 157, strand 1, 114953:115427

However I'm unsure if these open reading frames contain the ribosome binding site(consensus sequence of AAGGAGGTG and occurring between 6 to 9 nucleotides upstream from the positive 5' strand). Is there anyway to add to my code to include this? or would it already be included with Biopython?

My specification says to find the single open reading frame however since none of these begin the the M start codon i'm a bit lost how they can be so, I have no information on where my sequence is from (something I have to find out) so using the standard translation table 1.

You should start by finding all instances of AAGGAGGTG and then look for ORFs that begin 6 to 9 nucleotides downstream of them. By begin, I mean find an AUG codon, and then extend until the next in-frame stop codon.

There's also no point in translating, just look for long stretches of sequence between AUG and UGA, UAA or UAG, those will define your reading frames. You seem to be starting the reading frame at any non-stop codon instead of starting from AUGs. This may well make biological sense in some contexts (there are proteins that don't start with M) but is probably beyond the scope of what sounds like a simple assignment.


Find indels between two short sequences I have two sequences, say AAAGCTCGAGG and AAAGCGAGG. I need a convenient tool which shows me insertions or deletions between these, i.e. in this case something like

AAAGCTCGAGG
|||||  ||||
AAAGC  GAGG
2bp deletion starting at base pair 5

and also insertions such as

AAAGCTC GAGG
||||||| ||||
AAAGCTCTGAGG
1bp insertion at base pair 7

Does something like that exist? Blastn finds the alignment of course, but only outputs the number of gaps, not their position. I guess I could parse a pileup, but do not want to create a bam file from one of the sequences first, so samtools seems to be out. Any further ideas?

Smith-Waterman is guaranteed to give you the best alignment given your choice of substitution matrix.  For a few such tiny sequences, it's what I would use.  Many implementations exist.


How to efficiently compute the exact percentage of non-unique k-mers in a genome for given k? I'm looking for some software that can "efficiently" (time and memory) compute the exact percentage of non-unique k-mers in a genome for given k. I don't need the k-mers or the abundances itself, I just need the percentage.

Alternatively, the result could be:

the number of different k-mers observed in the genome
the number of unique k-mers, i.e. k-mers that occur only once in the genome, or the number of non-unique k-mers, i.e., k-mers that occur more than once in the genome.

which can be easily used to compute the requested percentage.

Requested features:

don't count k-mers with ambiguous nucleotides
count a k-mer and its reverse complement together
k upto 50 or 100
for large genomes (upto some GB)
efficiency

For the study in mind, I like to vary k. Hence, I will run the software several times with different k.

Any help is welcome.

You can do all of that with khmer. For example, abundance-dist-single.py produces a file with columns: k-mer abundance, k-mer count, cumulative count, and fraction of total distinct k-mers. So for question 1 you would sum column 2. For question 2 you would just get thek-mer countassociated with ak-mer abundance` of 1.

That package also provides a python API (see the read the docs link above) if you need to customize things more.


What are the different kinds of bioluminescent genes? I know of the common green glow gene but I forgot the name and I also know that some algae glow blue. There are so many types of bioluminesent organisms, so I am wondering what species have which genes and what other genes are associated with their metabolism and the systems that regulate these genes.

There are three processes wherein light is emitted.

Bioluminescence: a chemical reaction releases light. The enzyme that does this is luciferase, while its substrate is luciferin, a small molecule. There are half a dozen organisms that have analogous (evolved independently) enzymes and totally different substrates, such as fireflies, Alivibrio bacteria, diatomes, deep sea fish, some mushrooms etc. The wiki article has a list of these and their small molecules.
Fluorescence: light is absorbed by the molecule and is emitted at a lower energy (wider wavelength, redshifted). Green fluorescent protein of Aequorea victoria jellyfish is the classic example of this. It is excited with UV and emits green light. Here is an interactive model, I threw together. Variants exists for biotech reasons that span the whole of the visible spectrum. Other species have a homologous beta barrel, such as corals (red).
Phosphorescence: light is absorbed by the molecule and emitted later.  Ignoring the picosecond delay in fluorescent small molecules (cf. fluorescent lifetime microscopy), this is exclusively a properties of lanthanide elements (e.g. Europium). Lathanide binding domains have been designed, but do no exist in nature.


where to find pathogenicity metadata for E.coli genomes I have a list of about 900 E.coli genome ids (Genbank plus NCTC ids), e.g.:

NCTC8196_1
NCTC11105_1
CP027464.1
CP027459.1
CP027363.1
CP027388.1
CP027355.1
CP027675.1
CP027445.1
CP027325.1
CP027219.1

and I would like to find pathogenicity metadata for them so that I can classify the molecular data I found with the metadata.

Is there a way to get a table of metadata from these genomes that I can use to associate with my molecular data?

NOTES:

After some searching online, I found Enterobase (http://enterobase.warwick.ac.uk), but it seems the site is mostly dedicated to MLST typing, thus they use different naming conventions than the Genbank id. Maybe there is a way to link the info in Enterobase with Genbank ids somehow?

It seems like Enterobase contains project ids, e.g. PRJNA218110. These project ids can also be found on the other side in the GenBank entries, for which we could associate the GenBank chromosomal ids.

Searching for PRJNA218110 in Enterobase gives a very long list of entries that seem to be from Illumina sequencing data, but in this list, the 'Path/Nonpath' column is always empty.

Enterobase seems to contain lots of entries with useful information, such as: the strain is pathogenic, and produced diarrhea in humans. So a table with entries like the one below would be very useful:

Genbank_id: CPXXXXXXX.1
Path/Nonpath: Pathogenic
Disease: Diarrhea

Genbank_id: CPYYYYYYY.1
Path/Nonpath: Nonpathogenic
Disease: NA

you can use BacDive. This is the largest database for bacterial metadata.:
https://bacdive.dsmz.de/
e.g. you can search for the sequence accession "NCTC8196" and you would get the metadata of the according strain:
https://bacdive.dsmz.de/search?search=NCTC8196

What kind of data is available, is highly dependent on the strain. Usually type strains are more rich on data that others.

For further questions don't hesitate to contact us:
contact@bacdive.de


What is the meaning of the Hmmer profile? I want to understand what the hmmer Profile used in hmmer.

consider the profile

HMMER2.0  [converted from 3.1b2]
NAME  aRNH_profile
LENG  121
ALPH  Amino
RF    no
CS    no
MAP   yes
NSEQ  34
DATE  Tue Apr 26 21:24:57 2016
XT       -8455     -4  -1000  -1000  -8455     -4  -8455     -4
NULT        -4  -8455
NULE       656  -1722     98    419   -333    475  -1125    239    250    947  -1073   -271    -50   -338    114    451    113    430  -2131   -717
HMM        A      C      D      E      F      G      H      I      K      L      M      N      P      Q      R      S      T      V      W      Y
         m->m   m->i   m->d   i->m   i->i   d->m   d->d   b->m   m->e
           -6      *  -7932
     1    -214   -916   -565    694    862  -1448    512   -838    115  -1042    170   1348  -1481    369   -158   -290   -159   -625    -52   1675     1
     -    -210   -336    220    -37   -342    377     73   -669    210   -511   -732    404    423     70     27    439    205   -414   -161   -176
     -     -18  -6890  -7932   -894  -1115   -701  -1378    -18      *
     2     285    990  -2213  -1678   2305  -1819   -698    219  -1289   -140    545  -1167  -1835   -973  -1322   -424   1025    762   -559    377     2
     -    -210   -336    220    -37   -342    377     73   -669    210   -511   -732    404    423     70     27    439    205   -414   -161   -176
     -     -18  -6890  -7932   -894  -1115   -701  -1378      *      *
[....]

what exactly does each component mean? By definition a HMM $\lambda$ consist of an alphabet $V$ in my case twenty amino acids. Different states $S$ in our case main, insertion and deletion (beginning and end). And two probability matrices of state transitions and emission in certain state.

now if I look at the profile I can see that it has the allowed transitions
m->m   m->i   m->d   i->m   i->i   d->m   d->d   b->m   m->e
and all the amino acids, but what does each row in the table mean. I checked the manual from the documentation and it says rather ironically

  If you were to look at this file (and you don’t have to – it’s
  intended for HMMER’s consumption, not yours), you’d see [part of my output]

are these number translations of the probabilities in $[0,1]$  into another interval $[-x,+y]$? what is the meaning of the three rows under A 1

-124
-210
 -18

Is this the transition matrix? It would be great if someone could explain this  or point to the correct reference.

The tutorial in the documentation does indeed state what you quoted but a full description of the file format is given in Section 8 of the documentation, starting on page 106.
To summarise, as you noted, the number columns A-K is the alphabet size (the number of amino acids). The number of rows underneath is the number of nodes, which corresponds to the number of match states. The three numbers (-214, -210 and -18 in your example) are: the match emission, insert emission, and state transition values.
Caveat: The HMMER profile HMM format changes frequently between HMMER versions (even between minor versions) so there is no guarantee this information is correct for all the older file formats.

Quoting the relevant part of the documentation on page 109:

main model section
All the remaining fields are mandatory. The first
two lines in the main model section are atypical.‡ They contain
information for the core model’s BEGIN node. This is stored as model
node 0, and match state 0 is treated as the BEGIN state. The begin
state is mute, so there are no match emission probabilities. The first
line is the insert 0 emissions. The second line contains the
transitions from the begin state and insert state 0. These seven
numbers are: B → M1, B → I0, B → D1; I0 → M1, I0 → I0; then a 0.0 and
a ’’, because by convention, nonexistent transitions from the
nonexistent delete state 0 are set to log 1 = 0 and log 0 = −∞ = ‘’.
The remainder of the model has three lines per node, for M nodes
(where M is the number of match states, as given by the LENG line).
These three lines are (K is the alphabet size in residues):
Match emission line
The first field is the node number (1 . . . M). The
parser verifies this number as a consistency check (it expects the
nodes to come in order). The next K numbers for match emissions, one
per symbol, in alphabetic order. The next field is the MAP annotation
for this node. If MAP was yes in the header, then this is an integer,
representing the alignment column index for this match state
(1..alen); otherwise, this field is ‘-’. The next field is the CONS
consensus residue for this node. If CONS was yes in the header, then
this is a single character, representing the consensus residue
annotation for this match state; otherwise, this field is ‘-’. The
next field is the RF annotation for this node. If RF was yes in the
header, then this is a single character, representing the reference
annotation for this match state; otherwise, this field is ‘-’. The
next field is the MM mask value for this node. If MM was yes in the
header, then this is a single ’m’ character, indicating that the
position was identified as a masked position during model
construction; otherwise, this field is ‘-’. The next field is the CS
annotation for this node. If CS was yes, then this is a single
character, representing the consensus structure at this match state;
otherwise this field is ‘-’.
Insert emission line
The K fields on this
line are the insert emission scores, one per symbol, in alphabetic
order.
State transition line
The seven fields on this line are the
transitions for node k, in the order shown by the transition header
line: Mk → Mk+1, Ik, Dk+1; Ik → Mk+1, Ik; Dk → Mk+1, Dk+1. For
transitions from the final node M, match state M + 1 is interpreted as
the END state E, and there is no delete state M + 1; therefore the
final Mk → Dk+1 and Dk → Dk+1 transitions are always * (zero
probability), and the final Dk → Mk+1 transition is always 0.0
(probability 1.0).


How to make work programs from the $PATH? I am trying to analyze my RNAseq reads for defective genomes and I use this program (http://www.di-tector.cyame.eu/) that is suggested by Beauclair et al (https://pubmed.ncbi.nlm.nih.gov/30012569/)  for that.

After downloading a .py script and running it with the following command:

sudo python3 DI-tector_06.py /mnt/e/nastya/SLX066-02/sequence.fasta /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz

It says at the beginning that

=================================
Program: DI-tector
Version: 0.6 - Last modified 25/05/2018
=================================
Requirement: (must be in your $PATH)
-bwa
-samtools

Optional: (must be in your $PATH)
-bedtools

I have tried to install bwa using:

sudo apt-get update -y
sudo apt-get install -y bwa

and samtools using

cd /usr/local/bin
sudo wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
sudo tar -vxjf samtools-1.9.tar.bz2
cd samtools-1.9
sudo make

I did not get exactly what does that mean: (must be in your $PATH) but I have searched here and there and discovered that it might mean /usr/local/bin

The program starts and, after some calculations it gives me the following error:

Input file:               /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz
Host reference:           None
Virus reference:          /mnt/e/nastya/SLX066-02/sequence.fasta
Remove segment:         < 15 nt
Remove reads with MAPQ: < 25
Allow InDel. length:    > 1 nt

=================================
Reading gziped input file...
=================================
Number of reads in input file: 18,035,544
=================================

=================================
Step 1/5 : Alignment against Viral reference...Filtering...
=================================
[E::bwa_idx_load_from_disk] fail to locate the index files
Traceback (most recent call last):
  File "DI-tector_06.py", line 183, in <module>
    Popen(["samtools", "view", "-Sb", "-F", "4"], stdin=out_file_onVirus, stdout=out_file_Virus).communicate()[0]
  File "/usr/lib/python3.6/subprocess.py", line 709, in __init__
    restore_signals, start_new_session)
  File "/usr/lib/python3.6/subprocess.py", line 1344, in _execute_child
    raise child_exception_type(errno_num, err_msg, err_filename)
PermissionError: [Errno 13] Permission denied: 'samtools'

I have no idea how samtools must be accessed, I tried sudo command with python3 but it did not work.

echo $PATH will show the directories in your search path. Either copy the samtools executable there or add its location to the $PATH variable. You would be well served with a unix command line course or tutorial.


Annotation with Prokka or RAST. I was experimenting Prokka and RAST annotation tools. So, I took a well-annotated swinepox virus genome from genebank (NCBI Reference Sequence: NC_003389.1).

I ran those sequences on Prokka and RAST Seed server at the same time. I can see that only a few (may be around 1%) of the genes were annotated. Most of them were predicted as hypothetical protein. And the results were comparable between Prokka and RAST.

I would assume that these tools look for similar sequences in NCBI and find the best-match protein. But looks like that is not the case. It should be able to find that well annotated swinepox virus genome in genebank and predict most of the proteins.

Also, if almost all the genes are predicted as hypothetical protein then there is not much difference between gene prediction tool like Genemark and genome annotation tool. Are there any better annotation tools ? Or this is what we get ? Or have I misunderstood the concept of annotation ? Please someone help me to understand this thing.

I have attached the image for comparison. Left one the Swinepox genome in .gb format and right one is the same genome annotated with Prokka.

Both of the annotation tools you used are designed for prokaryotic genome annotation. I would not expect them to work very well for viruses.


filter out FASTQ reads which are shorter I tried to filter out FASTQ reads which are shorter than 259 bp with

bioawk -cfastx 'length() >= 259 {print "@" " " "\n""\n+\n"}' good/SZ005_NoIndex_L002_R1_009.good.fq.gz \
  | gzip >  good-filtered/SZ005_NoIndex_L002_R1_009.good-filtered.fq.gz`.

Unfortunately, the output was wrong:

@

+

@

+

@

What did I miss?

Thank you in advance,

You need to use variable $seq. For example:

bioawk -cfastx 'length($seq)>=259{print "@"$name"\n"$seq"\n+\n"$qual}' test.fq.gz

There are also more convenient tools like seqtk and seqkit. With those, you may

seqtk seq -L260 test.fq.gz


Manta SV caller: I can't find a file in the path specified by the user guide I trying to locally reproduce a Manta run as seen in the log of an illumina basespace RNA-Seq Alignment tool v1.1.1 run I've done. The log reads as follows:
[LOG][2021-03-12 20:58:56.377]03/12/2021,20:58:56.193,Launch process:  \
/path_to_RNA_Tools/bin/python  \
'/path_to_RNA_Tools/Manta/bin/configManta.py'  \
--rna --normalBam='path_to_alignments/alignments.bam' \
--runDir='path_to_samples/CTR/replicates/CTR/MantaFusions/manta.run'  \
--referenceFasta='/path_to_genomes/genome.fa' \
--config '/opt/illumina/Isis/2.6.25.18/Workflows/WholeGenomeRnaSeqWorker/scripts/configManta.py.ini'

The file shown in the line above configManta.py.ini is however no were to be found.
The Manta user guide states in the linked section

Options listed in the file: ${MANTA_INSTALL_PATH}/bin/configManta.py.ini

The configManta.py is in the expected directory, only the configManta.py.ini is absent.
I am using Manta version  1.6.0 from bioconda.
Your insights will be much appreciated!

The --config file option for Manta can be used to override the defaults in the global config file. If you type configManta.py -h, you'll see where the global config file has been installed. For me, it's here:
/home/steve/miniconda3/envs/manta/share/manta-1.6.0-1/bin/configManta.py.ini

So the command in your log file takes a specific config file installed into:
/opt/illumina/Isis/2.6.25.18/...

I suspect 'Isis' is a reference to Illumina's analysis tool 'IsisRNA'. Papers that mention Illumina's 'RnaReadCounter' tool, say that it is part of the internal analysis tool 'IsisRNA'. Because it's one of Illumina's in-house tools, I think you'd need to contact Illumina support to get a copy. I would only do this though if you need to reproduce exactly as per the command in your log file. If you only need to roughly reproduce the output, the defaults in the global config file are probably fine.


Interpretation of the violin plots from sc-RNA-seq I'm confused about the meaning of the black dots and the red shape in the violin plots from the seurat tutorial:

The black dots represent the values for individual cells. The "nGene" plot (the first one) shows the number of detected genes for every cell.

The red shape shows the distribution of the data. For the "nGene" plot, you can see that the average number of genes per cell is about 900 and most of the cells have roughly around 700-1100 genes.


Why is SMALT better for microbial genomics than other mappers? SMALT seems to be one of the most used read mappers for bacterial data, see, e.g., this query. I do not say that it is not a great mapper, but I cannot easily see what are its main strengths compared to mappers such as BWA-MEM, Bowtie2, NovoAlign or GEM. Moreover, it is not even published.

Could you name some of its distinguishing features (e.g., user support by Sanger Pathogens)?

So far I have heard only arguments like "We use SMALT because everyone does it.", but this is not convincing enough for me.

I have used smalt a couple of times. Smalt is more flexible than bwa/bowtie and can be tuned to be more sensitive to divergent hits, which is useful to certain applications. I heard Sanger did evaluate a few mappers for some non-typical applications (not sure what). They found smalt to perform better. Also smalt was developed at Sanger. I guess they could modify smalt based on their needs, though I don't know the details.


Bootstrapping cell type identifications (scRNA-seq) I have scRNA-seq data held in a SEURAT object, and cell-types were identified with SingleR (essentially a vector of strings - Monocytes, Granulocytes etc.). I have two conditions - treated and untreated. I want to know if the difference in the relative frequency (%) of a given cell type between the conditions is significant. For this, I thought to perform bootstrapping, to obtain frequencies for many different 'simulated' datasets - and then perform a statistical test. However, I cannot understand how best to do this. Should I subsample the scRNA-seq data, with replacement and calculate the frequencies each time? Can somebody provide a step by step workflow on how to do this?

I would use permutation:
i) Get the difference that you are interested in, i.e. number/proportion of monocytes in one group vs the other.
ii) Get the relevant columns (condition, cell label, ...) from the @meta.data slot of your seurat objects and rbind().
iii) Shuffle the cell labels randomly and split the data frame into two (per condition).
iv) Compare the difference between the cell counts/proportions of interest from the two permuted data frames and save this value.
iv) Repeat iii and iv for 100s of times and calculate how many times have you observed a difference at least as extreme as your real observation (i). Once divided by the number of trials, this would give you a p-value stating how likely it is that you get the observed value.
I have to add that batch effects are severe in scRNA-seq and hence a linear model where you might be able to account for batches could be more helpful.


How are residuals used as new phenotype? I'm a novice trying to wrap my head around with fancy jargons. If I search residual, I understand that it represents the difference between the observed value of the dependent variable (y) and the predicted value (ŷ).

I'm unable to follow then, how are these residuals used as new phenotype?

I found similar question in the forum:  Adjusting phenotypes by regressing out covariates

The origin of my question is predixcan pipeline, where residuals are used as the new gene expression values after regressing input values on covariates.

The residual is the level of error in a regression model, the lower the residual the better the model. Residuals cannot equate to a phenotype, it is actually the opposite if your regression model is trying to investigate a given genotype/phenotype(s) model. The residual in this case is the amount of infomation which the genotype/phentype does not account for, i.e. there are other factors resulting in your observations which you have not captured in your model.

What I assume is that the regression model you are referring  is  investigating whether the residual is lowered by examining different gene combinations. The objective of a regression model is to lower the residual and I suspect that it is comparing the new residual against the old one. Therefore it is not being used as an input into a regression model but an assessment of its robustness.


How should I address batch effects in my experiment? Let's say I have an RNA-Seq experiment, where I'm interested in the significantly differentiated genes between pre-treatment and post-treatment conditions. "rep" == biological replicate.

Sample PreA (3 reps **all** in 2017) vs Sample PostA (3 reps **all** in 2018)
Sample PreB (3 reps **all** in 2017) vs Sample PostB (3 reps **all** in 2018)
Sample PreC (3 reps **all** in 2017) vs Sample PostC (3 reps **all** in 2018)

I also have a sample:

Sample PreD (3 reps all in 2017) and Sample PreD (3 reps again in 2018)

It is safe to assume the samples in sequenced in 2017 were done differently to 2018.

Q1: Is there problem with batch effects in my study? When I do a test, there should't be any way to split variation from batch effects and biological variation?

Q2: If Q1 is correct. How should I reduce batch effects?

Can I add a batch factor to my model in DESeq2/edgeR?
Should I use my PreD control samples to derive normalisation factor as in RUV? (assume PreD as a negative control)
Should I use sva package?

EDIT

Devon suggested for adding a batch variable, but how? That's how I'm doing (ignoring batch effects):

dds <- DESeqDataSetFromTximport(txi, data, ~Sample)

data holds a data frame with a Sample column. Sample would be "PreA" and "PostA" in the first test. Very simple.

I couldn't add Date to it, like:

dds <- DESeqDataSetFromTximport(txi, data, ~Date + Sample)

because the code crashed complaining of the 2017/2018 problem discussed above. I'm not sure how I should add PreD to the design.

Yes, though thankfully Your 2018 PreD samples will help you resolve this.
Simply add the batch effect to the design (~Batch + Treatment) and DESeq2 (or edgeR or Limma) will handle this for you.

You do not need SVA or RUV, thankfully, since you quite cleverly sequenced one group in both batches.

To clarify, your coldata will be something like:

Group    Time    Batch
A        Pre     A
A        Pre     A
A        Pre     A
A        Post    B
A        Post    B
A        Post    B
...
D        Pre     A
D        Pre     A
D        Pre     A
D        Pre     B
D        Pre     B
D        Pre     B

And the design is actually ~Batch + Group + Time. You can probably add an interaction term in if you then modify the model matrix manually to ensure it's full rank.


scRNA-seq,10x cellranger pipelines,low custom tdtomato gene content! looking for help！ I have a set of scRNA-seq samples expressing TdTomato, which has high content in microscope. I followed the 10x cellranger pipelines to finsh the work, my procedures are as follows:

added TdTomato on an additional chromosome:

added an entry in the gtf
cellranger mkref
cellranger count

I got the result, which includes:

But strangely, the tdtomato expression content in my three sample are all very low! And cells with tdtomato are also very low and disperse!! Take sample 1 as example, you can see the picture:

It's not normal, but I can't find the reason. Following the 10x genomics help documentation noticed that by executing the CellRanger filtering steps, maybe my tdtomato reads are filtered too much, but because I'm not very skillful at programming and I'm first time to do RNA-seq (and scRNA-seq), I just don't know how to analysis the bam file as suggested by the official website and find the real reason of the problem:

If you can show me the code need to do, I will be very grateful!

A few suggestions:

Manually inspecting the BAM file with a genome browser like IGV, as @GWW suggested and checking the qualities of the alignments in your BAM/SAM file (SAM is the human readable version of BAM, there are tools to convert one to another);
Using another single cell analysis tool like Seurat so that you can manually set thresholds for QC and not lose to many (tdtomato) expressing cells at the QC step (I don't know about your experiment but if transfected, cells could be at a worse condition than non-transfected ones and hence might not pass QC);
Double checking your tdtomato reference sequence, I don't know about tdtomato but for example GFP has flavors like EGFP, similar but different sequences, might explain low alignment score if you observe in step 1;
Taking the limitations of scRNA-seq (especially droplet-based) technologies into account. The figure below is from the CITE-seq paper (from supp info), Simultaneous epitope and transcriptome measurement in single cells (Stoeckius et al., 2017), showing that the protein - RNA measurements of a given gene might show poor correlation.


what percentage of the human genome is MAPQ=0? When doing Illumina 2x150bp sequencing of genomic DNA, and after aligning the reads to GRCh38, what percentage of the non-N fraction of the human genome is MAPQ=0? This is, what part corresponds to regions that can't be uniquely mapped with 2x150bp reads.

And, how many genes are affected by MAPQ=0 regions?

I presume the numbers will dance around depending on fragment sizes, read quality, etc. but I am happy with some starting numbers.

Mapping quality is determined by the repetitiveness of the genome, the sequencing error rate, insert size, the capability of the mapper and the nasty heuristics behind the mapper. MAPQ=0 to one mapper is not necessarily MAPQ=0 to another.

That said, I get what you mean. You want to know the uniqueness/repetitiveness. It is still hard if you want to get a useful answer. For 150bp reads, each reference position is covered by 150 reads. What if 50 of them have no other exact hits elsewhere, but the rest 100 have? Is this position a repeat or not? In addition, what if the 50 are unique only because one mismatch? If there is a variant at that mismatch, the 50 would become repeats or mapped elsewhere.

My preference is to say a position is "unique" under k-long reads if over k/2 reads overlapping the position have no other perfect or 1-mismatch/1-gap hits elsewhere in the genome. Under this definition, 79.3% of human genome are unique for 35bp reads. 92.4% for 75bp reads. I don't have the number for 150bp reads. I guess it will be around 95%. Empirically, 94-95% of human genome is callable with 100bp paired-end reads.

As to other measurements, the most common one is the fraction of reads that has an exact hit elsewhere. You can use Fabio's method. It gives a good enough estimate once you simulate over 1 million reads. This fraction is around 86% for 35bp reads and 95% for 75bp reads, as I remember. The problem with this approach is this fraction is not very informative to variant calling due to the issues I talked about. Another way is to use RepeatMasker. It is worse. RepeatMasker masks 50% of human genome, but excludes segmental duplications where short reads can't be confidently mapped.


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

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

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

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

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

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

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

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

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

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

f_out.write(out_str)

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


How to convert Bed file to fasta file? I have Bed file containing start and end of a sequence, and I need to convert it to fasta format, any recommendations?

bedtools getfasta can do this. See documentation.


Using cellranger for non-10x data I'm trying to use cellranger to mkfastq and then count and aggregate single-cell data. This data was not generated using 10x Genomics RNASeq.

When I run count, cellranger fails to auto id the chemistry and crashes out. Is it possible to wrangle the data into the cellranger method or should I look elsewhere. I'd ideally like to use the same method for all my single-cell data but not all of it is 10x.

Possible?  Sure.  Your biggest problem is that 10x is expecting barcodes from its whitelist.  I used to alter the python scripts to accept a 'new' kind of chemistry where the barcodes and umi lengths matched what I had, but I found that my modifications stopped working when the version updated,  and I couldn't figure it out, so I started modifying fastqs instead.


How to paste RSIDs in CADD output I want to paste RSIDs in CADD output as CADD does not give RSID column in its output. For this purpose I am using bedtools intersect to compare two files and have RSID column in my CADD file.

This is the command I am using.

bedtools intersect -a CADD.bed -b Input.vcf| cut -f 3

But this command is giving me error

ERROR: file CADD.bed has non positional records, which are only valid for the groupBy tool.

Is there any way to solve this error or any other way to get RSID column in CADD output file.

It would be useful if you show us some line of your CADD.bed and Input.vcf and how your desired output should look like.

bedtools is complaining about non positional records in CADD.bed, which means there are lines in this file where the first three columns aren't CHROM, START, END.

I guess you will have some header lines in the beginning that aren't introduced by an #. Just delete them (or add in # to the beginning of the line) and try again.

fin swimmer


Error in seq.default in chromPlot I am using chromPlot to visualise the genome of C. elegans.

library(chromPlot)

I have created the following data frame with the lengths of C. elegans chromosomes.

Chrom Start      End    Name
1     1     0 15072434 contigs
2     2     0 15279421 contigs
3     3     0 13783801 contigs
4     4     0 17493829 contigs
5     6     0 20924180 contigs
6     7     0 17718942 contigs
7     5     0    13794 contigs

When I use chromPlot(dataframe), I get a plot with 4 chromosomes and the following output/error:

   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
Chrom 1 : 15072434 bp
Chrom 2 : 15279421 bp
Chrom 3 : 13783801 bp
Chrom 4 : 17493829 bp
Error in seq.default(minval, maxval, length = howmany) :
  'from' must be of length 1

Does anybody know why I get this error?

I don't get this error when I use the data frame provided by the package with human segments. The structure of the data frame is:

'data.frame':   457 obs. of  4 variables:
 $ Chrom: chr  "1" "1" "1" "1" ...
 $ Start: int  124535434 121535434 3845268 13219912 17125658 29878082 120697156 120936695 121485434 142731022 ...
 $ End  : int  142535434 124535434 3995268 13319912 17175658 30028082 120747156 121086695 121535434 142781022 ...
 $ Name : chr  "heterochromatin" "centromere" "contig" "contig" ...

Compared to my data frame:

'data.frame':   7 obs. of  4 variables:
 $ Chrom: chr  "1" "2" "3" "4" ...
 $ Start: int  0 0 0 0 0 0 0
 $ End  : int  15072434 15279421 13783801 17493829 20924180 17718942 13794
 $ Name : chr  "contigs" "contigs" "contigs" "contigs" ...

I am on MacOS X.

Update

The 20th line of the package-provided hg_gap data frame is:

   Chrom Start   End     Name
20     1     0 10000 telomere

Note that the Start is 0.

I've also tried changing the start from 0 to 1. The error stays the same.

Update 2

Sorting the data frame by chromosome name does not improve the situation.

The traceback of the error is:

6: stop("'from' must be of length 1")
5: seq.default(minval, maxval, length = howmany)
4: seq(minval, maxval, length = howmany)
3: pretty_ticks(minlab, maxlab, 2, ...)
2: draw.scale(y = 0.8 * chr.length[[chrom]] + 0.2 * xylims[3], minval = margin,
       maxval = maxGeneCount, minlab = annot1_plot_range[1], maxlab = annot1_plot_range[2],
       lwd = 4, col = colAnnot1, cex = cex, title = scale.title)
1: chromPlot::chromPlot(gaps)

I have emailed the maintainer of the package Karen Orostica with a link to this post.

Using this data frame:

gaps
  Chrom Start      End    Name
1     I     0 15072434 contigs
2    II     0 15279421 contigs
3   III     0 13783801 contigs
4    IV     0 17493829 contigs
5     V     0 20924180 contigs
6     X     0 17718942 contigs
7 MtDNA     0    13794 contigs

and this command:

chromPlot(bands = gaps)

I could obtain this plot:

and then, if you want to add some segments for a subset of bases, you can do:

gaps2 <- gaps[1:2, ]
gaps2$End[2] <- 7000000

chromPlot(bands = gaps, segment = gaps2)

and obtain the following plot:

A more satisfactory answer could put segments of different colour on the chromosome, so bands but this is the best I've got so far.


What changed wwith Python in WDL change from draft 2 to version 1? I'm updating my WDL from draft 2 (whichs works fine) to version 1. I'm stuck with a Py snippet
task getReadGroup {

    input {
        String inputFile
    }

    command {echo input = ${inputFile}}

    # command <<<
    # python <<CODE
    # print ("${inputFile}".strip().split('/'))
    # CODE
    # >>>

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

This command gives the expected result:
input = /mnt/lustre/working/genomeinfo/data/20191015_LongRead/20191015_LongRead_1/data/LXBAB201271/LXBAB201271/LXBAB201271/20190417_0612_1-A7-D7_PAD44174_42dcde1a/fast5/test/

But the commented out python code:
print ("${inputFile}".strip().split('/'))

Gives an empty list - "${inputFile}" is empty in the code command for some reason!
Any ideas? I've tried using a File instead of String to make it just like the example in the docs (below) other than that I see no diff:
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("${in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}

Complements of Phil Shapiro;
Use ~{} instead of ${}.
So this example from the current docs is wrong:
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("${in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}

It should be :
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("~{in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}


Public access to genomics databases I'm a statistician, and am interested in applying the theory of algebraic Markov models to genomics. Here's one paper I'm interested in: Large-scale genome-wide association analysis of bipolar disorder identifies a new susceptibility locus near ODZ4. I'm also interested in other diseases (i.e. not just bipolar).

Are there any genomics databases that I can access publicly?

There are some links in Human Connectome Project, but I'm not sure if these are genomics databases (or which papers have used their research).

Most genomics data that is published should be available in GEO (or SRA). Usually in the paper they mention the accession number. However, not every study ends up in there, and sometime studies that are in GEO are not reproducible (poor description of the samples, or design). But if you look at the original papers, most of the times you'll find an accession number.


How to retrieve the subcellular location info from uniprot? I have a list of proteins with their gene names/symbols. I would like to identify which of them are extracellular proteins or plasma membrane proteins.

I can get this information in UniprotKB by searching the proteins one by one. Because I have hundreds of proteins, I wonder whether I can do a systematic search of the UniprotKB in RStudio. I tried the UniProt.ws select method. But I couldn't even get the example code in the package work. I downloaded the uniprot_sprot.dat file. But I don't know how to read this file correctly in RStudio. I used tab-delimited format. It turned out there is only one variable in the file.

Can anybody tell me what I should do to solve the problem? I am pretty new to any computer languages.

Since you are already on UniProtKB, the simplest way will be if you add this information to the output table and download it as a that you can import and filter in R.

To get additional annotations, click the button labelled "Columns" on top of the UniPortKB website.

This allows you to add a wide variety of annotations such as subcellular location to the output. Tick the checkmark for subcellular location and save on the top right.

Then you can download the proteome-wide annotation for example as a tab-separated list using the "Download" dialog next to the "Columns" selection. This list can be imported into R.


What are the optimal parameters for docking a large ligand using Hex? I'm looking to dock a large ligand (~90kDa) to a receptor slightly larger receptor  (~125kDa) using Hex. If anyone is familiar with docking large structures, are there any recommended parameters for finding the best docking solution?

Parameters in particular:

Number of Solutions
N order of correlation for initial and final
searches
Receptor Range
Ligand Range

While I have no experience with this specific question you have, according to Feinstein & Brylinski (2015), a "fully automated procedure", i.e. a Perl script, can be used to optimize the box size itself, and it can be found here.


Negative binomial modeling of RNA-Seq data A common way to model RNA-seq data is using a negative binomial distribution, where each sample-gene pair is modeled by a different negative binomial distribution with mean $\mu_{ij}$ where $i$ and $j$ are indices for genes and samples, respectively (see article). My question is, given we only have a single observation ($X_{ij}$) from each sample-gene pair, how come we can learn the $\mu_{ij}$ from that single observation? We need multiple observations to estimate a mean, right? Maybe I am missing something critical.

The $\mu_{ij}$ in that section of the DESeq manuscript is the expected value of sample $j$ in gene $i$ given its group association (with expected value $\mu_{ij}$). This is why $K_{ij}$ (this is presumably what you meant by $X_{ij}$) is a function of $\mu_{ij}$, rather than $\mu_{ij}$ itself. I think you're just swapping $K_{ij}$ and $\mu_{ij}$ in your head.


Why use Needleman-Wunsch if there is no way to evaluate the statistical significance I thought that Needleman-Wunsch is the best approach to align sequences. However, I read that it is impossible to evaluate the statistical significance of the alignment if you do global alignment. So I would like to ask, what is the point of global alignment if we cannot evaluate statistical significance?

It's best to flip the question on its head and ask yourself instead, "Do we care about the statistical significance of an individual alignment?" The answer is "yes" if you're in a context like homology modeling, where you're aligning two protein sequences. The answer, however, is a very definite "no" in the context of NGS, where you just want the most likely alignment and a vague assessment of its likelihood of being wrong (for filtering purposes).


How to correct alpha, and not p-values themselves, for visualization purposes I have a set of differentially methylated/expressed/whatever entities with p-values attached (example below).

entity_name    p-value    magnitude
entity1        0.04459    0.68
entity2        0.02283    0.99
...
entity_n       0.78       0.025

Typically, I apply the p.adjust function in R with the "fdr" (Benjamini-Hochberg) approach to leave me with p-values adjusted to control the FDR.

adjusted <- p.adjust(mydata,"fdr")

However, I am interested in showing a volcano plot with the unadjusted p-values, and two alpha levels: 0.05 and one that corresponds to the correction. What is the best way to get this alpha? Is it appropriate to set the "corrected alpha" to the lowest original p-value that doesn't pass FDR correction?

The only way to get the alpha levels is to determine what they will be with p.adjust(), since they will depend on the distribution of your unadjusted p values. The general steps you should follow will be:

Add a column of adjusted p-values to your dataframe (mydata$padj = p.adjust(mydata, method="BH"), which is the same as FDR and saves a character).
Use which and max to determine your two alpha threshold (e.g., max(mydata$pvalue[mydata$padj < 0.05])

Then you can adjust your plots however you like (presumably with some horizontal lines at the various alphas). Whether you take the smallest non-significant value or the largest significant value is up to you, just describe what "dots on the line" represent.


Plotting distance tree from blastn output I'm trying to plot a simple distance tree of my blastn output with nj (like the tree view on NCBI). From what I understood, what I think I should do is

extract all the hsps from each alignment
re-align them using MUSCLE
generate a distance matrix from the msa file
plot the tree using nj

However, there is a couple problems I am not sure about:

1) When there are multiple hsps in one alignment then which one do I include in the tree?

2) If two sequences map to two different parts of the query sequence and do not overlap, then what should I do?

Specific questions

When there are multiple hsps in one alignment then which one do I include in the tree?

You can included them all because NJ is not computer intensive, so running loads is fine. You can remove identical sequences if you want either programmatically (these days I use pandas Python) or something like RAxML. Its honestly not critical for an NJ tree (and in fact the more the better) providing you can read the labels. A tree viewing program such as FigTree allows you to collapse clades to make the tree more readable to users.

If two sequences map to two different parts of the query sequence and do not overlap, then what should I do?

Hmmmm.... good question. Do 3 alignments one with one e.g. 5' fragment, the other with the 3' fragment and attempt a joint alignment. Under other non-pairwise tree buidling algorithms this would be disasterous, NJ gets around the problem because it is pairwise. HOWEVER, ask yourself whether your feel the resulting tree is meaningful, I would argue not because there are no homologous sites between the 5' and 3' partial sequences and phylogeny models point mutations. NJ will give you a tree however and will even bootstrap it.

Generic questions, you can actually do this as one of the new features of NCBI's Blast

Go to Blast here, https://blast.ncbi.nlm.nih.gov/Blast.cgi
Enter your sequence into the box (it doesn't accept PDB codes alone)
Enter the protein database - when I first did this calculation I used SwissProt, thinking there would be alot of sequences - I then used "nr"
Under the algorithm parameters enter "50" (default is too many)
Hit "Blast"
Once the search is complete at the top of the page are the hyperlinks: "Other reports:
Search Summary [Taxonomy reports] [Distance tree of results]"
Click on "Distance tree of results"
The tree page will load automatically, automatically aligning your sequences and producing, in this case a parsimony based tree, but there is also the option of a nj tree (recommended) ....
Click "Tool", "Download", "PDF" ...

You can generate a distance and matrix and/or perform a neighbor-joining phylogeny in MEGA10 and using Timura 3-parameter model is a good model. The bootstrap option is important. Its very easy to use, the format is a bit weird

# MEGA
Title: Stuff
# Sequence1
ACTAGACGT
# Sequence2
ACCTTAGGA

etc ...


Appropriate tool or algorithm for sloppy alignment of degenerate bases I have an optimization problem where I have a degenerate nucleotide sequence I want to align against subsets of a reference genome (exons, specifically, to make the problem more tractable).

The sequence itself contains subregions of variable length.

The sequence looks like this, containing degenerate bases:

           L           C           R
5' - SSS......SSS|NNN.....NNN|WWW......WWW - 3'

Per IUPAC notation, S denotes C or G, N denotes any base, and W denotes A or T.

I am marking each block L, C, and R to broadly denote "left", "center", and "right" regions, oriented from 5' to 3' ("forward"-strand oriented).

The widths of the L and R blocks can be from 20 to 25 bases long. The width of the C block can be 14 to 16 bases long.

One programmatic way I could think to do this would be to brute-force a solution, generating all combinations of sequences of varied lengths, blastn-aligning each of such combinations to the sequence for the exon, and then looking for the best alignments among such sequences.

The output would be the scores for regions from sequences that align, to be able to say, from point X to point Y of the genome, a sequence here aligns better than other candidate sequences.

To avoid reinventing tools, and avoid brute-force methods, I'd like to know if there is a tool or even a wrapper around tools, which does this, in case it is already written, or another approach that is more elegant/efficient.

I believe mafft can align degenerate sequences.  From the documentation for version 7 at: https://mafft.cbrc.jp/alignment/software/anysymbol.html

Acceptable nucleotide symbols in the default mode:
a, c, g, t, u
Ambiguous nucleotides (r, y, w, s, k, m, d, v, h, b; IUPAC-IUB codes) can be used and are scored as:
score(r,a) = ( score(a,a) + score(g,a) ) / 2
score(r,g) = ( score(a,g) + score(g,g) ) / 2
score(r,t) = ( score(a,t) + score(g,t) ) / 2
score(r,c) = ( score(a,c) + score(g,c) ) / 2
...
score(r,y) = ( score(r,t) + score(r,c) ) / 2
...
score(r,r) = ( score(a,a) + score(g,g) ) / 2
score(y,y) = ( score(c,c) + score(t,t) ) / 2

MAFFT is designed for multiple sequence alignment, so will attempt to find a local alignment of homologous segments (using FFT), and then attempt to make a global alignment between the input sequences.


how to change the PC use in the dimplot and feature plot I would like to know how to change the PC use in the dimplot and featureplot by using Seurat. for we can get the x-axis and the y-axis like PC-1 and PC-2, if I want to use PC-4 and PC-5. How to modify the code?

You can choose which PC dims to plot by specifying the dims argument within the functions. To use PC-4 and PC-5,

DimPlot(object = pbmc_small, dims = c(4, 5), reduction = "pca")

and

FeaturePlot(object = pbmc_small, features = "Xist", dims = c(4, 5), reduction = "pca")

Seurat documentation is available here.


Is there a straightforward way to get mismatches/indels in a BAM file using pysam? Ideally, I'd like to do something like the following:
f = pysam.AlignmentFile(filename, index_filename=index_filename)
read = list(f.fetch('Chromosome', 1, 10000))[0]
read.get_variations()

And have that output something like:
[
  (1, 'A', 'C'), # An A -> C mismatch at position 1
  (15, 'I', 'CGGG'), # An insertion of the string 'CGGG' starting at position 15 of the read
  (32, 'D', 4), # A deletion of 4 base pairs starting at position 32
]

I know that it's possible to figure out all of those values by parsing the CIGAR string and MD tags but I'm wondering if there's a more intuitive interface.

You're looking for pileup, which is the htslib (and thus samtools/bcftools) method for finding variants. See this section of the pysam documentation. You may want to peruse about the output as well. The gist is:
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            ...handle mismatches...


Biology behind PCA analysis based on SNP What is the biology behind SNP based PCA analysis to study population structure?  I am reading some articles where they used PCA analysis to compare isolates of drosophila that is collected from a different location. AND by PCA they found that SNP is different between all populations. What does that mean? (what it tell you about population structure/genetics?)

First of all, PCA is a technique for dimension reduction. Basically, the goal is to compare tens of thousands of SNPs in Drosophila. Now if you only have 2 SNPs, you can plot them on a 2D scatter plot. If you have 3 SNPs, you may try a 3D plot. But now imagine you have 30,000 SNPs, but you CANNOT plot a 30000-dimensional plot. To visualize this high dimensional data, what we can do is to perform dimensional reduction like PCA. PCA tries to find a set of orthogonal coordinations that explains most of the variation in the data (if there no variation, there is no information contained in the data, which essentially means there is no data). The idea is that PC1 carries most variation can be explained, and PC2 carries the second most... For lower PCs like PC50 or PC60, they probably only carry noise in the data. Therefore, the higher PCs (PC1, PC2 and so on) effectively summarizes the useful information in the data. So you can visualize the "structure" of the data in a 2D PCA plot.
By looking at the distance between points on a PCA plot, you can tell how similar the two data points are. But if you see two populations that are perfectly separated on PCA plot, it does not mean that the 2 population differ completely at every SNP, because PCA is a summarization of all SNP included.