Quality GL entries in vcf file I have a vcf v4.2 file with structural variants, the 8 columns describing the variant are clear, but in the genotype information I got confused by values of GL entries. In the vcf specification they mention that GL are likelihoods of the individual genotypes, but as far as I am concerned, likelihood should be a number between 0 and 1, therefore I am suprised to see negative numbers. Then I thought that it might be log-likelihood, but then it's strange that the higher the number is the better. Here is a sample of 5 variants with various called genotypes

...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-22.4949,0,-4.49485:45:PASS:20:1:11:0:1:4:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-3.1888,0,-43.9892:32:PASS:16:26:13:2:9:2:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-0.190635,0,-49.2888:4:LowQual:26:51:31:2:10:1:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/0:0,-1.01435,-75.899:10:LowQual:18:54:26:2:14:1:0:0
...       GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   1/1:-19.198,-1.5031,0:15:PASS:70:0:58:0:0:5:0:0

The vcf file was generated by delly genotyper. I followed their recommendations for germline SV calling. However, negative likelihoods are even in the vcf specification, therefore it seems that it's not something delly specific.

So, what exactly is GL in vcf files? And how do they relate to other entries GP (genotype posterior probabilities, not in the file by delly) and GQ (conditional genotype quality)?

The VCF spec clearly states that these are indeed log likelihood values:

  • GL : genotype likelihoods comprised of comma separated floating point log10-scaled likelihoods for all possible
  genotypes given the set of alleles defined in the REF and ALT fields. In presence of the GT field the same
  ploidy is expected and the canonical order is used; without GT field, diploidy is assumed. If A is the allele in
  REF and B,C,... are the alleles as ordered in ALT, the ordering of genotypes for the likelihoods is given by:
  F(j/k) = (k*(k+1)/2)+j. In other words, for biallelic sites the ordering is: AA,AB,BB; for triallelic sites the
  ordering is: AA,AB,BB,AC,BC,CC, etc. For example: GT:GL 0/1:-323.03,-99.29,-802.53 (Floats)

With log-likelihood values, the higher the number, the better the value. So -2 is better than -3.


plotMA function issue: can not circle genes I am trying to make MA-plot for my bulk-RNA-seq dataset and experiencing issues with it. I just copy-pasted the code that can be found here:

https://f1000research.com/articles/4-1070/v2

which is:

library("DESeq2")
plotMA(resLFC1, ylim=c(-5,5))
topGene <- rownames(resLFC1)[which.min(resLFC1$padj)]
with(resLFC1[topGene, ], {
  points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})

However, I am not getting the right plot with the top gene circled, but the following one:

Maybe it is a bug in the plotMA function, not sure. How can I reproduce the figure from the publication?

