How to show only my read on IGV? I have a large BAM file on my IGV. I'd like to visualize a read alignment. I have the read name, and I don't care anything else in the file.

IGV has an option "Select by name" (see screenshot). The problem is that IGV doesn't hide the alignments not matching the name. I'd still have to manually navigate on IGV with my mouse to find my read.

That was too cumbersome. What I'm  doing now is make a new alignment BAM file. For example,

samtools view -h L.V.51.bam | grep '@\|HWI-D00119:228:CB08RANXX:8:2116:4095:45403' | samtools view -bS - | samtools sort - > new.bam

I'd load the new alignment file on my IGV.

If I have 20 reads, I'd have to create 20 new BAM files... Is there anything I can do on IGV that I'm not aware of?

Partially a shameless plug... I've written a genome browser, ASCIIGenome, entirely controlled through the command line. As such, it makes relatively easy to deal with questions like yours. For example, you could do:

ASCIIGenome bam/WW00292a.bam  # Load bam file(s)

Filter for reads matching the regex passed to -i option:

[h] for help: grep -i '72789|57688|70785|56299|27116|12086'

(Optionally, show sam records as strings)

[h] for help: print

Example result:

For command help type grep -h and for more filtering options see also awk -h.

Hope this helps and feel free to post questions and open issues on github!


How can the cell line contribution be estimated from RNASeq data? Using a laser-capture microdissection of cells a group of cells stained with the marker of interest was sequenced. In another cohort of patients (this is all  human liver tissue) the whole tissue was sequenced (RNA-seq in both cases)

Can I estimate the contribution of the cells marked in the whole liver ("weight of these cells" in the liver in words of my PI)?

My gut feeling is that it can't be done this way, it would require both single cell sequencing and whole tissue sequencing to estimate the contribution of each cell line. But perhaps there is a tool that given the cell lines or the main expression of other cells lines it can be compared to using GSVA or some similar tool.

There are couple of computational methods which try to do this (I never used them, so no experience):

CellMix, based on sets of marker gene lists
Subset Prediction from Enrichment Correlation, which is based on correlations with subset-specific genes across a set of samples.
Cell type enrichment, which uses our highly expressed, cell specific gene database
Cell type-specific significance analysis using differential gene expression for each cell type

You might have to get some reference expression levels from public databases or papers for some of the methods.

One thing to keep in mind: you cannot really compute the cells proportion, only RNA proportion. If you have a good reason to assume that RNA quantity
per cell is very similar, this is a good proxy for the cells proportion in a tissue.


DESeqDataSetFromTximport invalid rownames length I am trying to use DESeqDataSetFromTximport function from DESeq2 package to construct dds object:

dds <- DESeqDataSetFromTximport(txi, sampleTable, ~Group)