It's likely that the fold-change of the topGene is greater than 5. Consequently, the circle is being drawn and the text produced, but they're outside the bounds of the plot (the bounds don't update if you add a circle and text label). Look at the log2FoldChange column of resLFC1[topGene, ] and increase the bounds to encompass that.


How to filter intervals (reads or genomic coordinates) that have the exact same 5' or 3' ends? Say I have reads that overlap some genes that produce small RNAs, but I want only those reads that start at exactly the TSS of the loci. In other words, reads whose 5' end match the 5' end of a genomic feature.

   5'....3'
...+++++++... Gene
...0000000... R1
...00000000.. R2
..000000000.. R3
..00000000... R4

The output should be R1 and R2.

I believe this is a fairly common bioinformatics operation but somehow it doesn't seem to be an option in the general use tools I looked at (py-bedtools, samtools, bedops).

A solution I thought of would be:

reduce the gene coordinates to its 5' end, and
using bedtools bedClosest to annotate the distance between reads and genes, and finally
select reads that overlap (d=0) and write them to a (bam) file.

This involves quite a bit of wrangling data with pybedtools or awk/bash, and I wonder is there is a more elegant solution to this. Solutions in R/Bioconductor are also welcome.

Presuming you have deepTools installed:

#/!usr/bin/env python
from deeptoolsintervals.parse import GTF
import pysam

anno = GTF(["some/gtf/or/bed/file", "or/more/than/one/if you like/"])
bam = pysam.AlignmentFile("something.bam")
for read in bam:
    s = read.reference_start
    if read.is_reverse:
        s = read.reference_end
    e = s + 1
    overlaps = anno.findOverlaps(read.reference_name s, e, matchType=4)
    # For a GTF, parse the results to see if this is the correct type.
    # For a BED3/BED6 file there will be no exons, so anything other than
    # None indicates an appropriate overlap

    do something

deeptoolsintervals in the GTF/BED parsing and custom interval tree library used by deepTools. anno holds the interval tree and the associated information (gene/transcript name or ID, exons, strand, etc.). the matchType=4 is the pivotal part since that specifies that only overlaps having the exact same start position should be included. Assuming you have a stranded library, you might instead do something like switch between matchType=4 (the same start position) and matchType=5 (the same end position) depending on the orientation of the alignment and the the strand= of the gene.

Note that you'll need to make a set() to hold read names to find the mates (that or jump around the file to get them, which is probably slower).


How to export web NCBI tBLASTn results in table format with many queries? Context

I'm an MSc student working on writing up my thesis (back home now) from my laptop and, therefore, unfortunately don't have access to a workstation/server capable of doing the tBLASTn search that I wanted to do. As a result I have been trying to use as many online resources as possible.

Question

Is there a way to export the tBLASTn results from the NCBI website to a simple text file? I want the output to be something similar to:

FASTA_Name Max_score Total_score Query_cover E_value Identity Accession Gene_ID

Question Details & requirements

Essentially, I have a yeast  database with gene expression information (where the gene names are named in the S. cerevisiae convention eg YMR105C or YDL204W) and I'm trying to relate each of those expression profiles to their respective homologous gene in my species of yeast.

I've exported the systemic file names (and their respective protein residue sequences) for S. cerevisiae, turned that into a big fasta file (~6,000 sequences) and then run a tBLASTn against my yeast genome to give a single Report page where you select each submitted FASTA sequence from the "Results for:" dropdown and it brings up the alignments for that specified sequence as normal. To manually curate this would take forever and I was wondering if there's a way to generate a table containing all the results where the E-value is below a certain threshold?
Ideally it'd be something along the lines of:

"FASTA Name"  "Max score"  "Total score"  "Query cover"  "E value" "Identity"
 "Accession" "Gene*"

*where "Gene" is what you get when you get the title of the page if you were to click that specific "Accession" hyperlink... appreciate that might be beyond the scope of this question though.

Is this possible (or indeed is there a completely different strategy) or should do I need to download the relevant databases and run something along the lines of this?

Thanks for all your help in advance!

It sounds like your easiest option is to just install blast on your local machine and set up the databases there. It really will only take a few minutes to blast 6,000 sentences against a yeast genome on a modern laptop. Any solution interfacing with web results or parsing the output will be painfully convoluted. There are executables for blast that will work on Mac, Linux, or Windows.

Some beginning instructions for installing it on your laptop can be found here: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=Download . Make sure to pass the --num_threads option to speed up the search.


How to convert data in gmt format to dataframe? I downloaded c5: gene ontology gene sets file from http://software.broadinstitute.org/gsea/downloads.jsp

I opened the "c5.all.v6.2.symbols.gmt" file in csv format and It looks like below:

I want to convert the .gmt file into a dataframe. And it should look like below:

can anyone say how I can do that. Thanq

Ok. Answering my own question.

I did this way. And I got the output I want.

install.packages("msigdbr")
library(msigdbr)

m_df = msigdbr(species = "Homo sapiens", category = "C5")
head(m_df)

Output:

gs_name gs_id gs_cat gs_subcat human_gene_symb~ species_name entrez_gene
  <chr>   <chr> <chr>  <chr>     <chr>            <chr>              <int>
1 GO_14_~ M184~ C5     MF        AANAT            Homo sapiens          15
2 GO_14_~ M184~ C5     MF        AKT1             Homo sapiens         207
3 GO_14_~ M184~ C5     MF        ARRB2            Homo sapiens         409
4 GO_14_~ M184~ C5     MF        BAD              Homo sapiens         572
5 GO_14_~ M184~ C5     MF        DAB2IP           Homo sapiens      153090
6 GO_14_~ M184~ C5     MF        DDIT4            Homo sapiens       54541
# ... with 2 more variables: gene_symbol <chr>, sources <chr>

found this solution from here https://cran.r-project.org/web/packages/msigdbr/vignettes/msigdbr-intro.html


Mutations in sequences relative to reference (histogram) I have a reference sequence of amino acids that I compare other sequences to using Bio.pairwise2.align.globalmx of Biopython package (one to many compare - one by one [is there a better way?])

It's works fine but only returns strings of the alignment with very little metadata.

Now I want to extract the indices of the amino acids that are different than the reference sequence and plot it in some kind of histogram according to protein features:

For example:

pos 1-36: 3 mutations.
pos 37-202: 4 mutations.
pos 203: 1 mutation
pos 204-660: 0 mutations

Is there any builtin package to do so?

How can I get the data of diff between 2 sequences in python? I mean which AA changed to which AA and at each position. Bio.pairwise2.align.globalmx returns some data structures of strings (see below) that can be printed but does not contain the locations in a convenient way.

LLGHKVNPFKALYNRIDMK-QVEALVEASKEE
||||||||||||||||||||||||||||||||
LLGHKVNPFKALYNRIDM-RQVEALVEASKEE

You can use AlignBuddy to process alignments in GenBank format, which allows you to keep the metadata. From the command line:

$: alignbuddy input.gb --generate_alignment mafft | \
   alignbuddy --extract_feature_sequences 'PROT_FEAT_1' | \
   alignbuddy --percent_id

If you want to do a large one-to-many comparison and loop through a bunch of features, I'd recommend importing AlignBuddy into a custom python script and calling the relevant functions directly. You can then collect the relevant data and build a histogram with matplotlib.


You really need to learn some python :)
with open('/path/list2.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip())

with open('/path/list1.txt', 'r') as fin:
    with open('/path/list3.txt', 'w') as fout:
        for line in fin:
            if line.strip().split(' ')[0] in to_keep:
                fout.write(line)


What does an FDR value of 1 in RNA-seq mean? I am looking at the supplemental data from the paper "An allelic series of miR-17 ∼ 92-mutant mice uncovers functional specialization and cooperation among members of a microRNA polycistron" which lists the genes that are differentially expressed between a particular knock-out mouse and the wild-type control. In addition to showing the logFC changes for each of the genes, the table includes their FDR value. In many cases, the FDR value is 1. What does this mean?

FDR stands for False Discovery Rate. It is a statistic tool used in multiple hypothesis testing.

As you may know, when you use a p-value cutoff (usually 0.05) for your experiments, it means strictly speaking that "is there was actually no signal, there would be a probability of 0.05 to observe this kind of extreme values". This can be understood as "there's a 5% chance that what you call a hit is actually not signal, just some extreme value that happened by chance, and really there's nothing to see here".

This definition has deep consequences if you run many experiments: imagine that you run 100 experiments (say, you test the effect of 100 different molecules), if none of them has any effect, on average you could still get 5% of false positives. So if you run these 100 experiments, and you get, say, 10 hits with a p-value below your thresold of 0.05; how can you tell which ones are real and which ones are just chance in the absence of biological signal?

This is where you can use FDR.

If you control for a false discovery rate of, say, 0.2 (20%), it means that after a few computations (the Benjamini-Hochberg procedure described on wikipedia is simple and commonly used) you will lower your p-value cutoff so as to be sure that, among the N hits (probably less than the 10 you initially selected above) that you got above, a maximum of 20% are actually false positives (i.e. no biological effect, but surprisingly high/low signal values by pure chance).

In the article you mention they have a different FDR value for each sample, which is unusual. The values are not round numbers, so it seems likely that they are numbers obtained after applying the FDR correction (i.e. they are "transformed p-values") on which they based their decision to include the samples or not. The FDR cutoff that they chose should be written somewhere in the main text.


How to assess if this strategy is good to have some phylogenetic information? I have several organisms and I want to calculate for each one of these their codon composition implementing a pairwise distance matrix between these codon vectors to build a phylogenetic tree. Therefore, the matrix  will have on the principal diagonal  0 values (since each organism will have a distance of zero from itself) and the other matrix members will represent the euclidean distance among one organism and another one. Thus we will have a square matrix.  To assess to robustness of the various phylogenetic relationships I would bootstrap this matrix building the tree on these data.
Is this a good strategy?

The method you are proposing is called a distance method, because you are constructing a distance matrix and then using a clustering algorithm to make a tree. The things to watch out for are:

Multiple substitutions, traditionally this is take care of using a Jukes-Cantor correction. You could argue this is less of a concern because this is a 64 substition matrix model rather than just 4 or 20 substition matrix for nucleotides or amino acids.
Use the neighbour-joining clustering algorithm. Do not use a standard hierarchical clustering algorithm such as UPGMA. I'm sure R must have the nj in something like R's ape, pretty sure scipy will have something like that
Code for bootstrap resampling. That is a method for measuring robustness

The codon model for phylogenetic trees has been implemented in MrBayes, I suspect it will be in phyloBayes. This method will be Bayesian and you should run this is parallel with your own solution.
Concerns
Codon models are very prone to mutation saturation (which is complicated) so you must only use it for closely related sequences. If e.g. 3rd codon distances start exceeding 70% codon models go kaput, because it will select a different codon, which is a completely separate part of the matrix - completely by chance. So whilst mutation saturation in a nucleotide model (4 substitution matrix) would invalidate the tree in theory, they would be sufficient 'good' signal for key phylogenetic signals to be identified, this wouldn't happen for a codon model.
I recognise the 'concerns' might sound a bit confusing at present, but just to reiterate, codon models need sequence data that is really closely related. For this data they excel and give excellent bootstrap support.


Batch Effects in RNA Seq Sample I am using the R (using EdgeR) for the RNA Seq analysis, I had few batch effect samples like Control vs treatment. Could anyone tell me the best way to remove the batch effects.

I have looked into the cluster dendrogram and MDS post (using EdgeR), the Control and Treatments are not grouping together, and their is differences in the different batch within Control and Treatment.

I tried this syntax to remove the batch effects for the rest of the analysis to get the DE Genelist. I feel it didn't able to remove the batch effects from the rest of the analysis.

# define the experiment
Batch <- factor(c(1,2,3,1,2,3))
Treatment <- factor(c(0,0,0,1,1,1))
design_noI <- model.matrix(~ 0 + Batch + Treatment)

# define our DGEList
San <- DGEList(counts=x,genes=gene_info)

A <- aveLogCPM(San)
y2 <- San[A>1,]
y2 <- calcNormFactors(y2)
logCPM <- cpm(y2, log=TRUE, prior.count=5)
plotMDS(logCPM)
logCPMR <- removeBatchEffect(y2, Clutch = Clutch)
plotMDS(logCPMR)

I can see the difference in MDS plot on removing the batch effects but I couldn't able to carry the same data flow into my rest of analysis.

# filter the lowly expresed genes and normalize the data
keep <- rowSums(cpm(San/logCPMR)>1) >= 3
keep2 <- removeBatchEffect(keep, Batch)
y_group <- San/logCPMR[keep2, , keep2.lib.size=FALSE]
y_group <- calcNormFactors(San/logCPMR)

I have mention the San/logCPMR, as San take as data without remove the batch effects, whereas I think logCPMR as batch effected removed data frame.

Any suggestion on this more useful to me.

Please let me know, any other information needed to deal with this batch effects.

The help page for removeBatchEffect() explicitly states that it should not be used in conjunction with differential expression testing, so don't do that. You already have the batch effect included in your model, so you're already correcting for it when testing for differential expression. If your MDS plots suggest that samples are variably affected by the batch effect then look into the combat() function from the sva package.


Extract LINE-1 methylation status from WGBS data I've been trying to figure out the correct way of extracting methylation status of particular elements from WGBS data, and found several possible approaches, but would like to discuss which one could be the best:

using software dedicated to the extraction of these elements on WGS (like TraFiC or MELT and then quantify methylation (like Bismark)
using bedtools to scan bam file with a filtered version of UCSC RepeatMasker table in bed format) and then quantify methylation
I've seen the use of custom perl scripts on fastq files like in this paper, in which the Supplemental Methods Section states: "A custom Perl script was written to assay methylation of CpGs within a previously published motif of LINE-1 elements (Yang et al., 2004). Methylation levels for each of the four CpGs within the consensus LINE-1 element were determined globally for each WGBS sample prior to alignment using Line1_FASTQ.pl. Dup15q and control were compared and tested for significance using an unpaired t-test."

Any thoughts on this?

After a more extensive search, I found a useful answer that I'll share in case anyone else wonders how to do this. I'm using MethylKit package for differential methylation expression and it is quite easy to integrate with Genomic Ranges and rtracklayer in R to extract particular element methylation differences:

element.of.interest <- import.bed("FILENAME.bed") #download custom track from Genome Browser
element.meth <- diffmeth[as(diffmeth, "GRanges") %over% element.of.interest]

where "diffmeth" is the resulting MethylKit object from the differential methylation analysis.


median normalization for proteomics I am using the data from a proteomics study were the data was log2 transformed and then a median normalization was applied. The data was normalized by groups of conditions (normal, mutant), not for all the arrays at the same time, so the median for the normal group is slightly higher than for the mutant group.
I am wondering if this is correct or should I normalized the data en general, I mean taking into account all the arrays and normalizing all at the same time so the 2 groups have the same median.
Regards

The kind of normalization depends on what you what to explore. So, there is no absolute answer to this. Different normalizations highlight different aspects of your dataset. There is a nice paper for a higher throughput quantitative study that I would recommend reading:
Nusinow, David P., John Szpyt, Mahmoud Ghandi, Christopher M. Rose, E. Robert McDonald 3rd, Marian Kalocsay, Judit Jané-Valbuena, et al. 2020. “Quantitative Proteomics of the Cancer Cell Line Encyclopedia.” Cell 180 (2): 387–402.e16.
What you want to aim for is a normalization that removes batch effects and that clusters together your quality control samples that you know are very similar to each other.


Analyzing proteins based on sequence similarity Using NCBI Gene database, I retrieved list of 7015 protein which are implicated in autism pathogenesis. This list includes tax-ID, organism name, symbol, aliases, map location, chromosome, genomic nucleotide accession version, start position on genomic accession, end position on genomic accession, orientation, exon count and OMIM). Among these 7015, 1144 proteins are involved in neurodevelopmental functions. I want to know that how many of these 1144 proteins are functionally similar to 7015 proteins? Is there a method to query 1144 proteins against 7015 proteins and group them based on structural similarity/ function similarity/ gene family?. Since I have large set of proteins, I am trying to understand if sub-set of proteins share similar function with this large set of proteins.

Thanks

Priya

I suggest you give DAVID a try. Specifically, their Functional Annotation tool. Just enter your list of protein IDs, and it will return groups of proteins where particular GO functions are overrepresented.

This should at least help you start categorizing your list into functional groups.


Searching for start and stop codons for protein sequencing of contigs I need to convert contigs into their respective protein sequences given a reference genome (i.e. I need to take a substring, whose position is already known on the string, and I need to locate the nearest start and stop codons).

This is tricky because sometimes the first position of triplet codon within the contig will not be a multiple of three and could be located in the intergenic region (i.e. between genes). The goal is to separate both coding and non-coding DNA.

This is my code so far:

from Bio.Seq import Seq
from Bio.Alphabet import generic_dna, generic_protein

start_codons = ['ATG']
stop_codons = ['TAG', 'TAA', 'TGA']

string = 'GG*TAG*CCAATT*ATG*AACGAA*TAG*GAC' #remove '*', just for visual
contigs = ['CCAA', 'TGAAC', 'GAA', 'GGAC']
positions = [5, 12, 17, 22] #position indices for each contig on string