And somehow it is giving me the following error:

  Error in rownames<-(*tmp*, value = c("ENSMUSG00000000001", "ENSMUSG00000000003",  :
    invalid rownames length

I used the same code for other datasets but for the current one it is somehow not working. My sampleTable looks like this:

       ID      Group
1    de_1 denervated
2    de_2 denervated
3    de_3 denervated
4    de_4 denervated
5   den_5       sham
6    ho_1        IgG
7    ho_2        IgG
8    ho_3        IgG
9    sh_1       sham
10   sh_2       sham
11   sh_3       sham
12   sh_4       sham
13   sh_5       sham
14 VLA_10  anti-VLA4
15 VLA_12  anti-VLA4
16  VLA_7  anti-VLA4

txi was formed in the following way:

k = keys(txdb, keytype = "GENEID")
df = AnnotationDbi::select(txdb, keys = k, keytype = "GENEID", columns = "TXNAME")
tx2gene = df[, 2:1]  # tx ID, then gene ID
txi <- tximport(names, type = "salmon", tx2gene = tx2gene)

Any suggestions would be greatly appreciated

I got exactly the same problem not long ago (see my other question what I tried to do). We used the same function!

My code ran fine on my HPC but the same code/data failed on my laptop. I later found out I had to upgrade my R to the latest version (3.5.1) and upgrade all my bioconductor packages to the latest. My HPC was on a recent bioconductor docker, so there was no problem.

Please make sure everything in your environment latest (R, R packages, Bioconductor).


GMAJ error: Exception in thread "main" java.lang.ClassCastException: java.lang.StringBuffer cannot be cast to java.lang.String I have installed GMAJ using:

wget http://globin.bx.psu.edu/dist/gmaj/gmaj.zip  -O output/software/gmaj.zip
unzip -u output/software/gmaj.zip -d output/software/gmaj/
rm output/software/gmaj.zip

I ran lastz on two chromosomes from two different assemblies of C. elegans:

lastz c_elegans.PRJNA13758.WS263.genomicI.fa I.fa --format=maf > MAF.maf

and got a MAF file.

Then, I tried to visualise the file using:

java -jar output/software/gmaj/gmaj.jar MAF.maf

but I got the following error:

Exception in thread "main" java.lang.ClassCastException: java.lang.StringBuffer cannot be cast to java.lang.String
    at edu.psu.bx.gmaj.MajGui.setDefaults(MajGui.java:272)
    at edu.psu.bx.gmaj.Maj.<init>(Maj.java:53)
    at edu.psu.bx.gmaj.MajMain.main(MajMain.java:86)

Can someone help me resolve this error or suggest an alternative way to visualise my alignment?

I work on Mac OS X.

Cheers.

It's a programming error that should look like this:

StringBuffer sb = (...)
(...)
String s = (String)sb;

while it should be

StringBuffer sb = (...)
(...)
String s = sb.toString();

( see also https://www.biostars.org/p/177992/ ) . You cannot do anything about it but asking the authors of the software.


Quantification of EST/cDNA based retained intron transcripts Note: this question has also been asked on Biostars

We are seeing a retained intron transcript event for some RNAseq samples, and we want to assess at the sequence level which intron retention events are actually occurring so that we can do protein domain predictions. For example, we have transcript level expression support for DROSHA-203, which according to Ensembl is a retained intron transcript that potential gives rise to multiple splicing variants via the cDNA sequence information:

In theory, the orange sites indicate where a mutation might occur which could obliterate a splicing recognition site, thereby causing an retained intron event. So in total there are 6 potential mutational events which might give rise to 4 possible retention events. Given that 4 intronic regions might be retained, the number of unique retention events is a combinatorics question equal to 2^4, e.g., intron1 retained only, vs intron1 + 2 retained, all possible unique combination, etc, correct?

Is there a better way to more granularly quantify these possible transcripts? Or do we have to PCR validate all of them to figure out which one's are actually occurring? Is there a way that I can look at the reads mapping to the corresponding intronic regions to make even a probabilistic assessment of which transcripts are present and to what extent they are expressed? The libraries are prepped for total RNA with a ribozero depletion, so I should still be able to quantify any enrichment to these intronic regions, correct?

You can do full-length whole-transcriptome sequencing using either the direct RNA kit or a cDNA kit (with or without PCR) on a MinION. This will allow you to quantify isoforms directly, rather than estimating by looking at read pairs.

However, it'd be a bit expensive to do it that way if you're looking at a single transcript (and aren't barcoding to multiplex samples on the same run).

If you've already done Illumina whole-transcriptome sequencing (e.g. strand-specific paired-end sequencing), tools like Kallisto, RSEM, or Salmon can be used to estimate isoform expression based on the mapped reads. These tools use the information about transcript placement and/or similarity to alter the observed read counts in an attempt to better represent reality.

Unfortunately, if you want these tools to quantify intronic sequence (i.e. the $2^4$ possibilities you have suggested), then those sequences need to be included in the sequences provided. They won't estimate frequencies for non-existent sequence.

As an alternative to including every possible gene model, something like a SuperTranscript model might work. I'm not familiar with Alicia Oshlack's method, but it at least gives the appearance that it might reduce the complexity of searching for multiple unannotated isoforms within a gene.


How can I compute gene expression for a set of RNA reads? I'm trying to compute a gene expression profile for an organism. I have gene nucleotide sequences of the mentioned organism stored in a fasta file and a set of paired reads stored in two separate files with the same fasta format. Now I want to compute the coverage of every gene by the given reads. That is, on average how many reads cover a given segment of the sequence.

Is there any tool for this specific purpose?
I'm trying to use the velvet sequencer but I have encountered a few problems.

It's unclear if your paired-end reads are actually in fasta format, I'll presume that they're in fastq instead.

The easiest tool to use is salmon, which nicely deals with things like multimapping. If you're trying to judge the quality of an assembly, then I recommend having a look at transrate, which uses some related methodology for assessing contig quality.


Is there a resource where I can donwload E. coli assemblies grouped by their pathotypes? I will preface this by saying I'm not a microbiologist, so I apologise for any haziness and I'll be happy to expand.

I would like to download the assemblies of different E. coli genomes and group them by their pathotypes (e.g. EPEC, ETEC, EIEC).

I have used enterobase previously to download salmonella genomes broadly by their subspecies, and it does have a set of E. coli genomes. However, they are mainly annotated by their serotype, and not with any broad nomination of pathogenesis. I think I would be able to use the API to grab these serotypes and map them to a type of pathogenicity but that would probably be quite messy.

Does such a database exist? Most I can find are broadly annotated or a are very detailed, but don't allow access to the raw annotation.

Thanks in advance.

If ENA archive contains such information, then you can access it programatically through the EBI api.

Also you can access the NCBI database programatically through their eutils.

In both cases you can access annotations, though not sure if they contain what you are looking for in all the cases.


How to import a phylogenetic tree (in nexus format) with node support values in R I have a phylogenetic tree in nexus format that was generated with RAxML. In the tree file, there are node support values. I want to plot this tree in R with the support values showing. However, when I import the tree file in R with the read.nexus function in ape, the node values do not appear in $node.labels. Can someone recommend another option in R?

I recommend you this tutorial about ggtree. You will find the solution.
[https://guangchuangyu.github.io/ggtree-book/short-introduction-to-r.html][ggtree]


Annotating gene names or gene IDs to a dataframe containing SNPs? I have a large data-frame (excel file) of SNPs with genotyping data. I need to filter the data to get SNP information of a specific gene alone. The list is too long to do it manually. I was wondering of annotating the SNPs with the gene names so that I can filter using the gene name. Does anyone know of a good R package to perform the task?

I assume your SNP file contains genomic positions - please clarify if not. For example, you could have an rsID column, chromsome, and position column that contain "rs3088379", "chr1", "54225281".
Approach #1: use R, as you suggest.
ensembldb is an R package that may be what you're looking for: https://bioconductor.org/packages/release/bioc/html/ensembldb.html.
You could map the genomic positions of your SNPs to gene names by adapting their workflow: https://bioconductor.org/packages/release/bioc/vignettes/ensembldb/inst/doc/coordinate-mapping.html.
After you've mapped the genomic position, you could filter for your gene name.
Approach #2: use bedtools.
You can convert your Excel file of genomic positions to a tab-delimited BED format. https://bedtools.readthedocs.io/en/latest/content/overview.html has an overview. So, your first SNP would turn into
chr1 54225280 54225281 rs3088379

--> make sure you turn 1-based to 0-based coordinates!
You could then download the coordinates of your gene in BED format from your favorite source, like UCSC Table Browser.
Then, you run bedtools intersect with your SNP BED and your gene BED. What remains is just that gene's SNPs.


Aggregate sequencing/mapping/etc. metrics from cellranger across Illumina samples I have a number of single-cell projects processed with cellranger from 10x Genomics. The pipeline produces a number of handy metrics that are summarized for each Illumina sample in web page. These include things like the estimated number of cells, the fraction of reads in cells, and mapping rates. What I would like is to aggregate this information from multiple Illumina samples into a single file. MultiQC is the obvious choice for such a task, but it appears that no one has written a module to process cellranger data for it. Are there any other prepackaged ways to perform this aggregation?

If there's no convenient tool for this I'll just write a MultiQC module, but I'd obviously like to avoid the extra effort unless it's really needed.

I wrote a simple script to do it. Not sure whether that's what you want.
#!/bin/bash
#CellrangerSummary.sh Sample1 Sample2
csvtool transpose $1/outs/*summary.csv | csvcut -c 1 > summary.csv
for sample in "$@"
    do
        csvtool transpose $sample/outs/*summary.csv > tempfilesforcsvtool.$sample
        csvjoin -c 1 summary.csv tempfilesforcsvtool.$sample > temp2
    mv temp2 summary.csv
rm tempfilesforcsvtool.$sample
    done

csvformat -T summary.csv > summary.tsv
rm summary.csv
```


Extract features from fasta sequences and train the classifier I am new to the bioinformatics field. I have positive and negative protein sequences for acetylation PTM. Now, I want to train a classifier, say SVM. What will be the next step? How can I convert these sequences into usable features? Any information or links would help.
>P31327_55|1|testing
KAQTAHIVLEDGTKMKGYSFGHPSSVA
>P31327_57|1|testing
QTAHIVLEDGTKMKGYSFGHPSSVAGE
>P31327_119|1|testing
APDTTALDELGLSKYLESNGIKVSGLL
>P31327_157|1|testing
LATKSLGQWLQEEKVPAIYGVDTRMLT

Just a recommendation: Be careful when making your data sets.
P31327_57, P31327_119, P31327_157 are all annotated as "By similarity" to the mouse entry in Swiss-Prot https://www.uniprot.org/uniprot/P31327#ptm_processing, which means there is no direct experimental evidence in the human entry. If you do have publications that report experimental evidence for these PTMs, please don't hesitate to share them with the UniProt team.
I also recommend reading this UniProt help page about using UniProt data in negative datasets if you haven't already done so: https://www.uniprot.org/help/negative_datasets


BioMart returns one more item I am supplying 30 ensembl ids into BioMart, but somehow it returns 31 result. What is going wrong here?
Here you can see the code:

library('biomaRt')
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
genes <- rownames(mat)
G_list <- getBM(filters= "ensembl_gene_id", attributes=
c("ensembl_gene_id","entrezgene",
"mgi_symbol","description"),values=genes,mart= mart)

genes list is the following one:

[1] "ENSMUSG00000097971" "ENSMUSG00000097582" "ENSMUSG00000061808" "ENSMUSG00000097206"
[5] "ENSMUSG00000022425" "ENSMUSG00000094708" "ENSMUSG00000095937" "ENSMUSG00000032484"
[9] "ENSMUSG00000060550" "ENSMUSG00000073409" "ENSMUSG00000026822" "ENSMUSG00000056071"
[13] "ENSMUSG00000035299" "ENSMUSG00000045932" "ENSMUSG00000006154" "ENSMUSG00000032496"
[17] "ENSMUSG00000040809" "ENSMUSG00000032690" "ENSMUSG00000055322" "ENSMUSG00000066443"
[21] "ENSMUSG00000005089" "ENSMUSG00000073489" "ENSMUSG00000005800" "ENSMUSG00000034855"
[25] "ENSMUSG00000029417" "ENSMUSG00000029561" "ENSMUSG00000021843" "ENSMUSG00000020108"
[29] "ENSMUSG00000037849" "ENSMUSG00000066861"

As you can see I have here 30 genes, but BioMart returns:

> G_list$mgi_symbol
[1] "Slc1a2"  "Mmp8"    "Eps8l1"  "Ddit4"   "Ktn1"    "Enpp2"   "Lcn2"    "Cxcl9"
[9] "Oasl2"   "Ngp"     "Ltf"     "Oas2"    "Cxcl10"  "Mid1"    "Ifi206"  "Chil3"
[17] "Ifit2"   "Tns1"    "S100a9"  "H2-Q7"   "H2-Q7"   "Ttr"     "Gm10163" "Oas1g"
[25] "H2-Q6"   "Ifi204"  "Gm10359" "Gm12671" "Gm26699" "Gm26527" "Gm26917"

If you print G_list, you'll see this half way down:

1 ENSMUSG00000060550      15018      H2-Q7      Mus musculus histocompatibility 2, Q region locus 9 (H2-Q9), mRNA. [Source:RefSeq mRNA;Acc:NM_001201460]
2 ENSMUSG00000060550     110558      H2-Q7      Mus musculus histocompatibility 2, Q region locus 9 (H2-Q9), mRNA. [Source:RefSeq mRNA;Acc:NM_001201460]

So there are two entrez ids that map to the same Ensembl id. One of them maps to an alt locus.

If you exclude "entrezgene" from the query, you'll get 30 results.


Does BLASTClust guarantee that proteins in different clusters are dissimilar? I need to find dissimilar proteins. Looking through the PDB I found the weekly BLASTClust results of proteins that are 30% similar. However, I do not know if protein A in cluster 1 is guaranteed to be less than 30% similar to protein B in cluster 2.
The best and only documentation I can find for BLASTClust is here (link).
It states:

The program begins with pairwise matches and places a sequence in a cluster if the sequence matches at least one sequence already in the cluster.

That's good, but consider a set of three proteins (A, B, and C). A is 40% similar to B and 40% similar to C. B is 20% similar to C. Say that the algorithm places B in cluster 1 and then C in cluster 2 as they are dissimilar. What happens when it tries to place A? Because A is similar to both B and C will they be joined into one cluster (this would be the best result for me) or is A just placed into one cluster even though it is similar to a protein in another cluster?

Based on the documentation (that I linked to in comment), it seems that the clustering algorithm used by blastclust is single-linkage hierarchical clustering.
SO what you can expect is that A will be placed with the protein to which it is most similar. And then the last protein will be added as an outgroup to that cluster due to its similarity to A (assuming that there are no other proteins).
So if your goal is "I want clusters that are guaranteed to contain dissimilar proteins", then you might not be able to naively take blastclust output and interpret it as such. You could either postprocess it to find the set of clusters in which no protein is more than 30% similar to a protein in any other cluster (which may just be one big cluster!), or find another clustering algorithm that is more suited to your task.
Update
It does look like blastclust uses a minimum similarity (not identity) threshold (from blastclust docs). So they use an implementation of single-linkage clustering that stops before merging everything into a single cluster:
If the coverage is above a certain threshold
 AND
the score density is above a certain threshold,

these two sequences are considered to be neighbored.

Thus determined neighbor relationships is considered symmetric and provides
the base for clustering by a single-linkage method (which puts a sequence
to a cluster if the sequence is a neighbor to at least one sequence in the
cluster)."

I don't think that this is a general property of the generic clustering algorithm though (see the "Naive algorithm" section).


How can I download genePred file from UCSC? I need a gtf annotation file. UCSC doesn't give us a proper gtf file with distinct gene_id and transcript_id. It asked us to get a genePred file to convert to gtf.

But where can we get genePred files?

It turned out that using table browser on UCSC we could download an output with all fields. That output actually contains a GenePred file but has one extra column on the left end.

The genePred file format definition is here.


What is the standard way to work with a diploid reference genome? Complementary strands? At the moment, the standard reference genomes (e.g. hg19, hg38) are haploid genomes. We know that the human genome is diploid. Naturally, the latter would be the respectively correct representation of the human genome.

More and more biologists are using emerging technologies into order to capture the diploid nature of genetic information, e.g. phasing SNPs between mother and father chromosomes.

What is the standard way that bioinformaticians generated a standard diploid reference genome?

Actually, reference genomes are not truly haploid (by my understanding). Given that reference genomes are 5'-3' coordinated, in order to create a complementary strand, one would need to take the 3'-5' complement. In order to have a diploid genome, you need two reference genomes and two 3'-5' complements.

More importantly, how have large-scale genomic studies dealt with the fact that the haploid reference genome is a consensus-based "half" of a human genome?

At the moment, the standard reference genomes (e.g. hg19, hg38) are haploid genomes. We know that the human genome is diploid. Naturally, the latter would be the respectively correct representation of the human genome.

The premise of the OP's question is false.  The natural reference representation of the human genome is not diploid.

Think of a reference genome as a map, and not as a specific example of a human being's DNA.

Not only is the human genome reference haploid, but it is also a composite genome. This means that the human genome reference sequence is composed of sequences from multiple individuals.  In other words, the human reference does not correspond to any one human sequence.

Any particular read from a DNA sequencer will be a read from a human genome that diverges from the reference genome. So an algorithm that tries to match the read to the reference genome will always need to handle potential discrepancies. Adding a second map against which to match a read would not change that fact.  Therefore, there is little value in providing a second haploid reference genome.

Side note:
There are parts of the human genome that "are too complex to be represented by a single path", and the Genome Reference Consortium provides "alternate loci" for such regions of the genome.

Answer to the OP's second question

  More importantly, how have large-scale genomic studies dealt with the fact that the haploid reference genome is a consensus-based "half" of a human genome?

I interpret this question as "How do large genomic studies represent genome diversity that cannot be represented by a haploid reference genome?"

One standard approach that has been very popular in genome-wide association studies is the use of a haplotype reference panel such as for example from the 1000 Genomes Project.  Modern phasing programs exploit the shared ancestry of the samples in a study and reference haplotypes to phase a study's samples. One popular phasing program is Impute2


16S species-level taxonomic assignment--what is the current state of the art? We currently use DADA2 for picking ASVs and the assignTaxonomy funciton for assingment to genera. Google does bring up various recent articles on species-level assignment with 16S.

Does anyone have an opinion on the current best-practice?

More simply, are there any current candidates for reasonably accurate species-level assignment?

I have been using dada2 and phyloseq (phangorn for trees) combination to analyze a V4 region 16s rRNA data set. I find their approach and explanations very understandable and usable since I prefer R.
I was under the impression that even full-length 16s rRNA gene sequences were inadequate to classify up to species level since closely related species within a genus can have identical 16S rRNA gene sequences.
But, turns out, intra-genomic variants of this gene within a species or even strain may allow their unambiguous identification according to this paper:
Johnson, J.S., Spakowicz, D.J., Hong, B. et al. Evaluation of 16S rRNA gene sequencing for species and strain-level microbiome analysis. Nat Commun 10, 5029 (2019). https://doi.org/10.1038/s41467-019-13036-1) available at: https://rdcu.be/b5cgB
Another relevant article regarding using black box ML algorithms (random forests) versus a straightforward logistic regression is discussed in this one:
Topçuoğlu BD, Lesniak NA, Ruffin MT 4th, Wiens J, Schloss PD. A Framework for Effective Application of Machine Learning to Microbiome-Based Classification Problems. mBio. 2020;11(3):e00434-20. Published 2020 Jun 9. doi:10.1128/mBio.00434-20
Since the OP has probably already moved on long since the original query, I am writing this to document relevant current references for myself and anyone researching this topic on StackExchange Bioinformatics.


Separating peaks of chip-seq with specific length My data contain several chip-seq results. I have the peaks called by MACS2.I wanted to only look at those peaks that their size is e.g 500bp to 1000bp. How can I separate those peaks efficiently?
I know one way is to take the peaks.xls file sort them based on the length column and take out the peaks less than 1kb. But I was looking for a more efficient way to do that. Any suggestion?
Many Thanks

You can do it by awk.
awk -v OFS="\t" '{ if ( ($3-$2) >= 500 && ($3-$2) <= 1000) { print } }' in.bed > out.bed

Assuming there is no header.


How to interpret contig-alignment.psa produced by velvet I'm using velvet to align given reads of RNA to given CDSs (i.e. coding areas and genes) of an organism, so I can generate gene-expression profiles. But after using velvetg out-dir/ -alignment yes, velvet produces contig-alignment.psa which is a file in a strange format.
The file for every contig contains zero or more records of 7 integers, and I cannot understand their meaning.

The following is an excerpt of the file, and you can see the entire file here:

>contig_1
1    193    1939    1939    6   194 2
195 244 1   1   6   22  71
>contig_2
1   84  1   1   6   59  142
86  170 1935    1935    6   22  106
172 285 1935    1935    6   108 221
>contig_3
1   98  2   2   6   1   98
100 321 2   2   6   100 321
334 415 1204    1204    6   1   82
>contig_15
1   23  3   3   6   84  106

It looks like you are trying to do a transcript isoform expression analysis. Velvet is a genome assembler, and really not designed for this task. The tools you are looking for are kallisto and salmon.


ScaleData() from Seurat causes crash on RStudio Cloud Using the sample data in the 2,700 PBMC clustering tutorial, the session crashes at the ScaleData() step. The progress bar remains at 0% and hangs, followed by an R crash. I've tried reducing the size for number of genes to scale at in a single computation with the argument block.size with no change. The Seurat S4 object is only ~70mb so I can't imagine I'm exceeding the RStudio Cloud 1gb RAM limitations.

If you'd like to check out the cloud instance I can add you to the workspace temporarily. There is a limit on members allowed in a workspace.

I've been trying to run this tutorial on another remote cluster in RStudio inside of a container (Singularity here vs. Docker in RStudio Cloud) and am hitting the same snag.

I get receive no errors from RStudio Cloud, but in the other evirnment I get this error:

ERROR session hadabend; LOGGED FROM: rstudio::core::Error {anonymous}::rInit(const rstudio::r::session::RInitInfo&) /home/ubuntu/rstudio/src/cpp/session/SessionMain.cpp:563.

On this issue thread someone mentions it they receive this error message when hitting their RAM limits.

The source for ScaleData() can be found in the Github repo.

It seems RStudio server doesn't know it's limitations. Currently, Projects are limited to 1GB of RAM. Using the R package ulimit which mirrors the function of ulimit and allows implementation of memory limits in linux I was able to set my own limitations within the R session. After setting the limit to 1GB I see the ulimit limitation being tripped when invoking ScaleData(). Looks like a case of exceeding available system resources.


Ensembl id to GeneSymbol with biomart I have 3224 Ensembl id's as rownames in a dataframe "G". To convert Ensembl ids into Genesymbols I used biomart like following.

library('biomaRt')
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
genes <- rownames(G)
G <-G[,-6]
G_list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id","hgnc_symbol"),values=genes,mart= mart)

Now in G_list I can see only 3200 ensembl ids showing Genesymbols / No Gene_symbols. Why the other 24 ensembl ids are not seen in G_list? If there are no gene_symbol for those 24 ensembl ids it should atleast show "-"

Examples of problematic IDs are: ENSG00000257061, ENSG00000255778, ENSG00000267268. These are not at all shown in G_list (biomaRt). So, I gave them in biodbnet, which seems to handle them.

what is the problem here?

It looks like you were using an old annotation. The problematic IDs you posted existed in the GRCh37 annotations, but don't in the most recent GRCh38 annotation. For that reason they were excluded. The IDs that have - as symbols don't have associated symbols, but are present in the database.

To use an archived version in biomart:

mart = useDataset("hsapiens_gene_ensembl", useEnsembl(biomart="ensembl", version=84))

That's an example for release 84.


tool to visualize a collection of phylogenetic trees I have a collection of phylogenetic trees produced from clustering DNA sequences, mostly a few nodes in each cluster/tree, with several clusters of size 1. And I would like to visualize them all together to get an idea of the relative sizes of the trees in the collection of clusters (cluster set).

E.g.

(ABC:0.0,('CL-118070':0.02010925,'CL-117615':0.00526463):0.00443474);
(DEF:0.0,('CL-117360':0.00967856,'CL-117181':0.00967856):0.00990975);
(GHI:0.0,('CL-123443':0.00123443,'CL-123443':0.00132446):0.00123435);
(JKL:0.0,('CL-134123':0.00123441,'CL-545443':0.00343446):0.00134335);
((((CL-114948:0.00908659,(CL-114989:0.00324875,CL-114959:0.00536289):0.00037038):0.00457230,(CL-114851:0.01291841,CL-114831:0.00590321):0.00077034):0.03411459,XYZ:0.00251364):0.00207594,(CL-114950:0.03074481,CL-114699:0.02256162):0.00733986,CL-116426:0.01643848);

I am thinking there must be a tool where I give it an input file with a tree per line and shows them as all in one go.

Any suggestions:

EDIT: looking at the options below:

https://en.wikipedia.org/wiki/List_of_phylogenetic_tree_visualization_software

FigTree: in this tool, one can load several trees, one per line, and there is a 'Prev/Next' menu to iterate over them, but no apparent option to display all of them at once.

Phylo.io: in this tool, one can load several trees, and display one at a time or two of them, with colouring options that compare the common nodes/branches. There isn't an apparent option to visualize many trees in one go.

PhyBin: this tool (https://cs.indiana.edu/~rrnewton/projects/phybin/) seems to cluster the sets of trees, but not visualize all of them together.

GrapeTree: this tool () can visualize only one tree at a time. It displays the last tree in the list if more than one are copy+pasted in the input text box.

There is another group of tools under Graph drawing:
https://en.wikipedia.org/wiki/Graph_drawing

Cytoscape: this tool is designed for network visualizations, which could be an option. Unsure if it visualizes newick trees, but there is are tools to convert trees to networks that could potentially work in combination with Cytoscape, e.g. https://github.com/MonsantoCo/Newick2Network

I recommend using densitree: https://www.cs.auckland.ac.nz/~remco/DensiTree/

And rwty: https://github.com/danlwarren/RWTY

Densitree plots the superimposed set of trees in a heatmap-like format:

And rwty, a very powerful R package that allows visualisation of treespace in an NMDS format (among other methods):


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

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

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

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

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

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

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


Finding gene length using ensembl ID I want to find the length of a list of genes, of Homo sapiens, that is reported in the GEO database. I have gathered the ensembl id's of those genes.

I understand this information can be parsed from the start and end position given in GTF file available in ensembl.org.

I would like to ask for advice on how to programmatically download the gtf files with ensembl id's.

On the ensembl webpage a couple of options are listed, API and ftp access. I would like to know if sample scripts are available.

Update: I found a script here, it shows how to parse the data from gtf file. Could someone help in understanding how the code can be modified to obtain gene length from start and end position?

Edit 1:
By start and end position, I refer to the start and end values given in the "fields" section here and the file from which I want to parse the values is the second last file available in this link.
For instance, the end - start value of the gene with ensembl id ENSG00000131482,
gene :G6PC

Edit 2: I get the gene names and the expression values from the SOFT files documented in GEO, link

The following is a sample output from SOFT file. It gives me the probe_set(ID_REF) , Gene name and expression values from samples (GSMXXXXXX). Here we see fro G6PC, values are reported from two probe_sets. Can we get the transcript id's from the probe_sets?

Obtaining the transcript id will help me in finding the length of the specific transcript.

> Table(gds)[2449,0:5]
           ID_REF IDENTIFIER GSM388749 GSM388750 GSM388753
2449 1555612_s_at       G6PC   4.39472   4.52235   4.26869
> Table(gds)[16399,0:5]
         ID_REF IDENTIFIER GSM388749 GSM388750 GSM388753
16399 206952_at       G6PC   2.76693   3.45107    2.9815

If you really just want the start and end positions of all genes in a list of genes from Homo_sapiens.GRCh38.93.gtf.gz, you could do this in awk. For example, if your list of genes is:

$ cat genes
ENSG00000131482
ENSG00000223972
ENSG00000228794
ENSG00000187961

You can do:

$ awk '{
        if(NR==FNR){
            wantedGenes[$1]++
        }
        else{
            if($3=="gene"){
                gsub(/";*/,"",$10);
                if($10 in wantedGenes){
                    print $10,$4,$5
                }
            }
        }
      }' genes Homo_sapiens.GRCh38.93.gtf
ENSG00000223972 11869 14409
ENSG00000228794 825138 859446
ENSG00000187961 960587 965715
ENSG00000131482 42900797 42913369

Explanation:

if(NR==FNR){ wantedGenes[$1]++ } : NR is the input line number and FNR is the current line number of this input file. The two are only equal while reading the first file. So, if this is the first file, save the first field of each line (the gene name) in the array wantedGenes.
if($3=="gene"){ : if this line is a 'gene' entry.
gsub(/";*/,"",$10); : the gene names in the GTF file are quoted and have an ; after the quote. For example: "ENSG00000187961";. Since we don't have this in the list, we need to remove it so we can match the gene name. This gsub command will remove all " followed by 0 or more ; (the ;* means 0 or more occurrences of ;) from the 10th field of each line. So, it will replace "ENSG00000187961"; with ENSG00000187961.
if($10 in wantedGenes){print $10,$4,$5} : if the 10th field (the modified gene name from above) is in the wantedGenes array, then print the 10th field (the gene name), the 4th (start) and the 5th field (end).

You can easily modify this to print the length instead of the start and end positions:

$ awk '{
        if(NR==FNR){
            wantedGenes[$1]++
        }
        else{
            if($3=="gene"){
                gsub(/";*/,"",$10);
                if($10 in wantedGenes){
                    print $10,$5-$4
                }
            }
        }
      }' genes Homo_sapiens.GRCh38.93.gtf
ENSG00000223972 2540
ENSG00000228794 34308
ENSG00000187961 5128
ENSG00000131482 12572

As requested in chat, here's a python version of the same thing:

#!/usr/bin/env python3
import re
import sys

wantedGenes = []
with open(sys.argv[1], mode='r') as f:
    wantedGenes = f.readlines()

wantedGenes = [ x.rstrip() for x in wantedGenes]

print('# Gene\tstart\tend\tlength')
with open(sys.argv[2], mode='r') as f:
    for line in f:
        if not line.startswith('#'):
            fields = line.split('\t')
            annotations = fields[8].split(' ')
            geneName = re.sub('[";]', '', annotations[1])
            if geneName in wantedGenes and fields[2] == 'gene':
                print('%s\t%s\t%s\t%d' %
                      (geneName, fields[3],fields[4],
                       int( int(fields[4]) - int(fields[3]) )
                      ))

Save the file as foo.py and run with the file names as arguments:

python foo.py geneList.txt Homo_sapiens.GRCh38.93.gtf


Building small molecules de novo Can anyone recommend a protocol for building the 3D structure of a small molecule de novo, given only a 2D rendering? So far my idea has been to build the molecule using the fragment builder in PyMol, then make conformers using Frog2. But I have no idea how realistic the starting structure is from PyMol and how exhaustively the conformers are sampled in Frog2.
The downstream application of this would be to dock the ligand into a protein binding pocket. If anyone has a more rigorous protocol for making sure all the conformational space is fully sampled, it would be greatly appreciated!

Most universities have a ChemDraw site lincence, so you may want to check that for drawing a non–3D-embedded molecule. Alternatively there are many other programs. PyMOL is not really a good choice. Also there are online drawing tools everywhere, where you can download a SMILES (see below).
Chemical formula is the wrong way to go, unless you are "enumerating".
The chemical formula of glucose is C6H12O6. This is the same for all hexoses. What you want is a SMILES. This string gives the connectivity of the atoms. For D-glucose it's OC[C@H]1OC(O)[C@H](O)[C@@H](O)[C@@H]1O, the @ symbols are chiral centres.
In Cheminformatics making a list of all possible compound with n-carbons is a big deal, this is called enumeration and you likely don't have to reinvent wheel. Furthermore, there are datasets that are filtered by the Lipinski rule.

GDB datasets (if you don't care if you cannot buy it)
Zinc-DB  (if you care if you can buy it)
Enamine Real (if you care if you can buy it reasonably)


What is a proper way for random subsampling of metagenomic data? Let's say we have a metagenomic sample that is paired-end FASTQ files including 10,000,000 DNA reads collected using shotgun sequencing.
How would one make a random subsample of the mentioned metagenomic sample with for example 1,000,000 reads? I know about the seqtk function that could be used in such cases, but my question is more about whether it is meaningful to represent a metagenomic sample with a small subsample generated using this way or it has to satisfy certain criteria, e.g. in terms of the number of reads and diversity to be a valid representative of the metagenomic sample?
Knowing that some species might be very rare in the sample, would it ever be proper to do subsampling on metagenomic samples?

I am not NGS expert, I known ML. Essentially you are performing data augmentation which is essential for ML. What I question is whether a replicate of e6 from a population of e7 would really be sufficient for ML because you need sample sizes (replicates in this case) of >100000.

The answer is really simple you bootstrap the data (sample with replacement), to obtain your replicates. Numpy has a random module that will do this,

p.random.seed(345)
pop = np.random.randint(0,10000000 , size=1000000)

Use pop to select from fastaq. The other way in is to define each fastq sequence within pop then simply,

sample = np.random.choice(pop, size=1000000)

Or a combination of the two. If you can do ML you can do Python coding.


How to correctly call a VCF file using damaged DNA? (IonTorrent & FFPE) EDIT: I am updating this question to make it more specific to my issue.

For context - original question prior to edit:

How do I obtain a deamination metric when doing the variant calling using the IonTorrent variant caller, and secondly, how do I correct my called variants for deamination to ensure that these don't provide false positives for the downstream analysis?

Background

I'm calling variants using the IonTorrent TorrentSuite on DNA which has been sequenced from formalin-fixed paraffin-embedded (FFPE) tissue. This has a major issue in that without addition of uracil-N-glycosylase, some of the Ts in the original DNA are deaminated to Cs, which upon sequencing and calling variants can show up as mutations, either as T>C transitions or G>A (from the opposite strand, due to PCR in the library prep). I do not have any idea how long these samples were stored without UNG before sequencing.

TVC (Torrent Variant Caller) gives a deamination metric (essentially, sum of T>C and G>A variants over all variants called), and for our samples, the highest value seen is ~0.92. Naively postprocessing the variants show that for these samples, C>T/T>C transitions overwhelm the remaining variants among my samples.

Application

My use case is this: these are medical samples, which have been inspected by a pathologist (hence the FFPE treatment), and I want to determine which variants are predictive of outcome, hence I have two potentially contradictory goals: reduce false positives and capture the rarer variants which may hold predictive power.

Question

Given the IonTorrent variant calling pipeline (sequencing > BAM file > TVC > VCF file with deamination statistic):

(Q1) Is there a way of correcting the output VCF to remove transition and deamination errors?
(Q2) Alternatively, is there a set of filters to use in bcftools to reduce the effect of the errors on how variants are called?

Sites that were deaminated are often called abasic sites.

Q1: How to obtain DNA damage metrics?

You might want to check out mapDamage - this software is used to perform these sorts of analyses and produce plots of types of DNA damage along sequencing reads.

Q2: How to correct for abasic sites and other DNA damage errors?

You can use your choice of kmer-based or alignment based error correction software. The errors will be random and if the sequencing depth is sufficient the erroneous kmers at the ends of reads will be corrected. However, depending on what you are trying to do downstream with your reads you may be better served by not correcting the reads, then dealing with transitions/damage later.


How to generate Genbank format out of paired end fastq data? I've unmapped cleaned paired end sequence data in fastq format of a bacterial genome. I want to get a sequence data in Genbank format in the end. What are the exact steps that I have to follow in galaxy webserver and if in some point necessary, also using linux based apps? I would appreciate a lot of details.

The steps needed to assemble a genome in Galaxy are layed out in the Assembly training material from the Galaxy Training network. The software for this should be available on the main galaxy servers (e.g., usegalaxy.eu). Since you have a prokaryote, I suggest you follow the "Unicycler Assembly" workflow, which is also available on usegalaxy.eu.


Biopython code refactoring I want to make function get_abstract, but i don't know what this can look like if the function must return idlist.

from Bio import Entrez, Medline, SeqIO

list_of_ids = []
Entrez.email = 'ski89@gmail.com'

handle_one = Entrez.esearch(db="pubmed", term='e coli',
                           retmode="text", retmax=100000)
record = Entrez.read(handle_one)
idlist = record["IdList"]
for ID in idlist:
    list_of_ids.append(ID)

    if len(list_of_ids) == 100000:
        handle_one = Entrez.esearch(db="pubmed", term='e coli',
                           retmode="text", retstart=100000, retmax=100000)
        record = Entrez.read(handle_one)
        idlist = record["IdList"]
        for ID in idlist:
            list_of_ids.append(ID)

            if len(list_of_ids) == 200000:
                handle_one = Entrez.esearch(db="pubmed", term='e coli',
                           retmode="text", retstart=200000, retmax=100000)
                record = Entrez.read(handle_one)
                idlist = record["IdList"]
                for ID in idlist:
                    list_of_ids.append(ID)

so the args will be the start and TERM:

get_abstracts(start,TERM):
handle_one = Entrez.esearch(db="pubmed", term=TERM,
                               retmode="text", retstart=start, retmax=100000)
                    record = Entrez.read(handle_one)
                    idlist = record["IdList"]
                    for ID in idlist:
                        list_of_ids.append(ID)

Can someone help refactor this? cheers

Here is a sample of code that should get you started.

I created a separate get_abstracts function.  Then that function get called repetitively until less than 100K ids are returned.  If you wanted to get fancy this whole process could be called as a recursive function but that seemed overly complex.

BUFFER_SIZE = 100000

def get_abstracts(start, term):
    """Get N number of records from Entrez
    :returns list: List of ids returned from search
    """
    handle_one = Entrez.esearch(db="pubmed", term=term,
                       retmode="text", retstart=start, retmax=BUFFER_SIZE)
    record = Entrez.read(handle_one)
    idlist = record["IdList"]
    return idlist

list_of_ids = []
start_at = 0

# Iterate the search returning 100000
while true:
    ids = get_abstracts(start_at, 'e coli')

    list_of_ids.extend(ids)

    # Break if less than 100K results were returned
    if len(ids) < BUFFER_SIZE:
        break
    start_at += BUFFER_SIZE

# Do something with list_of_ids


Using signal peptide and the expression levels of signal recognition particle in secretome analysis I have not found any work which investigates assessment of differences in levels of secreted proteins by taking advantage of differential expression of the genes which mediate the secretory pathway.

For example, suppose the gene whose product binds specifically to a certain signal peptide is down-regulated. I would assume that this implies a down-regulation in all the levels of secreted proteins to which it binds.

Is such an analysis possible to date? If so, where?

Thanks in advance.

This is something that Qiagen's Ingenuity Pathway Analysis is meant to solve. Unfortunately, the process of solving pathway dependency problems like these is difficult, and frequently leads to the discovery of super-linked proteins that are seen everywhere (and consequently fairly uninteresting from a biological perspective). The challenge is in filtering out the useless noise of "statistical significance", and identifying the biologically-significant signal. That's not something that computers are particularly good at doing at the moment.

My own suggestion is that you would be better off finding someone who has domain knowledge in the area of secreted proteins, and trying to present results to them in such a way that it makes it easy to interpret the results and ask further questions. If you are that person, perhaps there's a more specific question that you could ask to further investigate the problem that you have.


Combine structures in the same session - PyMOL I'm a student currently looking at antibody responses against a viral target protein of interest.

I have my own, annotated PyMOL session of my protein and I also have .pdbs of crystallised antibody structures that bind to this protein and rougly know whereabouts it binds. I would like to know if it is possible to combine these two structures into a single session to observe or at least overlay the two structures together. Many thanks!

Open the PyMOL session.
Use the load command to load your PDB file(s): load file.pdb.
Use the align or super commands to overlay structures, e.g. align prot1, prot2.


1000 genomes missing SNPs on chr X in GRCh38 I am trying to get a reference together for imputation using the 1000 genomes data set. The GRCh38 version of the data set was made available earlier this year, but seems to be missing a significant number of SNPs on the X chromosome (~96%).

Is there any resource that can explain why? I have tried searching on their home page with no luck.

The data sets in question are GRCh38, GRCh37.

As I told you on 1000 genomes helpdesk:
Because of the ploidy settings we used for our variant calling, we only called variants in the the pseudoautosomal region on the X chromosome. We are working on settings to correctly analyse X variants in both males and females and plan to analyse X properly soon.


Programmatically retrieve graph features from Cytoscape Is there a way to programmatically retrieve graph features (a networkx graph) from Cytoscape in Python 3?

EDIT: To add more details, I had a dictionary of proteins, each protein was associated with a list of proteins that it interacts with. I wanted to create a protein-protein interaction graph so that I could extract graph features from it, including things like k-core value of each node.

I was able to extract the graph features that I wanted just using networkx, without going through the intermediate step of Cytoscape.

EDIT: Initially, I planned to import the data into Cytoscape and extract graph features from Cytoscape. However, I noticed that a much easier solution would be to just use the methods of the networkx library. I initially had a few graph features in mind that I wanted to extract like k-core value of each node, and I was able to do this using methods from networkX. To calculate k-core for instance, I used the k_core method. For more information on networkX methods, check the docs


What sequencing data metrics should I record? I'm new to bioinformatics and Im starting a new microbial sequencing project and I'd like to record all of the qc data correctly. This is research and I'm not sure what I'll need later. Is there any consensus on useful qc metrics to record in a database for each sequencing run?

No rnaseq, running FastQC, assembling with a defined reference for a single species.

I would say that the following metrics are useful:

Sample source (Extraction place, person, date)
Details about pre-HTS preparation (e.g. amplicon PCR, rRNA removal)
Reads per sample (+ per barcode / lane, if available)
Sequencing kit (e.g. TruSeq v3), useful for adapter trimming
Sequence length (e.g. 125bp)
Paired end / single end
If paired end, target fragment length in sample prep (e.g. 650bp)
Reference sequence metadata

Scientific name
Location of reference genome (preferably a stable web address)
Version of reference genome
Read mapping rate (per sample / barcode / lane), as a percentage of total reads
Mapping program + version

If available, potential contaminating organisms

Name / description
Location of contaminant sequence (preferably a stable web address)
Estimated mapping rate (bearing in mind that there may be some overlap with the target)
Mapping program + version

More generally, ENA / EBI has an excellent reference in sample checklists, which indicate the minimum acceptable metadata for attaching to sequence archive data:

https://www.ebi.ac.uk/ena/submit/checklists


Why is it necessary to add hydrogen and delete water before protein-ligand docking? What is the reason for adding hydrogen and removing unnecessary water molecules from the protein structure before protein-ligand docking? FYI, the tools I used for docking is GOLD.

My answer is mainly based on molecular dynamics,but in this context it shouldn't make much of a difference.

The main reason has to do with the force-field-water model used during the simulation. Examples of water models are SPC, TIP3P or TIP4P; same molecule (water) but different topology.

A pdb usually comes from different sources and it can contain molecules that are not standard and therefore not recognised by the model of your simulation. Usually, once the structure is loaded, your system is re-built so that it conforms to the model topology. There are exception though, for example if your water is part of a functional active site. In that case you can keep it, but its structures (topology of the water molecule) has to be recognised by the model.

Hydrogens are usually added to "saturate" and stabilise the bonds.


Selecting part of an extracted ligand I have the PDB 6LU7.

I extracted the ligand using extract lig, org  and then the protein using extract prot, poly.  The ligand has three parts:  02J, PJE, and 010.

I want to select the 02J and 010 IN A SCRIPT and delete them.

I know how to do it by simply clicking, but want to incluse this in a script.

I know that structure way too well —SARS-COV-2 MPro.
In PyMOL command line
 remove resn 02J or 010

or
 cmd.remove('resn 02J or 010')

In python
 import pymol2
 with pymol2.PyMOL() as pymol:
      pymol.cmd.fetch('6LU7')
      pymol.cmd.remove('resn 02J or 010')
      pymol.cmd.save('6LU7.mod.pdb')
      #or pdbblock = pymol.cmd.get_pdbstr()

From the sounds of it you'd like to make the peptidomimetic a peptide. If you wanted to do the horrid thing of making PJE, the capped side chain glutamine ketamide, into glutamine, the command to rename atoms is a multitude of cmd.alter('resn PJE and name <atomname>', 'name="<newname>") followed by a cmd.alter('resn PJE', 'resn="GLN"') and lastly cmd.sort() to commit the changes. Not pretty, but is all in PyMOL.


What is "Possible polymerase run-on fragment"? I used gffcompare recently to compare a newly assembled transcriptome and a published transcriptome. The class code "p" transcripts are "possible polymerase run on".

What does that mean? I tried very hard but still cannot find the answer.

On most coding genes (with the exception of replication dependent histone genes), Transcripts (as opposed to transcription) are terminated by cleavage polyadenylation sites, where the growing transcript is cleaved by a protein complex and the poly-A tail. However, PolII continues to transcribe until it is terminated by torpedo termination (see here). Normally this extra RNA after the polyA site is degraded by Xrn2, but if you capture any particular snapshot of the RNA in the cell at any one time there will be a small amount that corresponds to this run on transcription.


Delly finds identical deletions with different consensus sequences I am using delly to call structural variants on a human WGS sample. I noticed that it finds certain variants that are essentially duplicates. For example:

chr4    177594033       DEL00005339     G       <DEL>   .       PASS    PRECISE;SVTYPE=DEL;SVMETHOD=EMBL.DELLYv0.8.1;CHR2=chr4;END=177594095;PE=0;MAPQ=0;CT=3to5;CIPOS=-3,3;CIEND=-3,3;SRMAPQ=60;INSLEN=0;HOMLEN=3;SR=5;SRQ=0.993243;CONSENSUS=GACAAAAACATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAAC;CE=1.75466     GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-51.3797,0,-42.791:10000:PASS:1115:1821:1820:1:0:0:23:27
chr4    177594033       DEL00005340     G       <DEL>   .       PASS    PRECISE;SVTYPE=DEL;SVMETHOD=EMBL.DELLYv0.8.1;CHR2=chr4;END=177594095;PE=0;MAPQ=0;CT=3to5;CIPOS=-3,3;CIEND=-3,3;SRMAPQ=60;INSLEN=0;HOMLEN=3;SR=12;SRQ=0.994186;CONSENSUS=ATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAACAATGTAATAAGTATGTTTGAATAGTATTATGAT;CE=1.72637    GT:GL:GQ:FT:RCL:RC:RCR:CN:DR:DV:RR:RV   0/1:-51.3797,0,-42.791:10000:PASS:1115:1821:1820:1:0:0:23:27

The two variants above have the same start position and the same end position (END=177594095). However, they differ in the SR, SRQ and CONSENSUS fields. These are defined in the VCF header as:

##INFO=<ID=SR,Number=1,Type=Integer,Description="Split-read support">
##INFO=<ID=SRQ,Number=1,Type=Float,Description="Split-read consensus alignment quality">
##INFO=<ID=CONSENSUS,Number=1,Type=String,Description="Split-read consensus sequence">

Of course, the consensus sequences are quite similar, but they aren't identical (I have simply added some spaces to the beginning of the 2nd sequence here to align them):

GACAAAAACATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAAC
         ATATAAACAAATACAAATTATTATTATATGAAGGGTGAAAAATTATGTTCAGCTTTAATATTATGTGGCCTGAAACAAATTTGTATAAGTCATATGATCCGTTTTGTATATTCATTAAAATAGTCATAAAAATATAAACAATGTAATAAGTATGTTTGAATAGTATTATGAT

So what am I looking at here, exactly? If this is a deletion on chromosome 4 from position 177594033 to position 177594095, then how can the length of the consensus sequence be different? This is the "split read consensus", so presumably the sequence that appears to have been deleted based on how the reads covering that region are split. This means that these reads align to the genome with a gap and this gap indicates the deletion (image source):

Therefore, if the start and end position are the same, the consensus should also be identical. And, in any case, if the start and end positions are the same, then this should be the same structural variant, so why is it being reported twice? How can I interpret this result?

It turns out that this was just a bug which has been fixed by this commit. The current master version of delly no longer exhibits this behavior.

If you are facing this issue, either build the current master version or wait a bit for the next release which should also include this fix.


Downloading a reference Genome for Bowtie2 How do I download a reference genome that I can use with bowtie2? Specifically HG19. On UCSC there are a lot of file options.

It’s a matter of preference I guess but I recommend the Ensembl builds. Decide whether you want the toplevel or primary assembly, and whether you want soft-masked, repeat-masked or unmasked files. The naming schema is very straightforward; the combinations are described in the README file, and all files reside in one directory.

For example, if you want the unmasked primary assembly, the file to download would be Homo_sapiens.GRCh37.75.dna.primary_assembly.fa.gz.

As for GoldenPath/UCSC, there’s no need to download and concatenate separate chromosomes (contrary to what the other answer said); you can download the whole (toplevel) reference from the bigZips directory; from the README:

  This directory contains the Feb. 2009 assembly of the human genome (hg19,
  GRCh37 Genome Reference Consortium Human Reference 37 (GCA_000001405.1)),
  as well as repeat annotations and GenBank sequences.

There are essentially three options here:

chromFa.tar.gz, which contains the whole genome in one chromosome per file;
chromFaMasked.tar.gz, the same with repeats masked by N;
hg19.2bit, which is the whole genome in one file, but needs to be extracted using the utility program twoBitToFa, which needs to be downloaded separately.

In any case, I always download the reference and build my own index for mapping, since this allows me more control; not everybody might need this much control, but then building the index once is fairly fast anyway.


Storing FindAllMarkers results in Seurat object I am currently working on multiple datasets where each is managed by a separate Seurat object. For each I'd like to also compute the marker genes using the FindAllMarkers function. Where would it be a good practice to save the markers returned from the FindAllMarkers function? I tried adding an extra slot to the Seurat object with no success.

I saw the calc.params slot and though I could add an item to that list. However, it doesn't seem like the right place as that slot is a list of the parameters used for the different functions executed on the data.

Currently I have some wrapping object that has a slot for the Seurat object and a slot for the results from the FindAllMarkers. Is there any better place to put it? Somewhere within the Seurat object?

You can stash anything you like in the misc slot (present in both v2 and v3 Seurat objects). I often use it for storing marker information to help with organization.


Plotting metadata from GenomicRanges object I have a GRanges object with metadata that I want to plot. Many of the online tutorials which manipulate these GRanges for plotting show how to plot the main columns (chromosome,start/end,strand) but nothing to plot the metadata as a function of the position in each row.

How precisely should I do this? Ideally, I would like to plot metadata from multiple GRanges objects on the same plot to compare them.

EDIT (following comment):

GRanges object with 6 ranges and 140 metadata columns:
      seqnames        ranges strand |
         <Rle>     <IRanges>  <Rle> |
  [1]     chr1 100000-199999      * |
  [2]     chr1 400000-499999      * |
  [3]     chr1 600000-699999      * |
  [4]     chr1 700000-799999      * |
  [5]     chr1 800000-899999      * |
  [6]     chr1 900000-999999      * |
                  eigen   domain         arm
              <numeric> <factor> <character>
  [1] -1.62475825405065     open          1p
  [2] -1.62475825405065     open          1p
  [3] -1.62475825405065     open          1p
  [4] -1.62475825405065     open          1p
  [5] -1.63354915136132     open          1p
  [6]  -1.6322733358426     open          1p
      exon.fraction       gc     short      long
          <numeric> <matrix> <numeric> <numeric>
  [1]       0.02968  0.45815         0         0
  [2]       0.04305  0.40192         0         0
  [3]       0.00959  0.42091         0         0
  [4]       0.04593  0.43337         0         1
  [5]       0.02425  0.44822         2         1
  [6]       0.16977  0.62123        27        56
         nfrags             ratio   short.corrected
      <numeric>         <numeric>         <numeric>
  [1]         0               NaN              <NA>
  [2]         0               NaN              <NA>
  [3]         0               NaN              <NA>
  [4]         1                 0 -3.50822623240367
  [5]         3                 2  2.81796238367926
  [6]        83 0.482142857142857  39.8824095079959
...

My question regards plotting a column like eigen as a function of the position (given by seqnames and ranges)

You can access the data.frame just using $. All of them can be accessed using values() . Something like this:
library(GenomicRanges)
gr = GRanges(seqnames=1,IRanges(start=1:10,width=1))
values(gr) = data.frame(eigen=rnorm(10),ratio=runif(10))

plot(start(gr),gr$eigen)

or
plot(start(gr),values(gr)$eigen)


Comparing two genome annotations I have one-dimentional array (human genome). Also I have two annotations for it, we can think about them as different peaks (it's nucleosome and secondary structures). How can we find correlation and/or causation between this annotations? Currently we are forming hypothesis like "one element is always before other" and statistically testing them. Can we be smarter and use known tools or at least automate this?

You can express the annotations of the genome with a dot plot, x and y represent the two annotations along the genome. The you can summary the result with Spearman correlation coefficient(s) based on whole set of annotations or particular gene families.


How do i map subcellular localizations to lower the number of classes to predict For a project, I want to predict the subcellular localization of some proteins. I am a student in software Informatics at the University of applied Science HEIG-VD.

Subcellular localization prediction is a common problem in the bioinformatics, but I have a problem

With the data I found. In this database: of human and mouse proteins, we can see that there are more than 30 different localizations.

At my level, and for the use of my exercise, 30 classes to predict is too much.

Where can I find a mapping to reduce those 30 classes to 4 or 5?

I tried QuickGO, to see the ancestor tree of the different classes, but nothing conclusive. Somebody might have done that, but I lack knowledge on bioinformatic web research.

There are the classes I wanted to map to reduce the number of classes :

Cytoplasm (GO:0005737)
Cytoplasmic Vesicles (GO:0016023)
Endoplasmic Reticulum (GO:0005783)
Endosomes (GO:0005768)
Extracellular (GO:0005576)
Golgi Apparatus (GO:0005794)
Lipid Particles (GO:0005811)
Lysosomes (GO:0005764)
Melanosome (GO:0042470)
Mitochondria (GO:0005739)
Nucleus (GO:0005634)
Peroxisome (GO:0005777)
Plasma Membrane (GO:0005886)
Synaptic Vesicles (GO:0008021)
Cellular Component Unknown (GO:0008372)
Apical Plasma Membrane (GO:0016324)
Basolateral Plasma Membrane (GO:0016323)
Centrosome (GO:0005813)
Golgi Cis Cisterna (GO:0000137)
Cytoskeleton (GO:0005856)
Early Endosomes (GO:0005769)
ERGIC (GO:0005793)
Inner Mitochondrial Membrane (GO:0005743)
Late Endosomes (GO:0005770)
Medial-Golgi (GO:0005797)
Outer Mitochondrial Membrane (GO:0005741)
Secretory Granule (GO:0030141)
Golgi Trans Cisterna (GO:0000137)
Golgi Trans Face (GO:0005802)
Tight Junction (GO:0005923)
Transport Vesicle (GO:0030133)

One way to do this would be to plot the first 10 GO categories of the the CC sub ontology (CC from cellular component). Or you can just visit the official webpage to see the list of children terms of them.

However there aren't only 4 or 5 categories, so you could arbitrary pick some of those.


too long header for fasta file I have a fasta file, like this:

>TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3::g.3039::m.3039 TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3::g.3039  ORF type:complete len:100 (-) TRINITY_DN100000_c1_g1_i3:1027-1326(-)
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6::TRINITY_DN100000_c1_g6_i2::g.84365::m.84365 TRINITY_DN100000_c1_g6::TRINITY_DN100000_c1_g6_i2::g.84365  ORF type:complete len:112 (-) TRINITY_DN100000_c1_g6_i2:379-714(-)
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

if I want to only keep the second header, like this (the amino acid sequence keep unchanged). which command could i USE?

>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

I know a command

> sed 's/::.*//' input > output

can only keep the first header >TRINITY_DN100000_c1_g1 if I want to keep the second header with the isoform information >TRINITY_DN100000_c1_g1_i3, how should I correct this command?

Here's one way:

$ sed -E 's/>[^:]*::([^:]*)::.*/>\1/' file.fa
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

The idea is to match a >, then as many non: as possible ([^:]*), and then capture the string between the first :: and the second ::. Finally, you replace the entire line with the captured string.

Once you're satisfied this does what you need, you can edit the original file with -i:

sed -i -E 's/>[^:]*::([^:]*)::.*/>\1/' file.fa

Alternatively, set awk's field separator to :: and print the second field only on lines that start with >:

$ awk -F'::' '/^>/ ? $0=">"$2 : 1' file.fa
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*


How to download the whole directory of an ensembl FTP page? I am interested in downloading all the files at:

ftp://ftp.ensembl.org/pub/release-91/embl/bos_taurus/

I know how to download a single file from that page using wget (Unix) or download.file() (R) but is there an easy way of downloading all of them?

You can use the mirroring function of wget:

wget -m ftp://ftp.ensembl.org/pub/release-91/embl/bos_taurus/


What is the best way to programmatically convert Ensembl ids from one release to other? e.g. transcript ids from release 75 to 100 3 solutions that I had to unfortunately rule out.
Potential solution #1
As far as I know, the conventional way is to do this is by using CCDS ids (e.g. using biomart). However, the issue is that, as I see on the CCDS project page, the last release of CCDS ids was made almost 3 years ago (!) and it seems to be compatible with Ensembl release 92. Then how should one convert ids to the latest Ensembl release e.g. 100?

CCDS Release
Species
NCBI Annotation Release
Ensembl Annotation Release
Assembly Name
Assembly ID
Made Public

22
Homo sapiens
109
92
GRCh38.p12
GCF_000001405.38
06/14/2018

Potential solution #2
I am also aware of Ensembl's ID History Converter. However, I am looking for a programmatic approach. Also, I noticed that if I try to convert ~>20K ids with ID History Converter, in my experience, there is always some issue. e.g. an error I got when I tried to convert ~70K ensembl transcript ids.

As a side note, I wonder why this ID conversion functionality is not available via Ensembl's REST API. To me, id conversion seems to be a basic enough task that almost every bioinformatician has to do. So it would be immensely helpful to be someday able to do this programmatically using Ensembl's mighty a REST API.
Potential solution #3
Finally, as another potential solution, over at Ensembl's FTP website, located within the directory of the release to which I am interested in converting my ids to (release 100, target release), I saw that according to a README file, "ensembl_archive_" directory should contain "Data on historical Ensembl releases". So, I wondered if this directory contains a map between ids from different releases. However, without any extra info about column names, contents, file formats etc, I couldn't really figure out how to use those files. So, I had to rule this one out too.
Index of /pub/release-100/mysql/ensembl_archive_100/
../
CHECKSUMS                                          26-Mar-2020 14:31                 132
ens_release.txt.gz                                 26-Mar-2020 14:31                1098
ensembl_archive_100.sql.gz                         26-Mar-2020 14:31                 943
release_species.txt.gz                             26-Mar-2020 14:31               39620
species.txt.gz                                     26-Mar-2020 14:31                6669

Of note, I saw a similar "mysql" directory over at the directory of the release whose ids I want to convert (source release): http://ftp.ensembl.org/pub/grch37/release-100/mysql/homo_sapiens_core_100_37/ . However, again, I do not understand its contents without any extra information.
So, finally, I am not sure what is the best way to programmatically convert Ensembl ids from between different releases. As a test case scenario, I want to convert all the protein_coding transcript ids from GRCh37 (release 75) to GRCh38 (release 100). If anybody has any suggestions about this, please kindly let me know.

The ID history converter is available as a Perl script.
This script accesses these MySQL tables, so you could query these tables directly.
Note that the mapping session is only one release at a time, so you would need to run 25 queries to go from release 75 to 100.


Simulate and test CNV workflow? I'd like to evaluate a CNV project.

My aim is to evaluate if the scripts are sufficient for calling reasonable CNVs. I know they have a paper, but their scripts may be buggy... and all papers are biased...

While there are CNV simulators on Google search, none seem popular in the community. How should I approach this problem?

Simulate CNV variants from a genome
Simple evaluation on the called variants against my simulation?

I cannot tell exactly from the post or the linked project what type of data these scripts are designed for, but I am going to assume that "Scripts involved in our workflow for detecting CNVs from WGS data using read depth-based methods" refers to short-read Illumina sequencing.

Simulations can be helpful, but they are usually imperfect. Nowadays we have long-read data available. I would suggest obtaining some long-read data for which short-read data is also available, doing CNV detection on both datasets independently, and comparing the results.

The problem with using long-read data is that any differences you find will be conflated with the biases of the technology. So, if you find a difference, then this approach will not tell you if the differences are due to the biases of the technologies, or due to bugs in the scripts. If you use simulations, then this conflation of biases remains.

Short of re-writing the pipeline, the only other suggestion I have is to borrow from software testing. Try constructing a few minimal examples that might also test some corner cases and see if the scripts behave exactly as expected.


What is a quick way to find the reverse complement in bash I have a DNA sequence for which I would like to quickly find the reverse complement. Is there a quick way of doing this on the bash command line using only GNU tools?

Thanks to Manu Tamminen for this solution:

echo ACCTTGAAA | tr ACGTacgt TGCAtgca | rev


Produce .vcf file of ALL mutations in .bam file I got tasked with a problem, that i thought would be quite simple to solve, but turned out to be quite tricky. Our lab is running targeted mutagenesis experiments in yeast using crispr base editors. What we have done is set up an experiment where the base editor selectively mutates a region inside of a gene of interest in a culture of yeast cells. We then extract the DNA of the whole culture and perform sequencing of our gene of interest (the size of which is around 1kb).
What i would like to do is plot the mutation frequency of say, C -> G edits, across the length of the entire gene (around 300bp), with the hope that i see a spike in mutations at the site where the base editor binds. Or at least, higher mutation rates at this site compared to the background mutation rate.
These edits are quite rare (we think), and since it's in a culture of multiple cells, they don't occur at the exact same position on the gene.
I initially thought i could do this by quality filtering the reads, aligning them to the gene sequence in order to produce a .bam file, and then piping this file through variant calling algorithms such as vcftools. However, i noticed that i get very few variants out of this, even when using extremely relaxed settings (e.g. p-value < 1 for the variant call).
Im therefore wondering if anyone knows how to extract from a .bam file ALL the mutations across all reads that do not match the reference sequence?
I'm aware that this will probably lead to quite a huge .vcf file, but since the reference sequence is only 300bp, it should still be manageable? It could also be the the experiment does not work, but i would at least expect some mutations in the sequencing reads compared to the reference gene sequence.
Any help would be greatly appreciated!

Ensure that the variant caller stringency is relaxed all the way. For example in freebayes use --min-alternate-count 1 --min-alternate-fraction 0 .


Phyre2 vs ITasser, completely different models generated Does anyone have experience generating pdb structures with Phyre and ITasser online tools. The results generated from each given the same amino acid sequence input are very different and I am wondering whether or not this is a usual experience. I know ITasser was the number 1 rated in the CASP trails, but should the results really be this disparate?

I'm less familiar with Phyre, but I-TASSER is a really sophisticated system that takes the results of a search using multiple threaders and plugs them into an ab initio simulation which tries to minimize the energy of the models by sampling many possible 3D conformations, which I don't think Phyre does.

https://en.wikipedia.org/wiki/I-TASSER#/media/File:I-TASSER-pipeline.jpg

Compare with a similar workflow schematic for Phyre:

http://www.nature.com/nprot/journal/v10/n6/images/nprot.2015.053-F1.jpg
http://www.nature.com/nprot/journal/v10/n6/fig_tab/nprot.2015.053_F2.html

Structure prediction still has a long way to go, and you'll always get better results if there are close homologues available in the PDB, but given the consistent high performance of I-TASSER in CASP I would treat those results as more significant. That said, it can't hurt to consider multiple answers.

Edit: included blmoore's link to second half of Phyre protocol schematic


Is it possible to filter contaminated reads for raw PacBio sequences (not HiFi reads) before assembly? De novo genome assembly for non-model organisms face the issue of bacterial contamination. For assembled contigs with mostly bacterial-like sequences (based on BLAST search), the entire contig can be dropped off. But for a long contig with small fragments of bacterial-like sequences, one way to solve it is to break the contigs into smaller pieces and drop the bacterial-like pieces. It would work better if the raw reads could be removed BEFORE assembly. This was commonly done for Illumina reads, but for PacBio long reads (non-HiFi reads) with a higher error rate, is it possible to filter contaminated bacterial sequences based on sequence similarity, like BLAST, Kraken? Of course, the potentially contaminated region should be long enough to drop the entire sequence, for PacBio long reads. Maybe 30%? Any suggestion of the processing details?

I don't think that in principle the issue is any different for PacBio and Illumina; it's just that the process might be noisier with higher error rates. Note that PacBio errors aren't randomly distributed substitutions and probably many (most?) short-ish k-mers (<30nt) will likely be intact. I would suggest doing k-mer filtering of your reads as suggested here or here or here.
Alternately, use your favorite alignment-based technique, if you have a good idea of what your contaminant is. In principle PB reads will be better filtered by this method than Illumina as local alignment will be able to figure out the errors quite well with longer queries.
If you go the route of trying to do alignment to some set of contaminant references, the obvious answer would be to use minimap2 as your aligner as this is exactly the task that it is designed for. You can then use samtools to filter out anything that maps well.
Here is a related question/answer that also addresses the contaminant issue and suggests using Blobology, which I'm not familiar with but may be helpful (granted it also seems to be Illumina-centric). Genome-size estimation using k-mers on PB reads appears to be pretty accurate too, so I think that it's likely ok to work with raw PB k-mers.


Reads count in metagenomics Background:  I am developing a pipeline for metagenomic studies of human gut microbiote. In particular, I am mapping the reads data originating from shotgun whole genome sequencing onto a gene catalogue (similar to IGC) and counting reads mapped on each feature/gene. Till now I have been using NGLess, which is a rather convenient tool. However, we have some questions about how exactly it counts reads and it might be not the best option for scaling up our pipeline.
Question:
I am looking for alternative tools for gene count. A brief search has shown that htseq-count and featureCount are two frequently used tools. However they mainly appear in RNAseq context, and I am not certain whether they could be used in my case. I would also appreciate comments on the possibility of using samtools mpileup.

I think that there is a good argument for using technologies adapted from RNA-seq for metagenomic abundance quantification. The tool that I am most familiar with (that I use for this purpose) is kallisto, which uses pseudoalignment instead and is therefore less compute intensive.
For some background on the conceptual linkages between metagenomic abundance and transcriptomic abundance you can see the kallisto metagenome abundance paper or the explanatory blog post.
I think that it is preferable to use a tool that is specifically adapted for metagenomics, as there are always some numerical issues that might have been accounted for by the authors. I think that for this purpose it is maybe best not to naively use tools like mpileup, though I have no specific objection to it. However, if you look around you can see various papers using e.g. Cufflinks and Cuffdiff for effectively the same purpose, and Cufflinks is definitely intended for RNA-seq.
Depending on your specific application, you might want to use a different tool like mOTUS2, if you are interested in estimating abundance of specific genomes or lineages in a metagenome mixture. This method is probably more accurate as it is based on averaging across the whole genome, but it is not very useful if what you are interested in is more e.g. functional profiling of the metagenome.


Expand granges object different amounts upstream vs. downstream I am attempting to get gene regions and their immediate neighborhoods using Bioconductor GenomicRanges packages.

It is very easy for me to extend the regions around genes if I want to extend it the same distance upstream and downstream as in this example:

library(GenomicRanges)
library(Homo.sapiens)

genes = genes(Homo.sapiens)
genes_5kb_flank = genes + 5000

And if I want to extend both upstream and downstream different amounts centered on the TSS I can use promoters(genes, upstream = 5000, downstream= 1000)

But I am having trouble finding a function that is strand-aware and also chromosome length aware that allows me to extend the gene region a different amount upstream vs. downstream... For instance if I wanted a 5kb upstream flanking region and a 1kb downstream flanking region as is used for genomic region annotations on GREAT.

Anyone know of a drop in function for this?  something like promoters() but for the full gene region instead of only around the TSS?

There doesn't appear to be a built in function for that, so here's a modified version of promoters():

expandRange = function(x, upstream=2000, downstream=1000) {
  strand_is_minus = strand(x) == "-"
  on_plus = which(!strand_is_minus)
  on_minus = which(strand_is_minus)
  start(x)[on_plus] = start(x)[on_plus] - upstream
  start(x)[on_minus] = start(x)[on_minus] - downstream
  end(x)[on_plus] = end(x)[on_plus] + downstream
  end(x)[on_minus] = end(x)[on_minus] + upstream
  x
}

You can then call expandedGenes = expandRange(genes, 5000, 1000).


Pepsin digest (cleavage) does not work using RE? Aim is to code the theoretical peptic cleavage of protein sequences in Python. The cleavage (cutting) rule for pepsin is: 1234^567. This is a position, ^ stands for the cleavage point.

There are rules that it mustn't have at all and rules where it must have at least one. The rule for mustn't is in position 2, it mustn't be R, K or H. In position 3 and 7, there can't be a P and in position 4, there mustn't be an R.

If it follows "the must not" rules, then it will check "the must" rules so that if any of one them is met, it will cleave. So if there is a P in position 1 or 2, or if there's an L, W, Y or F in position 4 or 5 then, it will be cleaved (use newline).

This is the code I have so far, but it doesn't work yet:

import re
print('INSERT YOUR PEPTIDE')
X= input();
r=re.findall('(?<=[P])(?<=[P,^R,^K,^H])(?<=([^P])(?<=[L,W,F,Y,^R])(?<=[L,W,Y,F])/D(?<=[^P])',X)
print(r)

My input is:

TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI

How can I make this work?

I am really not sure if I've understood your rules correctly, but I think you want this:

If all of these are true:

position 2: not R, K or H
position 3: not P
position 4: not R
position 7: not P

And, if at least one of these is true:

position 1: P
position 1: P
position 4: L or W or Y or F
position 5: L or W or Y or F

Then cleave.

If so, then you probably don't want to do it using regular expressions. First because it will be hard, and second because it will be absolutely awful to maintain if you need to come back to it. Instead, try something like this:

#!/usr/bin/env python3
# X = input();
X = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
seq = list(X)
for i, res in enumerate(seq):
    if i > len(seq) - 7:
        break
    if (
            seq[i+1] in ('R','K','H') or
            seq[i+2] == 'P' or
            seq[i+3] == 'R' or
            seq[i+6] == 'P'
    ):
        continue
    if res == 'P':
        start = i
    elif seq[i+1] == 'P':
        start = i+1
    else:
        continue
    if (
            seq[i+4] in ('L', 'W', 'Y', 'F') or
            seq[i+5] in ('L', 'W', 'Y', 'F')
    ):
        print("%s\n%s\n%s" %
            ("".join(seq[0:start]),
             "".join(seq[start:start+7]),
             "".join(seq[start+7:])
            )
        )

This will produce the following output when run:

$ foo.py
TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMV
PKCFNEL
LMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI

If you insist on using regular expressions, you could separate it into two expressions and use something like this:

import re
X = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
re1='P[^RKH][^P][^R][LWYF].[^P]'
re2='P[^P][^R][LWYF].[^P]'
r=re.findall('(?:P[^RKH][^P][^R][LWYF].[^P])|(?:P[^P][^R][LWYF].[^P])',X)
print(r)

Before you do, however, please consider the wise words of XKCD 208:


EMT/EMT-like processes in bioinformatics to study cancer progression First I am not sure why, but it seems EMT can be used to infer the status of cancer progression.

What kind of databases are there for the EMT process in the context of cancer?

I saw this book which mentions ArrayExpress, but I need some help to understand what kind of data it stores.

Are there similar/better resources? How can I use EMT in the context of cancer progression study?

Are there other process like EMT that are being used to study cancer progression?
Are there EMT like events in leukemia progression?

EMT describe cells' change in their state from being epithelial to the mesenchymal class. So, if a cancer cell line is gaining properties that allow it to move, it might reach the metastasis phase, were a cell can move to any point of the body and start a new tumor there.  As in the wikipedia page you linked:

  EMT has also been shown to occur in wound healing, in organ fibrosis and in the initiation of metastasis in cancer progression.

I am not aware of any kind of database for the EMT process, the most relevant resource might be the gene ontology, which has a term for it, which lists 185 (human) genes as related to this biological process.

ArrayExpress is a database of transcriptomics, it :

  stores data from high-throughput functional genomics experiments

You can evaluate the role of EMT in the data stored in the database.


Where can I get the population allele frequency vcf file? I want to use GATK to estimate cross-sample contamination for Whole Genome Sequencing data.

The specific tool is ContEst and it is run with:

java
 -jar GenomeAnalysisTK.jar \
 -T ContEst \
 -R reference.fasta \
 -I:eval tumor.bam \
 -I:genotype normal.bam \
 --popfile populationAlleleFrequencies.vcf \
 -L populationSites.interval_list
 [-L targets.interval_list] \
 -isr INTERSECTION \
 -o output.txt

The --popfile option requires a vcf file representing population allele frequencies. Where can I get it? Is it the one I can get from dbSNP?

I'm using the human_g1k_v37.fasta reference, a version of hg19, "for the phase1 analysis we mapped to GRCh37. Our fasta file which can be found here called human_g1k_v37.fasta.gz, it contains the autosomes, X, Y and MT but no haplotype sequence or EBV".

On the GATK forum they've recommended the population stratified VCF file for this purpose.


TPM or rlog(CPM) for comparing expression? I want to see the expression of a gene in a group of patient amongst the entire cohort using my RNA-Seq data. While I can do a differential expression analysis with limma or DESeq2, I want to see how much each sample from my cohort expresses the gene.

The plan is to plot a waterfall plot (as defined in this paper). However, when I used TPM and rlog(CPM) they gave me very different looking graphs (Not just the shape, which is expected, but the ranking of expression). My concerns are, for TPM, the batch effect with two or three lots in one cohort, and, for rlog(CPM), the results don't match with prediction but is adjusted for batch effect. Which one would be a right or better choice?

Upper one is with plot with. Lower one is plot with voom transformed counts (log(CPM)). I didn't include all cases since there are too many. This look more or less what it is like in the one including every case. Target is the condition group in which I am interested. Differential expression analysis definitely showed that it is upregulated.

                     logFC       CI.L       CI.R    AveExpr         t      P.Value    adj.P.Val         B
GeneA  1.6576257  0.9143490  2.4009024  5.8713219  4.391796 1.646627e-05 0.0002581676  2.484413

Working steps
TPM is basically from kallisto

log(CPM) calculation

#Import data from tximport
txi = tximport(files = read.hd5s,
                     type = "kallisto",
                     tx2gene = tx2gene,
                     countsFromAbundance = "lengthScaledTPM",
                     ignoreTxVersion = TRUE)
raw_data = DGEList(txi$counts)
# Filter low count genes
keep <- filterByExpr(raw_data)
raw_data <- raw_data[keep,,keep.lib.sizes=FALSE]
# Normalise counts
raw_data = calcNormFactors(raw_data)
# log(CPM) by voom
voom_data = voom(raw_data, design = NULL, normalize="none", plot = TRUE)
# Batch effect correction
voom_expr_psva = psva(as.matrix(voom_data$E),Batch)

Then I plot graph with TPM and log(CPM) respectively

For comparing the counts of different samples from DESeq2, Michael Love recommends using the variance-stabilized transform. It'd be great if you could provide some specific code examples in your question, but without that here's something that should work with the DESeq2 workflow as mentioned in the package documentation:

## [assuming "dds <- DESeq(dds)" has been run]
dds.counts <- assay(vst(dds, blind=FALSE));
## [This makes the scaled counts resemble more closely the actual read counts]
dds.quantile99 <- quantile(dds.counts[dds.counts > min(dds.counts)], 0.99);
dds.counts <- ((dds.counts - min(dds.counts)) / (dds.quantile99-min(dds.counts))) * (dds.quantile99);

The matrix dds.counts should now represent a model-normalised approximation of the log of counts, with column names as the sample names, and row names as the gene names. The default behaviour of the VST changes zero counts to an arbitrary positive number, which is why I've modified the output so that count values of zero will be changed to 0 in log space (i.e. a count of 1). To plot results for a single gene, subset based on the row name:

gene.VST <- sort(dds.counts["geneOfInterest",]);
plot(gene.VST);


What is a samtools mpileup reference skip? The samtools documentation for mpileup states:

  At this column, a dot stands for a match to the reference base on the forward strand, a comma for a match on the reverse strand, a '>' or '<' for a reference skip
  ...
  Similarly, a pattern `-[0-9]+[ACGTNacgtn]+' represents a deletion from the reference.

What is a reference skip?

That documentation has been expanded recently (see PR #1055) and now describes the characters seen in this column in more detail:

Forward   Reverse                    Meaning
---------------------------------------------------------------
 . dot    , comma   Base matches the reference base
 ACGTN     acgtn    Base is a mismatch to the reference base
   >         <      Reference skip (due to CIGAR "N")
   *         *      Deletion of the reference base (CIGAR "D")

(And so on — see the pull request or doc/samtools-mpileup.1 on the develop branch for the complete pileup format description. The documentation on that web page will be updated when this is released in a new samtools release.)

As noted in https://bioinformatics.stackexchange.com/a/421/134, SAM has always had separate D (deletion from the reference) and N (skipped region from the reference) operators. These are treated similarly computationally but emitted intentionally by aligners to distinguish introns and the like from other deletions.


Merging matched parts of two dataframes I have a big data frame like this

> dim(data1)
[1] 474303    128
> head(data1[1:2,1:6])
   Tumor_Sample_Barcode Chromosome Position
1:                   A1       chr1       15772155
2:                   A1       chr1       16054813
   End_Position Reference_Allele Tumor_Seq_Allele2
1:     15772155                G                 A
2:     16054813                G                 A
>

I have another dataframe like this

Chrom   Position    Ref Alt VAF
chr1    2716062 C   A   0.323
chr1    3095112 T   G   0.365
chr1    3342771 C   G   0.123
chr1    3526376 C   T   0.365
chr1    4121415 T   G   0.225
chr1    4215706 G   A   0.125
chr1    4304014 T   G   0.286
chr1    4304016 A   C   0.302
chr1    4388809 A   C   0.135

I want to create another column named i_TumorVAF_WU in my big data frame in which VAF column comes from small data frame. Please notice that VAF values should be matched with Position in big data because small data coming from big data

Can you help please?

    > unique(merge(data1,data,x.by=data1$Position,y.by=data$Position))

Returns a weird output with losing a big part of big data frame. I just want to add corresponding VAF values to the big data by keeping all the information in big data file

If you have issue with memory and dealing with large object, maybe data.table is the way to go (https://github.com/Rdatatable/data.table/wiki):

Let's take two dummy data.frame:

dfA <- data.frame(Chrom = rep("chr1",150),
                  Pos = 1:150,
                  ref = rep("A",150))
dfB <- data.frame(Chrom = rep("chr1",10),
                  Pos = 11:20,
                  Vaf = rep("VAF", 10))

The big data.frame dfA:

> head(dfA)
  Chrom Pos ref
1  chr1   1   A
2  chr1   2   A
3  chr1   3   A
4  chr1   4   A
5  chr1   5   A
6  chr1   6   A

and the small dataset:

> dfB
   Chrom Pos Vaf
1   chr1  11 VAF
2   chr1  12 VAF
3   chr1  13 VAF
4   chr1  14 VAF
5   chr1  15 VAF
6   chr1  16 VAF
7   chr1  17 VAF
8   chr1  18 VAF
9   chr1  19 VAF
10  chr1  20 VAF

Now, using data.table that is more efficient and less memory consuming for the manipulation of dataframe, you can do:

library(data.table)
setDT(dfA)
setDT(dfB)
dfA[dfB, on = c("Chrom","Pos"), vaf := i.Vaf]

And now, your data.frame dfA looks like:

> dfA[5:25]
    Chrom Pos ref  vaf
 1:  chr1   5   A <NA>
 2:  chr1   6   A <NA>
 3:  chr1   7   A <NA>
 4:  chr1   8   A <NA>
 5:  chr1   9   A <NA>
 6:  chr1  10   A <NA>
 7:  chr1  11   A  VAF
 8:  chr1  12   A  VAF
 9:  chr1  13   A  VAF
10:  chr1  14   A  VAF
11:  chr1  15   A  VAF
12:  chr1  16   A  VAF
13:  chr1  17   A  VAF
14:  chr1  18   A  VAF
15:  chr1  19   A  VAF
16:  chr1  20   A  VAF
17:  chr1  21   A <NA>
18:  chr1  22   A <NA>
19:  chr1  23   A <NA>
20:  chr1  24   A <NA>
21:  chr1  25   A <NA>
    Chrom Pos ref  vaf

As you see, we incorporate the column Vaf from the small dataframe into the big dataframe.

With your example, I will first rename the colnames of your small dataframe in order to have the perfect match with the big dataframe:

colnames(data) <- gsub("Chrom", "Chromosome", colnames(data))

then,

setDT(data1)
setDT(data)
data1[data, on = c("Chromosome","Position"), i_TumorVAF_WU := i.VAF]

Does it look what you are trying to obtain ?