extended_contigs = []
extended_position_contigs = []
intergenic_contigs = []
intergenic_position_contigs = []
for i in contigs:
    extended_contigs.append(#some code)
    extended_positions_contigs.append(#some code)
    intergenic_contigs.append(#some code)
    intergenic_positions_contigs.append(#some code)

I should get extended_contigs = ['ATGAAC', 'ATGAACGAA'] and extended_positions_contigs = [12, 17].  These are the contigs that are located within genes. In order to code them into peptides, I need to reach back in the string until I find the start codon and extend the initial contig (e.g. TGAAC -> ATGAAC and GAA -> ATGAACGAA)

I should also get intergenic_contigs = ['CCAA', 'GGAC'] and intergenic_positions_contigs = [5, 22].  When the missing code is run, the computer searches to the left of the string and finds a stop codon (e.g. TAG) before the start codon. Thus, the contig is located in between two genes and nothing needs to be added. These intergenic contigs are just stored in a new list.

My code continues:

prot_contigs = []
for i in extended_contigs:
    my_dna = Seq(i, generic_dna)
    my_prot = my_dna.translate()
    prot_contigs.append(str(my_prot))

Here, no new code needs to be added. After the above is run, prot_contigs = ['MN', 'MNE'].

The last step of code (which I need help with) will transform prot_contigs into new_prot_contigs = ['MN', 'E'].

How? If for any contig (e.g. 'TGAAC'), the beginning or end is part of another codon (not perfect multiple of 3), the extra codon on either end will be maintained (e.g. 'MN' stays 'MN').  Otherwise, if the contig (e.g. 'GAA') perfectly matches a codon, anything that was added onto it will be removed (e.g. 'MNE' becomes 'E').

I would try to solve the 2 portions of code myself, but I am unsure on how to take a position on a string (i.e. beginning points of contig) and look along the string to find the nearest start/stop codons, so I could determine the DNA's function and accurately sequence the protein coding contigs into peptides.

Any help would be appreciated!

If your string is 'GGTAGCCAATTATGAACGAATAGGAC' then I would assume that you could separate this into coding and non-coding regions such as:
GG*TAG*CCAATT*    (non-coding)
ATG*AACGAA*TAG    (coding region)
*GAC              (non-coding)

It starts with a start codon ATG and ends with an end codon TAG.
Your first event position 12 does not make sense to me. What does * mean? It seems of that your sequence starts with GG which is not a triplet? Perhaps did you mean GGG?
I'm the author of node based DNA visualisation AminoSee. It has a bug in that it only renders half right: because DNA can be read backwards, and my program only does it forwards. I actually came to BioInformatics today to try and get my head around reading backwards and how that might work.
AminoSee only transcodes 3 base pairs at a time, and discards triplets that aren't, well, 3 base pairs longs. I'm self taught and so I guess unaware if it's possible to have any sequence that is not a multiple of 3. Perhaps it's possible. I assumed not. My program is more of an art piece that a science tool, but I plan to use it one day to perform a "visual diff" between two species (Humans and Chimp). See AminoSee  DNA Visualisation over at aminosee.funk.nz


How to get the corrected matrix after SVA batch effect correction I ran SVA to remove batch effects for my bulk RNAseq experiments, but now I need to somehow correct my data matrix in order to run pca, mds. I am using DESeq2 for the analysis. Here is the code that I got now:

dds <- estimateSizeFactors(dds)
dat = counts(dds, normalized = TRUE)
idx = rowMeans(dat) > 1
dat = dat[idx,]
mod = model.matrix(~Group, colData(dds))
mod0 = model.matrix(~1, colData(dds))

nsv = num.sv(dat, mod)
svseq = svaseq(dat, mod, mod0, n.sv = nsv)

dds_sva = dds
dds_sva$SV1 = svseq$sv[,1]
design(dds_sva) = ~ SV1 + Group
dds_sva = DESeq(dds_sva)

So, I am wondering now, how to use svseq object to plot, say, the corrected version of MDS plot.

If you want to plot the "corrected" expression, you will need to remove the variation introduced by these surrogate variables. Removing the expression affected can introduce some bias too and it is usually not recommended (despite comBat doing so). You should apply linear algebra, you can look at here is an example how to do it:

library("corpcor")
s <- corpcor::fast.svd(t(scale(t(dat), center = TRUE, scale = TRUE)))
pcSds <- s$d
pcSds[1] <- 0
svdexp <- s$u %*% diag(pcSds) %*% t(s$v)
colnames(svdexp) <- colnames(dat)

But if your surrogate variables do have a big enough influence on your MDS plot to be observable (IMHO) there would be too much batch effect to work with the data.


How to highlight specific cells in Seurat 2.4 I used Seurat 2.4 on our scRNA dataset to obtain the following tSNE plot.I was able to successfully extract cell IDs from the different clusters, and generate gene expression profiles. The analysis, and the biology makes sense. I have been however stuck in trying to highlight specific cells we are interested in using the Cell IDs (barcodes). For e.g, in this plot I am hoping to highlight specific cells (~ 120 of them) within Cluster 5 (blue cluster on the top). It can be either in featureplot mode or in this plot itself by an overlay, it doesn't matter. All I have to show are the 120 cells within the cluster.
For eg. if cluster 5 has cells atag-1 atgc-2 atat-3, cacc-4 cat-5... i want to recreate this plot using atag-1 atgc-2, cacc-4

I tried using do.identify = TRUE but that is a visual selection and can select regions, but is not useful to select specific cells. I was thinking if there is a way to subset and replace cluster 5 to include just the cells (~120 odd out of the total 200) but was unable to crack it. I tried by using cells.use but then realized that doesn't work either.Any ideas will be greatly appreciated, to get me out of this rut.
Thank you

You can "highlight" any cells you'd like by creating a new column (or modifying an existing column) in the meta data slot within the Seurat object. From there you can modify the contents of that vector to your liking and plot the results.

It's not clear if you've been successful in obtaining the barcodes for the subpopulation to begin with, because that's a separate issue all together, and you need to clarify what attribute you are trying to subset the cells on because there are a number of ways to accomplish that. If you've already obtained the barcodes then I'll describe how you would plot the results below.

Let's assume you'd like to keep the original cluster IDs for all other cells. We'll start off by duplicating the cluster ID column (this isn't necessary, but it's good practice to preserve the original IDs).

n <- seurat_obj@meta.data # for convenience
n$new_column <- n$res.1.0

Now you need to modify the ids for cluster 5 to break them up into 2 populations.

n[rownames(n) %in% my_barcodes, "new_column"] <- "5b"

You can plot the new column using DimPlot with the argument group.by = "new_column", which will make the subpopulation stand out from the rest of the cells in cluster 5.


GWAS, MWAS, EWAS: what are the (in)dependent variables? I started reading some papers on X-wide association studies, where X can be metabolome, epigenome, etc... The authors usually describe which are the dependent and which are the independent variables of their models. Taking as an example a metabolome-wide association study (MWAS), the authors use as model the following: metabolite ~ condition. That is, each metabolite is considered the dependent variable and the condition (e.g., smoker vs. non-smoker) is considered the independent variable. I was a little bit confused, so I looked into an example study using GWAS: in this case the model can be written as condition ~ SNPs, which is what I would expect (knowing nothing about GWAS but something about Data Science) and quite the opposite from the MWAS described above.

Does the choice of the dependent variable depend on the questions asked in the study, or does it exist some sort of agreement? Moreover, regarding the GWAS model described above, I read that usually researchers fit a model for each SNP: what's wrong with fitting a model which consists of several independent variables (i.e., y ~ SNP_1 + SNP_2 + ... + SNP_N)?

To answer the first part of your question, the dependent and independent variable of X-WAS is kind of arbitrary and dependent on the question you asked. But it gradually becomes a convention in the field after the initial name and concept are accepted by the community. For example, GWAS from the very beginning is written as condition ~ SNP, and there is no doubt about it. Later, people did something like condition ~ gene expression and called it TWAS, which is also widely accepted now. In your case, MWAS is indeed metabolite ~ condition, and the only important thing is to know the convention in the research community so people are on the same page when they talk to each other.

As for why running GWAS using single SNP. In a GWAS, it is typical to use 500K to 7M SNPs. So this is a large number. It is obvious that we cannot include all the SNP in one single model, then the question becomes how many SNPs to include in each model, and how to choose those SNPs. This is very hard to decide. Also, due to link disequilibrium, many SNPs are correlated, which means co-linearity, a thing you want to avoid in linear models. It is not clear how to correctly interpret the effect sizes in a model with many SNPs as well. Therefore, due to feasibility, simplicity, and interpretability, it might be best to perform regression on every single SNP.

Now since we are running GWAS in an SNP-by-SNP fashion, there are a lot of work that needs to be done after GWAS in order to find the true causal SNP. Many fine-mapping methods have been developed to tackle this problem.

Interestingly, there are some multi-variate GWASs methods, but they do not come in the flavor of simple linear models. Instead, sometimes complex Bayesian statistics are used, for instance.


KeyError when getting features from a genbank file with biopython with some accessions but not others I'm very new to python but i've been using it to extract the sequence of a gene from a genbank file. The issue is is that sometimes i'll get the output i want (prints the sequence to a file) and sometimes it will return a key error. This depends on which accession i'm using. Does anyone know why it sometimes might give a key error? I thought it might be something to do with the genbank records themselves, but they look pretty similar and the gene is there (in the gene feature qualifier). EG works with HG738867.1 but not AP019703.1. Here's my code -
from Bio import Entrez, SeqIO

gi_genome = 'accession'
name = 'acrA'
Entrez.email = 'email'
handle = Entrez.efetch(db="nucleotide", id=gi_genome, rettype="gbwithparts", retmode="text")
record = SeqIO.read(handle, "gb")
handle.close()
element = 0
for feature in record.features:
    if feature.type == 'CDS' and name in feature.qualifiers["gene"]:
        report = 'record.features[%s]' % str(element)
        gene_sequence = feature.extract(record.seq)
        with open('output.fasta', 'a') as f:
            print('>' + gi_genome + ' ' + name, file=f)
            print(gene_sequence, file=f)
        break
    else:
        element = element + 1

Here's the traceback -
Traceback (most recent call last):
  File "/home/ubuntu/Documents/Git_Branches/Project_planning/Learning/In_progress/utils/data.py", line 11, in <module>
    if feature.type == 'CDS' and name in feature.qualifiers["gene"]:
KeyError: 'gene'

Process finished with exit code 1

Thanks in advance!

Good question.
feature.qualifiers is a dict, if the dict doesn't have that key it will throw a KeyError. The way your code works it expects feature.qualifiers to have 'gene' as a key every time feature.type == 'CDS'. However, it is possible to have feature.type == 'CDS' without a gene entry and in the case of AP019703.1 this comes up before you reach the gene acrA.
As simeple fix is to simply use get instead:
feature.qualifiers.get("gene",[])

Now returns a empty list if there is not key "gene", which also satisfies the 'name in' part of your if statement.
If you get stuck on stuff like this use more print statements, break up your if statements and look at the whole file (example below) to understand the structure.
gi_genome = 'AP019703.1'
name = 'acrA'
Entrez.email = 'email'
handle = Entrez.efetch(db="nucleotide", id=gi_genome, rettype="gbwithparts", retmode="text")
record = SeqIO.read(handle, "gb")
handle.close()
with open('AP019703.gb' ,'w') as fout:
    SeqIO.write(record, fout, 'genbank')


Is there a way to quickly verify the presence of some SNPs in Fastq files? Is there a tool that can scan fastq files without assembling them for a custom list of user defined snps?

Try this pipeline:

Clean raw fastq files (trimmomatic)
Align clean fastq files to the reference genome (dna-to-dna aligner of your choice)
Convert the alignment sam file to bam (samtools)
Call SNPs from bam file (samtools)

Here's a start tutorial on this: http://ged.msu.edu/angus/tutorials-2013/snp_tutorial.html


Extracting strand specific reads from MinION cDNA-PCR protocol I recently performed my first MinION run, and now I'm trying to analyze the data. Being pretty new to the bioinformatics field, I was hoping some of you could help me out.

As a bit of background, I'm trying to define transcript isoforms, and down the line, I'm hoping to quantify isoform abundances. As my first attempt, I used a 1D kit that first performs RT starting at the polyA tail and then template switching to enrich for full length transcripts. To initiate RT, a poly-dT primer is added and after template switching a second primer is added. Using those added primer sequences, I then performed PCR to amplify the full-length cDNAs. This seemed to work well, so I went ahead and sequenced the library.

Now, on to my problem. I have aligned my fastq file to my genome with both bwa-mem and minimap2. In both cases, I get alignment of reads to the genome where I expect transcription to be occurring, but half of the reads are in the sense direction and half are in the antisense direction. This is presumably due to the fact that the cDNA is double stranded and either strand can be recruited to pore to be sequenced. However, because I'm interested in mapping transcript isoforms which could include anti-sense transcripts, this isn't exactly what I hoped for.

I'd like to figure out a way to extract only the reads that correspond to the actual mRNA transcript's sequence and not the reverse complement. Does anyone have ideas as far as programs that would help with this?
I'm thinking that I should be able to extract out only the molecules that have the polyA or have the strand switching primer sequence or something like that, but I'm not sure how to properly implement this.

This has been set up as a bioinformatics protocol on protocols.io.

I've been using LAST to identify the adapter orientation relative to the genome, and then using that information to split BAM files up and recombine them to create two strand-specific files that are displayable in a genome browser.

I have written my own script to process LAST results into a CSV format, which makes it easier to do line-by-line data filtering.

Note: these scripts have been slightly modified from scripts that I have run. Consider them demonstrative: pay attention to the words rather than the script.

Read correction

I prefer starting off my data analysis with a read correction with Canu (ideally v1.7 now that it's out, because that attempts correction of all reads, but here I use canu v1.6). I use minimap as the mapper to speed this up. The genomeSize parameter should be approximately a tenth to a fortieth of the number of bases in your dataset to make sure that no sequences are excluded (bigger is better, as long as Canu doesn't freak out about memory consumption):

~/install/canu/canu-1.6/Linux-amd64/bin/canu overlapper=minimap \
  genomeSize=100M minReadLength=100 minOverlapLength=30 -correct \
  -p 4T1_BC06 -d 4T1_BC06 -nanopore-raw \
  workspace/pass/barcode06/fastq_runid_*.fastq

This creates a file 4T1_BC06/4T1_BC06.correctedReads.fasta.gz.

Chimeric read filtering

The next step I carry out is a basic read-level QC to exclude chimeric reads. Porechop can be used for this, although that removes adapters by default, which is not particularly useful in this case.

I use LAST to search for adapter sequences within the corrected reads, pass it through my conversion script, and extract out duplicated mappings (i.e. where the same read/adapter pair appears more than once in the mapping results):

lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | awk -F',' '{print $1,$2}' | sort | \
  uniq -d | awk '{print $1}' | uniq > reads_with_duplicated_adapters.txt

I use this file to filter out chimeric reads from the corrected dataset using another fastq filtering script I've created:

~/scripts/fastx-fetch.pl -v -i reads_with_duplicated_adapters.txt 4T1_BC06/4T1_BC06.correctedReads.fasta.gz | \
  gzip > 4T1_BC06.correctedReads.uniqueOnly.fasta.gz

Adapter filtering

Along roughly similar lines to the chimeric read filtering, I then look for the strand switch and VNT adapters in the sequences. For a forward-orientated query, I expect the strand switch primer to be in the forward direction, and the VNP primer to be in the reverse direction:

lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,+' -e 'ONT_VNP,-' | \
  awk -F',' '{print $1}' | sort | uniq > fwdQry_seqs_BC06.txt
lastal -P 10 barcodes_primerSeqs.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,-' -e 'ONT_VNP,+' | \
  awk -F',' '{print $1}' | sort | uniq > revQry_seqs_BC06.txt

And then more filtering to split the reads up into forward and reverse [genome-relative] subsets:

~/scripts/fastx-fetch.pl -i fwdQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz
~/scripts/fastx-fetch.pl -i revQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz

Mapping

Now that the reads have been oriented, the mapping can be done:

~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > fwd_4T1_BC06.CU_vs_mmus.bam
~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > rev_4T1_BC06.CU_vs_mmus.bam

Splitting and recombining

The mapped files are then split based on their mapping direction using samtools view flag filtering to exclude or include reverse-mapped reads:

samtools view -F 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > fwd_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > rev_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -F 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > rev_rev_4T1_BC06.CU_vs_mmus.bam

And finally the files are recombined to identify forward and reverse-encoded transcripts via samtools merge. I'm calling the "pos" strand the one that is encoded in the same direction as the genome, and the "neg" strand the reverse-complement direction:

samtools merge pos_4T1_BC06.CU_vs_mmus.bam \
  rev_fwd_4T1_BC06.CU_vs_mmus.bam fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools merge neg_4T1_BC06.CU_vs_mmus.bam \
  fwd_fwd_4T1_BC06.CU_vs_mmus.bam rev_rev_4T1_BC06.CU_vs_mmus.bam

Alternatively, it's possible to set the BAM flags for read 1 and read 2 for the reverse and forward-encoded adapters respectively, which allows everything to then be combined into a single BAM file (and treated in the same way as a strand-specific Illumina BAM file).

If this has worked properly, then mapping human or mouse to the mitochondrial genome should show most expression appearing on the positive strand, with a small scattering of negative-strand expression, a bit like this:

For a bit of a background on how/why this works, see this answer.


Count genomic ranges I have a set of genomic ranges that are potentially overlapping. I want to count the amount of ranges at certain positions using R.

I'm Pretty sure there are good solutions, but I seem to be unable to find them.

Solutions like cut or findIntervals don't achieve what I want as they only count on one vector or accumulate by all values <= break.

Also countMatches {GenomicRanges} doesn't seem to cover it.

Probably one could use Bedtools, but I don't want to leave R.

I could only come up with a hilariously slow solution

# generate test data
testdata <- data.frame(chrom = rep(seq(1,10),10),
                       starts = abs(rnorm(100, mean = 1, sd = 1)) * 1000,
                       ends = abs(rnorm(100, mean = 2, sd = 1)) * 2000)

# make sure that all end coordinates are bigger than start
# this is a requirement of the original data
testdata <- testdata[testdata$ends - testdata$starts > 0,]

# count overlapping ranges on certain positions
count.data <- lapply(unique(testdata$chrom), function(chromosome){
    tmp.inner <- lapply(seq(1,10000, by = 120), function(i){
        sum(testdata$chrom == chromosome & testdata$starts <= i & testdata$ends >= i)
    })
    return(unlist(tmp.inner))
})

# generate a data.frame containing all data
df.count.data <- ldply(count.data, rbind)

# ideally the chromosome will be columns and not rows
t(df.count.data)

GenomicRanges::countOverlaps seems to be what you’re after:

position_range = GRanges(position$chrom, IRanges(position, position, width = 1))
ranges_at_position = countOverlaps(position_ranges, granges)


Convert a BAM file from one reference to another? I have a set of BAM files that are aligned using the NCBI GRCh37 human genome reference (with the chromosome names as NC_000001.10) but I want to analyze it using a BED file that has the UCSC hg19 chromosome names (e.g. chr1). I want to use bedtools to pull out all the on-target and off-target reads.

Are NCBI and UCSC directly comparable? Or do I need to re-align the BAM/lift-over the BED to the UCSC reference?
Should I convert the BED file or the BAM file? Everyone here uses the UCSC chromosome names/positions so I'll need to convert the eventual files to UCSC anyway.

You're the second person I have ever seen using NCBI "chromosome names" (they're more like supercontig IDs). Normally I would point you to a resource providing mappings between chromosome names, but since no one has added NCBI names (yet, maybe I'll add them now) you're currently out of luck there.

Anyway, the quickest way to do what you want is to samtools view -H foo.bam > header to get the BAM header and then change each NCBI "chromosome name" to its corresponding UCSC chromosome name. DO NOT REORDER THE LINES! You can then use samtools reheader and be done.

Why, you might ask, would this work? The answer is that chromosome/contig names in BAM files aren't stored in each alignment. Rather, the names are stored in a list in the header and each alignment just contains the integer index into that list (read group IDs are similar, for what it's worth). This also leads to the warning above against reordering entries, since that's a VERY convenient way to start swapping alignments between chromosomes.

As an aside, you'd be well served switching to Gencode or Ensembl chromosome names, they're rather more coherent than the something_random mess that's present in hg19 from UCSC.

Update: Because I'm nice, here is the conversion between NCBI and UCSC. Note that if you have any alignments to patches that there is simply no UCSC equivalent. One of the many reasons not to use UCSC (avoid their annotations too).


Intuitive explanation of GSVA analysis I'm trying to understand the way the GSVA analysis is working behind the scenes.And I was wondering if there is any way to understand it more intuitively the whole process.

So at first according to paper it starts by evaluating whether a gene i is highly or lowly expressed in sample j in the context of the sample population distribution. They use these kernel estimations of the cumulative density functions to transform the initial values so not to be affected by the problematic intensities.

After this "transformation" and a following normalization, GSVA calculates the enrichment scores using the Kolmogorov-Smirnov (KS) like random walk statistic.

As I know, Kolmogorov-Smirnov checks for differences in distributions. Which distributions does it check? Gene-set's against all the others genes? And what is the role of the random walk?

So is there any intuitive way to understand this kind of Kolmogorov-Smirnov (KS) like random walk statistic? How does it actually work? Which one is the null and which the alternative hypothesis in that case?

The most intuitive explanation is  also explained in the background section:

  Conceptually, this methodology can be understood as a change in coordinate systems for gene expression data, from genes to gene sets.

Which I think it is explained in Figure 1 of the paper, which is also one of the most informative I found about the methods of an algorithm1.

I'll will use it to explain what's going on.

Fit a distribution function per gene

  For each gene expression profile ... a non-parametric kernel estimation ... function is performed.

I can't give more detail because I haven't understood it fully

Then they are normalized to make the ranks symmetric around zero

It is an inline formula between equations 2 and 3.
The Kolmogorov-Smirnov like random walk statistic is applied to those kernel estimations normalized distributions.

Which two distributions does it compare? The genes in a gene set and all the others (but I'm not totally sure)

The role of random walk is explained a bit later from your quote:

  [the random walk] produces a distribution over the genes to assess if the genes in the gene set are more likely to be found at either tail of the rank distribution

Which is later used to calculate the enrichment score by looking the maximum deviation from zero, or by summing the largest and the lowest deviations.

To address some other questions:

  GSVA calculates sample-wise gene set enrichment scores as a function of genes inside and outside the gene set, analogously to a competitive gene set test

1: Try finding which fittings and approximations are done inside limma


How to run a one-time 'setup' step for amrfinder in snakemake using conda I'm trying to add amrfinder running within conda to a snakemake workflow, working on a slurm HPC.

The first time you run amrfinder, you have to run amrfinder -u to update/download the database of genes/proteins which it searches for. Since I'm using the --use-conda flag of snakemake, the first time I run the workflow, it creates a new env for the conda package, so this command has to be run before that copy of amrfinder will work.

My question is: how do I run amrfinder -u within snakemake? If I do it as part of the standard amrfinder command, it runs as many times as there are samples, and as all the jobs are looking at the same database, when it updates for one job, it breaks the other jobs.

I can't figure out how to run this command in a 'snakemake-y' way, as it has no input, it's just literally amrfinder -u.

Of course, I can make it work with a hack/workaround, but I'm interested in a snakemake way of doing it.

This is the main snakemake rule for running amr_finder_plus, the new rule, to run amrfinder -u needs to run before this

rule amr_finder_plus:
    input:
        assembly = 'my_assembly.fasta'
    output:
        amr_finder_plus_results = 'my_amr_finder_plus_results.tsv'
    conda:
        '../../envs/amrfinderplus.yaml'
    shell:
        'amrfinder -n {input.assembly} -O Salmonella --output {output.amr_finder_plus_results} --threads 4'

There are two approaches you can take:

There's a --conda-create-envs-only option that you can use to exit after creating the envs. You should get information about where the envs are so you can run amrfinder -u.
Add a new rule called "amr_setup" or something like that that runs:

.

amrfinder -u
touch amrfinder.setup

Then put amrfinder.setup into your snakemake workflow as an input file into the amr_finder_plus rule. Either of these will work


biological fluorescence supervised machine learing Some general questions...

So it appears with most classification learners one must come up with a series of quantifiable variables to associate with each observation.  This might include mean and std deviation intensity values, texture, and zernike polynomials.  Does anyone have any other suggestions as to how one might quantify 3-dimensional fluorescence pixel data features?

Once appropriate quantification's have been acquired for a given class one can then train a classifier.  My next question is - does one then have to segment test images and create a matrix of the same feature variables from that segmented image to use the classifier? Or can one use the raw image itself?

What supervised machine learning techniques will evaluate the raw image data itself - either for training or testing? Can manually annotated data sets or clustered data sets provide a direct model for training data?  That is to say, are there supervised machine learning algorithms which accept pure pixel data as models for training?

Are there reliable pre-trained classifiers that are already implemented for these type of fluorescence data sets?

Any resources or advice is greatly appreciated.

Thanks in advance for any input!

Whether to quantify things at all with depend on the technique you want to use. An end-to-end CNN (convolutional neural network) for example would just be fed the 3D image for training and predictions. However the major downside to CNNs is that you need a LOT of data to train them, though perhaps you can get lucky and find an existing CNN that you can use and just retrain the last couple layers.

Other than that, you're going to need to come up with as many (hopefully useful) metrics as you can. There are then a lot of techniques one could try, ranging from simple things like regression to random forests to SVMs, among many others. In general I think simple regression and random forests are the easiest to get started with.


Detecting broad peaks in sRNA-seq data What kind of tool would be appropriate do detect "broad peaks" in small RNA-seq sequencing data?

MACS2 appears to be  developed for ChIP-seq data, but I see that there is a --nomodel option. Would that make this program usable in my case?

I suppose I should also set the --shift and --extsize options, but I'm not sure how:

  --nomodel             Whether or not to build the shifting model. If True,
                        MACS will not build model. by default it means
                        shifting size = 100, try to set extsize to change it.
                        DEFAULT: False
  --shift SHIFT         (NOT the legacy --shiftsize option!) The arbitrary
                        shift in bp. Use discretion while setting it other
                        than default value. When NOMODEL is set, MACS will use
                        this value to move cutting ends (5') towards 5'->3'
                        direction then apply EXTSIZE to extend them to
                        fragments. When this value is negative, ends will be
                        moved toward 3'->5' direction. Recommended to keep it
                        as default 0 for ChIP-Seq datasets, or -1 * half of
                        EXTSIZE together with EXTSIZE option for detecting
                        enriched cutting loci such as certain DNAseI-Seq
                        datasets. Note, you can't set values other than 0 if
                        format is BAMPE for paired-end data. DEFAULT: 0.
  --extsize EXTSIZE     The arbitrary extension size in bp. When nomodel is
                        true, MACS will use this value as fragment size to
                        extend each read towards 3' end, then pile them up.
                        It's exactly twice the number of obsolete SHIFTSIZE.
                        In previous language, each read is moved 5'->3'
                        direction to middle of fragment by 1/2 d, then
                        extended to both direction with 1/2 d. This is
                        equivalent to say each read is extended towards 5'->3'
                        into a d size fragment. DEFAULT: 200. EXTSIZE and
                        SHIFT can be combined when necessary. Check SHIFT
                        option

In my case the tags are the same things as the fragments, so this "extension" thing is probably not appropriate.

Basically, I would like to call peaks based on the plain real coverage of my small reads. Are there other tools that would be more appropriate?

The sure why any software should be confined to only broad peaks, but there is two packages that should do what I thing you have in mind. Just to clarify, for smallRNAs peak detection is sometimes also referred to as loci or cluster detection:

ShortStack which is developed for alignment, annotation, and quantification of small RNAs. The de-novo cluster finding option as copied from the manual might do the trick:

DE-NOVO CLUSTER FINDING
    Unless options --locus or --locifile are used (see below), ShortStack
    will de-novo identify clusters of small RNA accumulation genome-wide.
    Cluster definition is simple: First, all regions containing at least one
    primary alignment are found where the maximum distance between the ends
    of the alignments is <= option --pad (default: 75). Second, if the
    number of alignments in the cluster is >= option --mincov (default:
    0.5rpm), the cluster is kept. The mincov threshold can also be specified
    in terms of reads per million by using a value such as 3.2rpm (which
    specifies the threshold to be 3.2 reads per million). Using a rpm
    threshold allows the sensitivity of cluster discovery to be normalized
    between libraries of different sizes. Alternatively, reads per million
    mapped (rpmm) can be specified: A --mincov of 1.2rpmm indicates 1.2
    reads per million mapped is the threshold. rpm is a fraction of total
    library size, while rpmm is a fraction of only the aligned & placed
    fraction of the library.

segmentSeq
is an R-package useing a Bayesian approach for small RNA locus detection. It's however really slow and not developed any further by the main author.


Use paste function with apply families I am trying to assign different window sizes to my SNPs dataset to identify regions under selection.

this is the head of my data

head(snp_ids)
   snp_id   chr  pos
Chr01__912    1  912
Chr01__944    1  944
Chr01__1107   1 1107
Chr01__1118   1 1118
Chr01__1146   1 1146
Chr01__1160   1 1160

class(snp_ids)
data.frame

I have chosen 4 different window sizes, win_size <- c(15000, 30000, 50000, 100000).
I have assigned each of these different window sizes to my snp_ids dataset to identify how many SNPs are distributed within each window by looping through each window size

for (i in 1:length(win_size)){
  windows <- sapply(snp_ids$pos, function(x) (ceiling(x/win_size[i])))
}

Not clear what we are trying to do, but here is the solution using lapply that should replicate your existing forloop.

# example input data
snp_ids <- read.table(text = "
snp_id  chr  pos
Chr01__912   1  912
Chr01__944   1  944
Chr01__1107  1 1107
Chr01__1118  1 1118
Chr01__1146  1 1146
Chr01__1160  1 1160", header = TRUE)

# window sizes to loop through
win_size <- c(15000, 30000, 50000, 100000)

res <- cbind(snp_ids,
             data.frame(
               lapply(setNames(win_size, paste("window", win_size, sep = "_")), function(w)
                 as.numeric(paste(snp_ids$chr, ceiling(snp_ids$pos/w), sep = "."))
               )))

# result
res
#        snp_id chr  pos window_15000 window_30000 window_50000 window_1e.05
# 1  Chr01__912   1  912          1.1          1.1          1.1          1.1
# 2  Chr01__944   1  944          1.1          1.1          1.1          1.1
# 3 Chr01__1107   1 1107          1.1          1.1          1.1          1.1
# 4 Chr01__1118   1 1118          1.1          1.1          1.1          1.1
# 5 Chr01__1146   1 1146          1.1          1.1          1.1          1.1
# 6 Chr01__1160   1 1160          1.1          1.1          1.1          1.1

Edit: If we really want to use forloop, then try below:

for (i in win_size){
  windows <- ceiling(snp_ids$pos/i)
  snp_ids[, paste0("window_", i)] <- as.numeric(paste(snp_ids$chr, windows, sep = "."))
}


checks for spike-in sequence controls I would like to know what do people verify when designing/using spike-in controls, to be used in sequencing experiments (mainly Illumina). So far I came up with this list:

Does it align only to a given genome synthetic reference?
Does it contain G-quadruplexes? E.g. QGRS Mapper
Does it contain known Illumina error motifs? (Do this matter anymore for recent platforms/chemistries?) See this article
Does it contain other motifs known to be polymerase inhibitors? Ref

Anything else?

Disclaimer: I'm a developer for http://www.sequin.xyz.

Sequin is a new set of spiked-in controls for next-generation sequencing, and that includes Illumina. We design controls for RNA-Seq, genome sequencing, metagenomics etc.

Please study our papers if you want more details.

Reference Standards For Next-Generation Sequencing by Hardwick should be a good start.

Here's my list:

Measured counts against input concentration. Pearson's correlation, spearman correlation and regression slope.
PCA analysis
Whether the synthetic aligned to the synthetic reference
Spiked-in amount (e.g. dilution)
Diagnostic statistics
ROC Curve
Detection limit

(1) is important because the controls should give you correlation = 1, and slope = 1 in a perfect flawless experiment. This is a sample output (you can read more about it in our papers):

(2) This is our sample PCA (our paper has the details)

(3) is also important but you can only do that in simulations.

(5)

(6)

(7)


Are duplicate variants against the VCF standard? I have been trying to understand an error that the EBI's vcf_validator gives when run on my vcf file. Consider this minimal file:

##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  .   AACAC   AATAC   0.00    .   .   GT  0/1
chr12   529516  .   C   T   .   PASS    .   GT  0/1

If I then run vcf_validator on it, it complains about:

$ cat $(LC_ALL="C" ~/vcf_validator_linux  -i test.vcf |& awk '/written to/{print $NF}' )
According to the VCF specification, the input file is not valid
Error: Duplicated variant chr12:529516:C>T found. This occurs 2 time(s), first time in line 6.

This is because the two variants are indeed the same. Changing
AACAC to AATAC at position 529514 just means change C to T at position 529516.

OK, but how is that against the VCF standard? What makes that an invalid VCF file? I am trying to automatically validate input VCF files and stop the analysis if they're invalid, but this seems like a perfectly valid file. Why is  vcf_validator not letting it through? The vcf-validator perl script from vcftools doesn't complain about this. So is it just a bug in  EBI's vcf_validator or is it really somehow illegal to have duplicate variants in a VCF file according to the specs?

No. The only point the specification makes is that:

  It is permitted to have multiple records with the same POS

I suspect this is intended to let the file list different variants occurring at the same position, but this incidentally allows duplicate variants.

This said, you probably don't want duplicate variants, because it will mess with any statistics you apply over the variants, and will unnecessarily clutter the file and slow down processing, if this happens a lot. You might also get variant annotations on one variant but not another, which can become a problem depending on how you query your variants.


Use a specific python.exe instance with PyMOL/Windows? I want to use RDKit with a python script in PyMOL. RDKit requires a conda environment, so I think I need to be able to specify that the instance of python that I want to use would be in that environment, instead of the default for PyMOL, so that I can run it through PyMOL. The script I made is working for me through PyCharm when I test it, by using the RDkit environment "rdkit-env" as my python interpreter, with all of the additional modules I need installed there as well.

I would like to run it through PyMOL so I can add some PyMOL specific functionality as well. I've tried using shebang to specify the RDkit environment that I made, but that does not seem to be working for me with Windows.

Could someone please help me figure out what other approach might work for this? I'm new to this area, so I would appreciate any insights. Thank you!

So it is really hard to make a given python instance run a library from a different instance. It is not just one dist-packages library folder, but multiple that you would need to link, and all that assuming they use the same version of python.
So the best bet is to install a new instance of pymol the python module.

The module can be installed in a python that uses conda with

conda install -c schrodinger pymol

This is the incentive version of PyMOL 2 and does not require a licence.
You installed rdkit, so I am assuming you used conda as conda install -c rdkit rdkit. But just in case you went through the pain of compiling rdkit for a given python instance without conda and want the same for PyMOL. There are guides of instructions required to compile open source. Such as https://pymolwiki.org/index.php/Windows_Install. On a barebone OS like a Pi's Raspsian, I had to do a lot more than claimed say sudo apt-get install python3-dev libglm-dev freeglut3-dev libglew-dev libpng12-dev libfreetype6-dev build-essential libxml++3.6-dev and link a few things manually. So just like for rdkit, conda make life a breeze.

PS. RDKit and PyMOL? Add PyRosetta and you just my favourite modules!


What is the actual cause of excessive zeroes in single cell RNA-seq data? Is it PCR? First, sorry if I am missing something basic - I am a programmer recently turned bioinformatician so I still don't know a lot of stuff. This is a cross post with a Biostars question hope that's not bad form.

While it is obvious that scRNA-seq data contain lots of zeroes, I couldn't find any detailed explanation of why they occur except for short notices along the lines of "substantial technical and biological noise". For the following text, let's assume we are looking at a single gene that is not differentially expressed across cells.

If zeroes were caused solely by low capture efficiency and sequencing depth, all observed zeroes should be explained by low mean expression across cells. This however does not seem to be the case as the distribution of gene counts across cells often has more zeroes than would be expected from a negative binomial model. For Example the ZIFA paper explicitly uses a zero-inflated negative binomial distribution to model scRNA-seq data. Modelling scRNA-seq as zero-inflated negative binomial seems widespread throughout the literature.

However assuming negative binomial distribution for the original counts (as measured in bulk RNA-seq) and assuming that every RNA fragment of the same gene from every cell has approximately the same (low) chance of being captured and sequenced, the distribution across single cells should still be negative binomial (see this question  for related math).

So the only remaining possible cause is that inflated zero counts are caused by PCR. Only non-zero counts (after capture) are amplified and then sequenced, shifting the mean of the observed gene counts away from zero while the pre-PCR zero counts stay zero. Indeed some quick simulations show that such a procedure could occasionally generate zero-inflated negative binomial distributions. This would suggest that excessive zeroes should not be present when UMIs are used - I checked one scRNA-seq dataset with UMIs and it seems to be fit well by plain negative binomial.

Is my reasoning correct? Thanks for any pointers.

The question
How can we distinguish between true zero and dropout-zero counts in single-cell RNA-seq? is related, but provides no clues to my present inquiry.

It may be necessary to distinguish between methods that use unique molecular identifiers (UMIs), such as 10X's Chromium, Drop-seq, etc, and non-UMI methods, such as SMRT-seq.  At least for UMI-based methods, the alternative perspective, that there is no significant zero-inflation in scRNA-seq, is also advocated in the single-cell research community.  The argument is straight-forward: the empirical mean expression vs. dropout rate curve matches the theoretically predicted one, given the current levels of capture efficiency.

Examples

Svensson Blog

A couple of blog posts from Valentine Svensson argue this point rather pedagogically, and include citations from across the literature:

Droplet scRNA-seq is not zero-inflated

Count-depth variation makes Poisson scRNA-seq data negative binomial

bayNorm

There is a more extensive preprint by Tang, Shahrezaei, et al. (BioRxiv, 2018) that claims to show a binomial model is sufficient to account for the observed dropout noise.  Here is a snippet of a relevant conclusion:

  Importantly, as bayNorm recovered dropout rates successfully in both UMI-based and non-UMI protocols without the need for specific assumptions, we conclude that invoking zero-inflation models is not required to describe scRNA-seq data.  Consistent with this, the differences in mean expression levels of lowly expressed genes observed between bulk and scRNA-seq data, which were suggested to be indicative of zero-inflation, were recovered by our simulated data using the binomial model only.

Multinomial Modeling

There is also a very clearly written preprint by Townes, Irizarry, et al. (BioRxiv, 2019) where the authors consider scRNA-seq as a proper compositional sampling (i.e., multinomial process) and they come to a similar conclusion, though specifically for UMI-based methods. From the paper:

  The multinomial model makes two predictions which we verified using negative control data. First, the fraction of zeros in a sample (cell or droplet) is inversely related to the total number of UMIs in that sample. Second, the probability of an endogenous gene or ERCC spike-in having zero counts is a decreasing function of its mean expression (equations provided in Methods). Both of these predictions were validated by the negative control data (Figure 1). In particular, the empirical probability of a gene being zero across droplets was well calibrated to the theoretical prediction based on the multinomial model.  This also demonstrates that UMI counts are not zero inflated.

Furthermore, by comparing raw read counts (prior to UMI-based deduplication) and UMI counts, they conclude that PCR is indeed the cause of zero-inflation:

  The results suggest that while read counts appear zero-inflated and multimodal, UMI counts follow a discrete distribution with no zero inflation (Figure S1). The apparent zero inflation in read counts is a result of PCR duplicates.

I highly recommend giving this a read, especially because it nicely situates other common generative models (e.g., binomial, Poisson) as valid simplifying assumptions of the multinomial model.

It should be noted that this same group previously published a work (Hicks, Irizarry, et al. 2018), mostly focused on non-UMI-based datasets (SMRT-seq), where they showed evidence that, relative to bulk RNA-seq, there was significant zero-inflation.


GSEA enrichr with 10x genomics differential_expression ranks I am attempting to use GSEA enrichr with 10x genomics differential_expression rankings. Reading what people seem to be doing with GSEA, there seems to be a pre-/post-singlecell gene expression versioning of the tools associated, so I want to make sure what I am doing makes sense.
Explanation:
For each gene in columns 1 and 2, there are 3 columns in each of the clusters that the 10X Genomics cellranger count software produces, e.g.:
Cluster 3 Mean Counts   Cluster 3 Log2 fold change   Cluster 3 Adjusted p value
0.2392551160479023      -6.591931364165717           0.029177998608334674
0.1435530696287414      -5.608385347948343           0.0003833420249831818
1.107409394278862       -5.441829814465051           3.8059567638748296e-09
0.1948220230675776      -4.9155704967404             0.008739337761029243
0.20849374398460058     -4.65702722396748            1.0980978990257375e-11
0.36571853453036496     -4.6498340596466585          1.1896434946397147e-17
0.1811503021505546      -4.445327449243366           1.3416233207430882e-12
0.4409129995739914      -4.433552335619653           1.2400396349891802e-16
1.7226368355448967      -4.240367394494225           2.2072280062348652e-15
0.10937376733618391     -3.902215199856403           2.0315030904651526e-12
3.2948847410025404      -3.7580328996321235          0.028858857004902847
0.3588826740718535      -3.738237134795285           8.439629283346646e-13
[...]

In trying to give a tissue / cell type identity to each of these clusters, I've used the following method to feed a subset of genes to EnrichR:
1- Sort the Log Fold change column numerically so we get low negative values.
2- Filter out for the Adjusted p-value to have a value smaller than 0.05.
3- Filter out for the Means Counts value to be bigger than 0.01.
4- Take the top n genes from the list, and feed it to GSEA EnrichR to compare against Mouse_Gene_Atlas or ARCHS4_Tissues.
Given that the fold change of the clusters is relative to each other, does it make more sense to rank the top n genes numerically or in inverse numerical value? E.g. selecting the first 100 genes in each cluster that have the most negative Log2 fold change values or the most positive Log2 fold change values or both (the top 100 genes with highest absolute Log2 fold change)?
Looking at previous answers, I get that one shouldn't rank by p-value: e.g. see @gringer reply
https://bioinformatics.stackexchange.com/a/13465/180

You should be using the most positive fold changes, not the most negative. The positive ones are the ones that are most highly expressed in your cluster of interest. You probably want to use less than 100 genes. Depending on the experiment, you may not have that many good markers for a particular cluster.
The other question you mentioned refers to DESeq2 specifically, so it is not applicable here.


Analysis of differential transcript usage (DTU) Recent breakthroughs in bioinformatics tools for quantification (e.g.
Cufflinks/Kallisto/Salmon etc.) and tools which can identify differential transcript usage (DTU) (e.g. DRIMSeq, Cufflinks etc.) mean that from RNA-seq data we can now relatively easy obtain a genome wide analysis of transcripts that are differentially used betwen conditions.

What can you use these results for? What systematic analysis does this transcript level data enable?

I will take the liberty of giving one possible answers to my own question – but I’m very interested in other answers.

One analysis type that such data enables is the analysis of transcript switches with predicted potential consequences.

I myself have recently developed such a tool called IsoformSwitchAnalyzeR. IsoformSwitchAnalyzeR enables statistical identification (via DRIMSeq) of isoform switches with predicted functional consequences. The consequences analyzed can be chosen from a long list which includes gain/loss of protein domains, signal peptides changes in NMD sensitivity etc. The R package also enables easy visualization of isoform switches along with their consequences and it directly supports the output of Cufflinks/Cuffdiff, RSEM, Salmon and Kallisto.

Apart from enabling identification of interesting examples switch identification also enables systematic analysis of what genes are affected. For inspiration I will recommend one of my own articles (describing results obtained with IsoformSwitchAnalyzeR) as well Hector et al’s recent paper (which is not using IsoformSwitchAnalyzeR). Of particular interest and finesse is Hector’s analysis of how isoform switches can disrupt protein-protein interactions.

Looking forward to hear more ideas.


Why doesn't the mouse GRCm38/mm10 refGene genome annotation file contain non-coding transcripts? I am new to Bioinformatics and I am exploring the refGene.txt files from the UCSC genome annotation database for several species.

My question concerns the Dec. 2011 (GRCm38/mm10) assembly of the mouse genome. I have seen that the Human one (hg39) contains both coding and non-coding transcripts. But the Mouse (mm10) RefGene.txt only contains coding transcripts.

Why is that?

Source of the file: http://hgdownload.soe.ucsc.edu/goldenPath/mm10/database/refGene.txt.gz

It wasn't the correct file.
I downloaded the file again from UCSC, compared it with the file that I was using, and they were different. This new file contained both information from coding and non-coding transcripts.
To all the kind people that took interest in my problem, I am sorry.
Edited 30/12/2020:
To add a bit of clarity; the link in the question is the correct link from which to download the correct file. The file with which I was working from the beginning was the wrong one, but I don't really know where I got it from the first time.


Hisat2 : which option should mention for strand specific library read I was trying hisat2 I get confused by the strand options...
The question has been already asked on github here but didn't get any satisfactory answer.
I'm surprise how much this strand information could be confusing in different tools.

Looking at the hisat2 manual I get confuse by the two different options:

--rna-strandness

  Specify strand-specific information: the default is unstranded.
  For single-end reads, use F or R. 'F' means a read corresponds to a transcript. 'R' means a read corresponds to the reverse complemented counterpart of a transcript. For paired-end reads, use either FR or RF.
  With this option being used, every read alignment will have an XS attribute tag: '+' means a read belongs to a transcript on '+' strand of genome. '-' means a read belongs to a transcript on '-' strand of genome.

and --fr/--rf/--ff

  The upstream/downstream mate orientations for a valid paired-end alignment against the forward reference strand. E.g., if --fr is specified and there is a candidate paired-end alignment where mate 1 appears upstream of the reverse complement of mate 2 and the fragment length constraints (-I and -X) are met, that alignment is valid. Also, if mate 2 appears upstream of the reverse complement of mate 1 and all other constraints are met, that too is valid. --rf likewise requires that an upstream mate1 be reverse-complemented and a downstream mate2 be forward-oriented. --ff requires both an upstream mate 1 and a downstream mate 2 to be forward-oriented. Default: --fr (appropriate for Illumina's Paired-end Sequencing Assay).

Should both need to be specified ? Are they redundant ?
I didn't find so many person struggling with that, so I guess there is something obvious that I don't get. I hope someone may enlighten me about it !

--fr/--rf/--ff should rarely be set, since they refer to the relative orientation of reads and basically everything these days is --fr. Mate-pairs were --rf, but you don't see those much any more.

--rna-strandedness is a very different option, since it sets how reads are expected to align against genes. Most stranded protocols in use these days follow the dUTP-method, where read #2 in a pair has the same orientation as the transcript from which it arose. So either R or RF would typically be appropriate, unless the library is unstranded. In practice, I expect this is more useful if you plan to run stringTie downstream, since then the XS auxiliary tag is set appropriately.


Fast way to count number of reads and number of bases in a fastq file? I am looking for a tool, preferably written in C or C++, that can quickly and efficiently count the number of reads and the number of bases in a compressed fastq file. I am currently doing this using zgrep and awk:

zgrep . foo.fasq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'

The zgrep . will print non-blank lines from the input file and the awk 'NR%4==2 will process every 4th line starting with the second (the sequence).
This works fine, but can take a very long time when dealing with large files such as WGS data. Is there a tool I can use (on Linux) that will give me these values? Or, if not, I'm also open to suggestions for speeding up the above command.

It's difficult to get this to go massively quicker I think - as with this question working with large gzipped FASTQ files is mostly IO-bound. We could instead focus on making sure we are getting the right answer.

People deride them too often, but this is where a well-written parser is worth it's weight in gold. Heng Li gives us this FASTQ Parser in C.

I downloaded the example tarball and modified the example code (excuse my C...):

#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
    gzFile fp;
    kseq_t *seq;
    int l;
    if (argc == 1) {
        fprintf(stderr, "Usage: %s <in.seq>\n", argv[0]);
        return 1;
    }
    fp = gzopen(argv[1], "r");
    seq = kseq_init(fp);
    int seqcount = 0;
    long seqlen = 0;
    while ((l = kseq_read(seq)) >= 0) {
        seqcount = seqcount + 1;
        seqlen = seqlen + (long)strlen(seq->seq.s);
    }
    kseq_destroy(seq);
    gzclose(fp);
    printf("Number of sequences: %d\n", seqcount);
    printf("Number of bases in sequences: %ld\n", seqlen);
    return 0;
}

Then make and kseq_test foo.fastq.gz.

For my example file (~35m reads of ~75bp) this took:

real    0m49.670s
user    0m49.364s
sys     0m0.304s

Compared with your example:

real    0m43.616s
user    1m35.060s
sys     0m5.240s

Konrad's solution (in my hands):

real    0m39.682s
user    1m11.900s
sys     0m5.112s

(By the way, just zcat-ing the data file to /dev/null):

real    0m38.736s
user    0m38.356s
sys     0m0.308s

So, I get pretty close in speed, but am likely to be more standards compliant. Also this solution gives you more flexibility with what you can do with the data.

And my horrible C can almost certainly be optimised.

Same test, with kseq.h from Github, as suggested in the comments:

My machine is under different load this morning, so I've retested. Wall clock times:

OP: 0m44.813s

Konrad: 0m40.061s

zcat > /dev/null: 0m34.508s

kseq.h (Github): 0m32.909s

So most recent version of kseq.h is faster than simply zcat-ing the file (consistently in my tests...).


Difference between de novo transcriptome assembly methods I have been looking around (including read the original papers) to understand what is essentially the difference between StringTie in non-reference based mode (de novo) and Trinity de novo assembly. I understand that in the genome-guided nature of StringTie, we are supported by the genome information (i.e. taking into account the mapped reads). In addition to that, we also have an option to either run using the gene annotation file or not.

Simply put, am I correct if I say that Trinity is purely looking at the overlapping reads, then assemble them, whereas StringTie (non-reference based) is looking at the mapped reads that are either close or overlapping to each other? Is there any other thing I'm missing here (I understand it's not that simple, but I'm trying to make it as simple/intuitive as possible)?

Also if that's the case, would it be safe to say that generally non-reference based StringTie is preferred over Trinity and reference based StringTie (as gene annotation can lead to some "bias", hence preventing discovery of novel TSS, exon length or even make unnecessary duplicates with small differences)?

Your intuition is correct. stringTie is just looking at clumps of alignments and how they might relate to each other (either due to spliced alignments or proximity). Trinity is doing the more computationally difficult task of finding parts of reads that overlap with each other and trying to link them together into longer sequences.

Whether it makes sense to use something like stringTie or not depends mostly on the quality of the reference sequence. If you're working on an organism with a high-quality reference sequence (pretty much anything with a lot of papers using it) then there's no reason to bother with Trinity. If your reference genome additionally has a high-quality annotation (human, mouse, etc.) then you should use it, as anything "novel" that the program is likely to output otherwise is probably just an artefact. If you either don't have an annotation or it's known to be highly problematic then it makes sense to omit it when using stringTie.

Without seeing evidence to the contrary, I don't buy the argument that using a good annotation is actually biasing the output in a problematic way. For well-studied organisms, anything reported as "novel" these days likely doesn't exist, or is otherwise so fleeting that it's biologically irrelevant. The arguments of bias sound like something that was come up with once someone didn't get as many "novel" results as they were hoping.


How can I use pandas agg here to avoid iteration? can anyone help me improve this pandas code?
import pandas as pd

df = pd.DataFrame(
    [
        [
            'chr1', 222
        ],
        [
            'chr1', 233
        ],
        [
            'chr1', 2123
        ],
        [
            'chr2', 244
        ]
    ], columns = ['chrom', 'pos']
)

df2 = pd.DataFrame(
    [
        [
            'chr1', 221, 223
        ],
        [
            'chr1', 230, 240
        ],
    ], columns = ['chrom', 'start', 'end']
)

Gives me 2 dfs with genomic coordinates. The first one is an exact position:
    chrom   pos
0   chr1    222
1   chr1    233
2   chr1    2123
3   chr2    244

and the second is a range:
    chrom   start   end
0   chr1    221 223
1   chr1    230 240

I need to find the count of exact coordinates that are in one of the ranges (in the same chrom)
This works but is slow:
c=0
for chrom, data in df.groupby('chrom'):
    tmp = df2.query(f'chrom == "{chrom}"')
    for p in data.pos:
        for s, e in zip(tmp.start, tmp.end):
            if s < p < e:
                c+=1

Then c = 2
I think I can use agg to do this without iteration (and hopefully faster) but I can't get it working. Can anyone show me how?
PS I am also asking this on stackoverflow.


Finding homologs of a protein sequence I have a refseq ID of a protein from E.coli and I want to find homologs of this protein. I ran Blast against refseq database but I got a lot of sequences most of which were from Ecoli again. I decided to run PSI-Blast to get more divergent species, but I do not exactly know if my result are real homologs or false positives. what is your idea for finding homologous protein sequences from more divergent species? And what can I do for selecting the real (but not false positive) hits?

It depends on what you're looking for. If you're just looking for sequence homology, then you can simply pick the best hits from a blast search. If, however, you are referring to functional homology, if you are looking for the protein which has the same functions as your query, then it's more complicated.

Sequence homology is not enough to infer functional homology. For example, you can have cases of gene duplication and subsequent functional divergence. Such paralogs are still homologs (paralogs are a subset of homologs), but they don't necessarily have the same function. It is also often the case that the homolog (be it orthologous or paralogous) of a protein in species B has a completely different function than its homolog in species A despite a high level of sequence similarity. This is usually very hard to determine in silico.

To find the functionally homologous protein(s), you would ideally need to identify the essential residues that allow your protein to perform its function. This could be done using something like PFam which will identify protein domains. You can then check whether the homologs you find also have this domain.

This is essentially what PSI-blast does. Although it doesn't take domains into account, each successive iteration is used to build a model of your proptein. The model is built under the reasonable assumption that highly conserved residues are important. So it will consider more diverged sequences as homologous if those residues are conserved.

If you know how your protein works and what residues are important, you can use that knowledge to refine the results of your PSI-blast. If you don't, you'll have to use only "good" hits to make the model. One way to do this, for well studied proteins, is to only add proteins that are already annotated as homologs of what you are looking for to build your model, then use that model to search in un-annotated species.

If you don't know, you could try looking for recognizable protein domains in your query protein (use PFam) and then use the HMM (hidden markov model) of the domain to identify important residues. For example, this is the HMM logo for the zf-A20 zinc finger domain:

The huge cysteine (C) residues are shown that size because they are very conserved across proteins carrying this domain and, presumably, are functionally important for the domain. So, if you pass your protein through PFam and identify domains, find the important residues and make sure all your homologs have those conserved. If using PSI-blast, only include sequences where those residues are conserved in the results you keep.

Finally, another useful tool that works in the same way is HMMER. This takes a protein alignment as input, like PSI-blast builds an HMM model from it and then can use this model to query a protein database for more hits. Methods like HMMER and PSI-blast are far better than simple sequence similarity approaches when looking for homologs.


How do I re-name the headers of my Fasta file? I appologize for asking this, but I really am really bad with regex...

Can someone help me transform the headers of my fasta files from this:

>Agaricus_chiangmaiensis|JF514531|SH174817.07FU|reps|k__Fungi;p__Basidiomycota;c__Agaricomycetes;o__Agaricales;f__Agaricaceae;g__Agaricus;s__Agaricus_chiangmaiensis

to this:

>JF514531|SH174817.07FU;tax=d:Fungi,p:Basidiomycota,c:Agaricomycetes,o:Agaricales,f:Agaricaceae,g:Agaricus,s:Agaricus_chiangmaiensis_SH174817.07FU;

seems like I could use the command

sed 's/__/:/g;s/;/,/g;s/.*|/|/g'

to substitute : for __, and , for ;.  I'm not entirely sure how to proceed from there though...

Thank you

If awk is good as well:

$ awk -v FS="|" -v OFS="|" '/>/{
    gsub(";",",", $5);
    gsub("__",":",$5);
    print ">"$2,$3";tax="$5"_"$3";";
    next;
}1' file.fa
>JF514531|SH174817.07FU;tax=k:Fungi,p:Basidiomycota,c:Agaricomycetes,o:Agaricales,f:Agaricaceae,g:Agaricus,s:Agaricus_chiangmaiensis_SH174817.07FU;

This relies on the order of the parts separated by |.

I don't understand why k__Fungi becomes d:Fungi, so I ignored it for now.

Some explanation for what awkis doing.

| is defined as a field (column) separator for the input and output
If the line contains a > the code between { and } is executed.
the first gsub replaces ; with , in the 5th field
the second gsub replaces __ with : in the 5th field
then we print the new header line. The $2, $3 and $5 get replaced by the values of the second, third and 5th field of the input. , get replaced by the defined output field separator
with next we tell awk to go to the next input line and repeat the code
The 1 is in equivalent to print $0 for printing the whole input line. In our code this is only executed if the line is no header line.


Can a data file in VCF format be converted into FASTA? I'm considering purchasing the 'MyGenome' product by Veritas Genetics to analyze my genome for a project. I'd like the data to be in FASTA format, but Veritas only provides VCF data. Is it possible to convert this VCF data into FASTA format?

You can try gatk function FastaAlternateReferenceMaker

java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -L input.intervals \
   -V input.vcf \
   [--snpmask mask.vcf]


panic: division of zero by zero or infinity by infinity error when running GODON positive selection analysis I am using godon to find genes under positive selection with the branch-site test. However, after running for a little while I get the error "panic: division of zero by zero or infinity by infinity".
Below is the full output with the error report. I can supply the phylogenetic tree and fasta file I am using if that would be helpful. Thank you.

iter=33 lnL=-73104.072
iter=66 lnL=-72956.855
iter=99 lnL=-72912.956
iter=131 lnL=-72902.900
iter=162 lnL=-72900.147
iter=194 lnL=-72898.562
iter=226 lnL=-72897.970
iter=258 lnL=-72897.651
iter=289 lnL=-72897.435
iter=321 lnL=-72897.284
iter=353 lnL=-72897.199
iter=385 lnL=-72897.153
iter=416 lnL=-72897.108
iter=448 lnL=-72897.076
iter=17 lnL=-72890.551
panic: division of zero by zero or infinity by infinity
goroutine 1 [running]:
math/big.(*Float).Quo(0xc4273a8f40, 0xc4273a8f68, 0xc4273a8f18, 0xc42e88c118)
/usr/local/go/src/math/big/float.go:1639 +0x15c
bitbucket.org/Davydov/godon/cmodel.(*BranchSite).BEBPosterior(0xc422164120, 0xc4259bb900, 0x91, 0x91)
/go/src/bitbucket.org/Davydov/godon/cmodel/bs.go:416 +0x6e9
bitbucket.org/Davydov/godon/cmodel.(*BranchSite).Final(0xc422164120, 0x101)
/go/src/bitbucket.org/Davydov/godon/cmodel/bs.go:455 +0x112
main.computeFinal(0x1f263a0, 0xc422164120, 0xc4201522a0, 0x101, 0x6a80, 0x0)
/go/src/bitbucket.org/Davydov/godon/godon/optimize.go:150 +0x93
main.performSingleTest(0xc420162320, 0xc428266000, 0x9399, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, ...)
/go/src/bitbucket.org/Davydov/godon/godon/hyptest.go:335 +0x1b52
main.hypTest(0x1ee03a7, 0x1ee03a7, 0x4, 0xc4201aa501, 0x1, 0x1)
/go/src/bitbucket.org/Davydov/godon/godon/hyptest.go:88 +0x71a
main.main()
/go/src/bitbucket.org/Davydov/godon/godon/godon.go:316 +0xefa

Ah wait you mentioned a branch sites test, right that makes sense. You've got two or more identical sequences in the analysis. Ensure there are no identical sequences in your alignments and ... PRESTO! It will work :-)
Below are my original thoughts before I re-read your question. I've kept the notes because it helps explain the theory a bit. There is no need to supply the tree, it should be resolved.

The tree likelihood looks to be decreasing nicely, so it can see both the tree and the codons in the alignment.
At a guess the error relates to two or more indentical nucleotide sequences in your alignment. Thus zero non-synonymous sites (amino acid change) and zero synonymous (silent) sites. This error would only occur if you were calculating selection for each 'node' in the tree (internal branch).  Alternatively, possibly if the number of synonymous sites was zero, regardless of the non-synonymous mutations this error could appear, but this doesn't seem to fit the error output.
I don't know "Godon", but recommend PAML for this analysis. Admittedly it is complicated to deploy, but it is a very robust algorithm and hugely respected.


Post docking evaluation What to do after protein-ligand docking? Is there any necessary or enhancement step recommended after docking to validate the confirmation of highest score is the best? For your information, I am using CCDC GOLD for docking and GOLD score as my scoring function.

You can run a Molecular Dynamics simulation to see if the pose is stable and whether the protein changes/adapts it’s conformation.

You can calculate the ligand efficiency, the docking score divided by the number if heavy atoms, and rank the compounds accordingly.

You can also compare the scores of known ligands to the experimental values and see if the docking the scores are predicting known ligands or correlate with the experimental ones.

At some point you have to make a selection and purchase or synthesize the compounds. How many compounds can you test experimentally? How many can you purchase? What is the purpose of these ligands going to be? Drugs? Did you screen the correct set of ligands e.g. a library of lead compounds that are ready to purchase for cheap, that would be a good start.