Visually summarizing phylogenetic networks After having worked with Bayesian phylogenetic tree inference for some time, I am now trying out a method to infer network phylogenies. To get an idea of what it infers, I'd like to find a way of summarizing the results.

The output of the inference is a NEXUS file containing a stream of networks encoded in eXtended Newick as trees with duplicate node names, like this:

(((((((Kiraman:0.028228,Kui:0.028228)S1:0.0797301)#H4:0.35229)#H2:23.4168)#H1:15.8027,(Kamang:29.3316,#H1:5.45455[&&NHX:gamma=0.36825381342449554])S11:10.3482)S12:1.27331,((Klon:1.16727,((#H4:0.524833[&&NHX:gamma=0.33697036713230577],Kafoa:0.632791)S6:0.0615575)#H3:0.472921[&&NHX:gamma=0.6642547809404831])S2:35.6082,(#H2:22.2107[&&NHX:gamma=0.7544913491697585],((A-Petleng:0.421975,(A-Takalelang:0.168615,(A-Ulaga:0.145472,A-Atimelang:0.145472)S4:0.0231422)S5:0.253361)S7:0.902144,#H3:0.629771)S3:21.3468)S8:14.1045)S9:4.17761)S10:54.5301):0;

￼

There is no formal restriction on the shape of the network, and hybridization edges have a weight 0<γ<1. How do visually or formally summarize a stack of several thousands of these networks, and in particular the backbone tree structure plus a glimpse of the hybridization events they display?

I have played with some tree summarizing approaches and tried to first infer a summary tree and then glue some cross-edges in there, but the way I came up with did not look reasonable.

I am now considering to calculate distances between leaves and using network visualization tools like SplitsTree or NeighborNet to get a visual representation of that, but how shoud I represent the uncertainty of the distances in even a single network (not to mention thousands of them) in the input for these algorithms?

Are there any good approaches to summarize phylogenetic networks of the shape given above? Or to consistently simplify them to a shape that can be easily summarized?

Firstly you treat this as two separate independent calculations.

The traditional approach is simply to produce a consensus phylogeny using your Bayesian method and assess the equivocal bootstraps/Bayesian probability.
Feed your data set (not the trees) into NeighborNet within SplitsTrees. There are statistics that classify the level of departure from a bifurcating phylogeny.

Examine both outputs for orthogonal support for incongruence.

Note Phylogenetic network statistics are a relatively soft area of phylogeny (they are not used as hard evidence).

Networks from Bayesian trees

The final approach (which is what you are asking about), is to feed your Bayesian trees into Dendroscope, this is a GUI that represents the uncertainty between different bifurcating phylogenies as networks and you can use it to understand the incongruence between different trees within your Bayesian MCMC, because the phylogenetic networks directly corresponds to the level of incongruence.

Dendroscope is a nice Java API and is very easy to use.


Parsing the json file from gene ontology result into dataframe How do I parse the json file from ontology output.
My file json file
jdata <- read_json("analysis.json", simplifyVector = TRUE)

I want to convert it into a data.frame format. I can able to take the GO term and respective p value etc. But I would like to know the genes that are getting part of the enriched list which is under the upload_1 for example under DNA Replication there are 33 terms but 8 are from my input. So my objective is to get the 8 gene so on and so forth.
So I would like to see a data.frame where I want to have GO term in the first column gene mapped into the term would be second i.e upload_[#] third one fold enrichment and last one p value.
This is how i would like to see this one is the text output in the column upload_1 (450) instead of that number i would like to see then gene names which are present in the json file.
Some progress here is my code
l <- fromJSON("analysis.json")

m <- lapply(
  l[[1]]$result,
  function(x) c(x$input_list$mapped_id_list['mapped_id'], x$term['id'],x$input_list['pValue'],x$term['label'])
)
m <- do.call(rbind, m)

output:
mapped_id     id           pValue       label
[1,] Character,8   "GO:0006270" 0.007540801  "DNA replication initiation"
[2,] Character,388 "GO:0009987" 3.232843e-10 "cellular process"

Right now I can parse out first two row.
Things I would like to append to the code is

I would like to do for all the terms
The first column where i am trying to get all the genes now it is coming as character while number suggests gene present. That I wont need probably since i need only gene names which is being shown as character now.

dput(head(b))
    structure(list(`GO biological process complete` = c("DNA replication initiation (GO:0006270)",
    "regulation of blood coagulation (GO:0030193)", "regulation of hemostasis (GO:1900046)",
    "regulation of coagulation (GO:0050818)", "platelet degranulation (GO:0002576)",
    "cell cycle G1/S phase transition (GO:0044843)"), `upload_1 (450)` = c(8,
    13, 13, 13, 16, 15), `upload_1 (fold Enrichment)` = c(11.23,
    7.09, 7, 6.62, 5.75, 5.35), `upload_1 (P-value)` = c(0.00754,
    0.000627, 0.000716, 0.00135, 0.000347, 0.00228)), row.names = c(NA,
    -6L), class = c("tbl_df", "tbl", "data.frame"))

Any suggestion or help would be really appreciated.

I don't know how to do that in R but I think this is what you want:
I used http://json2table.com for view the data structure of the JSON file, then I used json_normalize to normalize the JSON data into a flat table using the key [result] to a list of DataFrames for each index in the JSON file. Then, I concatenated the list of DataFrames in a single one with df.concat()
import json
import pandas as pd

with open('analysis.json') as f:
    d = json.load(f)

dfs = []
for x, y in enumerate(d):
    data = pd.json_normalize(d[x]['result'])
    dfs.append(data)

df = pd.concat(dfs, ignore_index=True)

df.to_csv('analysis.csv', index=False)


Common QC parameters and threshold in human WGS pipeline I am trying to summarize a list of QC-parameters and its threshold that should be applied for human WGS in general to avoid bad quality of data.

Which qc parameters for WGS are you using and why?

Thanks!

I would suggest you look at the website for the excellent MultiQC tool (https://multiqc.info/).  We use it in a clinical setting to collate all the QC statistics from our pipelines into a single report (It displays results from multiple QC programs such as FastQC, Picard, bcl2fastq, VerifyBAMid etc).  The homepage has a selection of example reports - clicking on "Whole Genome Seq" shows a report with the metrics most labs would expect from their WGS pipeline.

If you dig around in their GitHub repo (https://multiqc.info/) you will get a good idea of the default thresholds set for flagging results - there should be a config file somewhere specifying all these values.


Identity plus Similarity Interpretation I am helping a colleague to interpret protein multiple sequence alignment results. My colleague used a third-party proprietary solution (I know...) to calculate the alignment, and one of the resulting tables is a pairwise distance matrix with identity % in the top triangle and identity plus similarity in the bottom triangle.

I am struggling to interpret identity plus similarity, despite the fact I found the metric in several papers.

I guess you already know what identity means. But for completeness, it is the % of pairwise alignment position which have an identical amino acid residue. For the identity+similarity it is important to know that amino acids can be grouped according to their physicochemical properties e.g. charge, size etc. While there are standard groupings for amino acids it may be that the proprietary algorithm uses different groupings. Either way, % identity+similarity would imply the number of pairwise alignment positions which contain either the same amino acid or amino acids from the same group (however they define groups).


vcf file to gene annotations and or genes Is there some tool that will convert a vcf/bcf file to some gene annotation, based on the fact that I have the gene annotation of the already exsisting reference genome. Or some tool that will give me the genes based on the most probable chance.

Creating a FASTA file by applying variants found in a VCF file is called "creating a consensus FASTA file". You can use bcftools consensus. Unfortunately, the tools only works genome-wide (at the FASTA level), so you might need to use a different tool to extract a gene sequence from the FASTA file using the gene's co-ordinates.

If the VCF file is annotated with HGVS notations for the variants, there might be a way to use those notations to apply the edit to the reference gene/transcript sequence directly. I know the hgvs python module works well with the HGVS notation, but I don't see a way to "apply" the variant to the reference sequence.


Modify how a script prints the output? I'm trying to print out all measured distances between a ligand and 2 residues (using PyMol), below is the script that I'm using. I only want to measure the distances between all the ligand's O & ND2 of N312, as well as all the ligand's N & OD2 of D113.
from pymol import stored
stored.a, stored.b, stored.c, stored.d = [],[],[],[]
cmd.iterate_state(1, "r55_01///UNL`1/N*", "stored.a.append(ID)")
cmd.iterate_state(1, "r55_01///UNL`1/O*", "stored.b.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASP`113/OD2", "stored.c.append(ID)")
cmd.iterate_state(1, "/r55_01//R/ASN`312/ND2", "stored.d.append(ID)")

f=open("dst4.txt", 'w')

for a in stored.a:
  for b in stored.b:
    for c in stored.c:
      for d in stored.d:
        f.write( "distance from %s to %s is %s\n" % (a, c,
cmd.distance( "id %s" % a, "id %s" % c,cutoff=5,mode=2)))
        f.write( "distance from %s to %s is %s\n" % (b, d,
cmd.distance( "id %s" % b, "id %s" % d,cutoff=5,mode=2)))

f.close()

The output file currently looks like this
distance from 9 to 1334 is 4.589359760284424
distance from 8 to 4204 is 3.1509475708007812
distance from 9 to 1334 is 4.589359760284424
distance from 15 to 4204 is 0.0
distance from 9 to 1334 is 4.589359760284424
distance from 17 to 4204 is 0.0
distance from 9 to 1334 is 4.589359760284424
distance from 20 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 8 to 4204 is 3.1509475708007812
distance from 23 to 1334 is 0.0
distance from 15 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 17 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 20 to 4204 is 0.0

I want my output to look more like this. Thanks in advance!
Lig_N D113 4.59
Lig_O N312 3.15

Let's simplify...
cmd.distance(xx,xxx) accepts two selector strings and is fine with and state 1. And best of all if you select multiple atoms it gives you the smallest distance!
So all you need to do is:
# measure
d_O = cmd.distance('UNL and element O', 'resi 312 and name ND2')
d_N = cmd.distance('UNL and element N', 'resi 113 and name OD2')
# write
f=open("dst4.txt", 'w')
f.write(f'Lig_N D113 {d_N:.2f}\n')
f.write(f'Lig_O N312 {d_O:.2f}\n')
f.close()

If you are using python 2 (if so do note pymol scripting is better in 3), then the second bit is
f=open("dst4.txt", 'w')
f.write('Lig_N D113 %.2f\n' % (d_N,))
f.write(f'Lig_O N312 %.2f\n' % (d_O,))
f.close()

Unsolicited coding tips
A few unsolicited tips... Your original script is rather over-complicated and obfuscated and uncommented, so for better legibility why not do:

cmd.iterate_state is always harder to follow than the for atom in get_model('xxxx').atoms synthax
unless you are stuck with python 2, use f-strings not the %s sprintf style formatting
the /x/x/x/x selector synthax is never as readable as the longer
resi/resn/name selection algebra.

Code golf is always tempting, but it can result in future-you getting stuck as happened here...


What are Approximate Read Counts (Library Sizes) and Lengths (Insert Sizes) for Next-Generation DNA Sequencers? I am performing a simulation study and am curious about the parameters of my simulated metagenome.

What are the library and insert sizes of some of the most "used" sequencers. Mainly, I am interested in the different Illumina platforms (MiSeq, HiSeq, HiSeqX), 454 Pyrosequencing, and Ion Torrent.

I do not need the long read technologies like Pac-Bio's SMRT and Oxford Nanopore since they are not typically used with metagenomics, yet.

Both of those depend on the application of the library.  For example whole genome assembly of a mammal would require a larger library size than targeted re-sequencing of a virus.  Also insert size varies if you are doing genome assembly a larger insert size is valuable to allow for long range linking and scaffolding.

I think that if you edit your question of be more specific people will be able to provided a more concrete answer other wise it is "it depends"


How to calculate statistical significance of sequence motifs I would like to show that certain types of RBP motifs are enriched in RNA editing islands (i.e. clusters of RNA editing). However, I am unsure about how to think about sequence motifs with respect to their occurrence in other genomic features.

I understand how to find the probability of a motif at a location. i.e.

P(sequence is at position i) = P(A)^[A] * P(C)^[C] * P(G)^[G] * P(T)^[T]

where (.) is the base pair and [.] is the number of bps.

Instead what I would like to find is:

P(sequence S is contained in a feature type T) = ???

where feature type T is an gene, intron, editing island, etc. I think I should incorporate length since I will mainly be comparing genes or introns vs. editing islands. Also, I am not sure what to do about the editing islands being located within genes. How can I keep from counting the same motif twice?

Any ideas would be greatly appreciated. Thank you for your time.

As far as I remember the exact probability computation is an open problem. The reason is that potential motifs can overlap, which makes probability computations for an arbitrary string non-trivial, and depends on the motif.

For example if he have a binary string of four digits, the probability of "01" will be 11/16, while the probability of "11" will be 8/16=1/2.

The simplest approximation is to assume that probability of the motif on every position of the sequence is equal and independent. In such case (this comes from this course):

(Probability of text of length N having at least t occurrences of a k-mer Pattern, A is the number of letters in the alphabet, n = N – t * k)

This approximation is very rough, there are better here.


How to add data to AnnotationHub I recently started using the Bioconductor AnnotationHub package for accessing lots of common datasets, but I've noticed that there are a lot of newer databases that have not made it into AnnotationHub's collection.

Has anyone had any luck getting databases added to AnnotationHub or adding them manually yourself (I guess via the package process?)?  Is it painful or relatively straightforward?

I love the annotationhub interface, but I'm wondering if it is way too inefficient to add new resources rather than just manually downloading and curating myself.  Thanks!

As the description of the package says:

  The AnnotationHub web resource provides a central location where
  genomic files (e.g., VCF, bed, wig) and other resources from standard
  locations (e.g., UCSC, Ensembl) can be discovered.

If the new resources aren't in the standard locations it won't be listed by AnnotationHub. The easy/common way to have the data accessible through AnnotationHub is make it available through these standard locations (before publications if you want to give early access).

Another method would be to follow this steps of the vignette to add new resources, but it might involve creating a new package.


RapMap: reference transcriptome for simulated reads I am very new to bioinformatics and trying to repeat the benchmark in the RapMap paper with an experimental tool working in a similar but different fashion.

In the paper (taken from their github) about RapMap, an alignment/mapping tool for rna-seq reads to a reference genome, the dataset used for their benchmarks is described like this:

  To test the accuracy of different mapping and alignment tools in a
  scenario where we know the true origin of each read, we generated data
  using the Flux Simulator (Griebel et al., 2012). This synthetic
  dataset was generated for the human transcriptome from an annotation
  taken from the ENSEMBL (Cunningham et al., 2015) database consisting
  of 86090 transcripts corresponding to protein-coding genes. The
  dataset consists of  48 million 76 bp, paired-end reads. [...]
  When benchmarking these methods, reads were aligned directly to the transcriptome, rather than to the genome.

I understand that they generated reads using flux simulator, but I don't understand which transcriptome they mapped them against. I can't find any other information in the paper about which annotation they used either.

I am assuming one would need to generate an entire transcriptome in .fa format from the GTF and genome? Or did they generate a transcriptome from only the simulated reads and then mapped the same reads against it?

GENCODE offers GTFs and FASTAs for gene-encoding transcripts here, might that be the kind of thing they used? Would the reads simulated from these GTFs fit to these transcript FASTA files?

Thanks for your interest in RapMap.

At that time we were using flux simulator for simulating read sequence data. We used the genome and gtf file together as an input to flux.

I dug into the scripts and got hold of the gtf link ftp://ftp.ensembl.org/pub/release-80/gtf/homo_sapiens/Homo_sapiens.GRCh38.80.gtf.gz,
and the genome file link ftp://ftp.ensembl.org/pub/release-80/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.

We used rsem-prepare-reference to create the reference transcriptome file later in order to run RapMap and other transcriptome based quantification tools. I hope this helps.

PS: I am a coauthor of the RapMap paper.


Details of ENCODE processing pipelines I am trying to find details for how the ENCODE ChIP seq samples were processed.

Linking to a ChIPseq samples on www.encodeproject.org (e.g. https://www.encodeproject.org/experiments/ENCSR051DXE/) offers a PDF "overview" of the pipeline which says basically nothing. The page will also offer you a map of the pipeline, which basically just has headings:

Clicking on a step gives you a popup with a link to the "Chip-seq read mapping" pipeline. This just contains the same overview PDF and a link to "ChIP-seq Data Standards and Processing Pipeline". This page contains again, a vague outline of the pipeline, a link back to the page we just came from and a list of quality control thresholds.

I can't seem to find anywhere the proper details. For example, if BWA was used, was it MEM or aln? What parameters were used for the mapping? One of the steps is read filtering: what where the filters? How were they applied? There is a link to a github page for running the pipeline on DNAnexus, but I can't make head nor tail of it.

The codes are usually available at DCC's GitHub.
You can always email them if you have more questions.


How to programmatically classify a protein according to its genbank feature Say I found an interesting protein in a genbank file, e.g.:
 gene            273792..275552
                 /locus_tag="ECs_0242"
 CDS             273792..275552
                 /locus_tag="ECs_0242"
                 /codon_start=1
                 /transl_table=11
                 /product="type IV secretion protein Rhs"
                 /protein_id="NP_308269.1"
                 /translation="[...]"

(from https://www.ncbi.nlm.nih.gov/nuccore/NC_002695.2)
I would like to programmatically get a classification of this protein.
The (not new) motivation is that my pipeline finds many proteins, and I want to identify enriched or depleted protein classes (e.g., enrichment for DNA-binding proteins).
By googling for "type IV secretion protein Rhs", I found a UniProt entry which has a link to a QuickGO page, which includes three categories: "transposase activity", "transposition, DNA-mediated" and "DNA binding". From this page I can also get an ancestor chart for each of these categories, e.g.:

Optimally, the code would work quickly (preferably locally), and for a description from a CDS feature (e.g.: "type IV secretion protein Rhs") return the GO categories it belongs to. It would be even better to also have the ancestry of each category in a computer-readable format.
What is the best way to do that?
Is there a specific UniProt subset that fits my need (finding GO categories given a protein description)?
Is there a specific GO subset that fits my need (extracting ancestor GO categories as explained here)?
(I figured that GO categories fit my need of classifying proteins, but if a better option exists, please let me know)

Using the UniProt website API (https://www.uniprot.org/help/api) this could be done with a query like
https://www.uniprot.org/uniprot/?query=name%3A%22type%20IV%20secretion%20protein%20Rhs%22&columns=id%2Centry%20name%2Creviewed%2Cprotein%20names%2Cgenes%2Cgo(biological%20process)%2Cgo(molecular%20function)%2Cgo(cellular%20component)&sort=score
or, in tab-separated format:
https://www.uniprot.org/uniprot/?query=name:%22type%20IV%20secretion%20protein%20Rhs%22&format=tab&limit=10&columns=id,entry%20name,reviewed,protein%20names,genes,go(biological%20process),go(molecular%20function),go(cellular%20component)&sort=score
(Remove the "&limit=10" part from the URL to obtain the full set instead of the first 10 entries.)
You can add colums as described here:
https://www.uniprot.org/help/customize,
https://www.youtube.com/watch?v=ado1r8IDm3U (video, at about 1min25)
or
https://insideuniprot.blogspot.com/2015/03
Note that while UniProt tries to improve protein name annotation, the unreviewed section of UniProtKB (i.e. UniProtKB/TrEMBL) is often far from perfect and its accuracy relies on the information provided by the submitter of the nucleotide entry, unless an automatic annotation rule is available to assign a more suitable name (see https://www.uniprot.org/help/protein_names).
Please don't hesitate to contact the UniProt helpdesk if you have any additional questions.


Metagenomic shotgun data with internal control For 24 human stool samples I have metagenomic shotgun reads from an Illumina platform. The reads were filtered and mapped against a bacterial species library and specific maps to species were kept and counted. Importantly, before DNA extraction from any sample it was spiked with a solution of a fixed DNA concentration of Thermus thermophilus, a species not occurring in the human microbiota. The reads from T. thermophilus allow me to convert relative abundance estimates back to abundance in the stool sample, to absolute abundance.

I want to test whether T. thermophilus is indeed a good internal control that removes confounding by converting to absolute abundance. To test this, I compare between-species correlations on the relative abundance scale with their correlations on absolute abundance scale. I would expect less correlation in the latter. I first take the raw read counts of a random set of species, I divide counts by a sample-specific size factor estimate (total number of reads r in sample i / geometric mean of r's) to adjust for library size, then I do arcsinh transform (to stabilise variance) and then I look at between species correlation. Next I take again the raw read counts, correct by the size factor, then correct for an estimate of the coverage ( estimated proportion of T. thermophilus in sample i ), and I do the acsinh transform. Again I look at between-species correlations.

I compare species-pair correlation before and after T. thermophilus correction in a scatterplot (see attached) to visualise the effect of the correction. I expected that the point cloud indicates a slope <1 reflecting that previously strong correlations have weakened on the absolute scale. But that is not what my result shows. Indeed, some weakly positive correlations have become stronger. How does this make sense? What other ways are there to test validity of an internal control?

The problem with correlations in relative scale is know on the literature, it has been suggested to change to other metrics. So I wouldn't use it for internal control.

As you can see on figure 3 of the linked article the correlation on the absolute scale tend to be higher. So I would say this is within expectations.

As other way for internal controls are blanks (Just water to see if there are any contaminants), and technical replications, to make sure that the measure is correct not some problem with picking some amplifications instead of others...

Also you don't mention it but the tool and library used to impute the microorganism is also important. You can check with other tools to see if the estimation is similar.


How to get a MSA fasta from BAM/SAM? I am working with paired end NGS miseq data from a viral genome with multiple timepoints.  I have filtered and trimmed this data for quality and adapter sequences.  I have then merged the filtered paired end reads and mapped them to a reference genome using bowtie2.  I then used samtools to filter for mapping quality and visualized the resulting filtered bam file on IGV.
My goal is to now scan along the genome using various window sizes (e.g. 300 nt) across the reads that have mapped to the reference. I will then chop up the genome and build phylogenetic trees based on each window.
How can I extract a multiple sequence alignment of these 300 nt windows from a sam or bam file using the coordinates of a reference viral genome?  I have tried using samtools view but this does not actually provide the subsection or alignment within the window but rather all reads that overlap the window in some way.
The virus is HIV1 and the immediate goal is to get an alignment of all reads that completely span a user specified coordinate range eg. 1-300 based on some reference sequence. Ideally reads that fit this criteria would be extracted and aligned. After ragged edges are removed, the resulting alignments subjected to phylogenetic analyses.

If need to get msa between every read and ref in the specific region ,
you can try this https://github.com/orangeSi/bam2msa.git

bam2sam ref.fa input.bam --regions chr1:200-300


Bedtools get fasta and ORF from a blastX run Hello I made a blastX research  (A query genome in nucleotide format translated into protein in 6 reading frames against a protein dabatase)

And here is a head of the result :

IDBA_scaffold_7517  493 YP_009316100    0.360   61  27  1   150 4   126 186 6.775E-03   37
IDBA_scaffold_29149 519 AIG51500    0.428   35  18  1   2792    2896    392 424 2.712E-03   40

I have 2 hits:

The first one is the query IDBA_scaffold_7517 that have a hit with a protein sequence YP_009316100 in the (-) strand because the start begins by 150 and ends with 4.

The seconde one is the query IDBA_scaffold_29149 that have a hit with a protein sequence AIG51500 in the (+) strand because the start begins 2794 and ends by 2893.

so with that I actually created a bed file in order to extract the fasta sequences from my genome ( I assume that since there is a hit with a protein in the Genomes the coordinates of this hit should be in a good ORF (without stop codon).

Here is the bed_file I made :

IDBA_scaffold_7517  4   150 Species 0   -
IDBA_scaffold_29149 2794    2893    Species 0   +

Then in order to extract the sequences from the coordinates in the Genome I used:

bedtools getfasta -fi Genome.fa -bed bed_file -s -fo test

with -s to force strandedness. If the feature occupies the antisense, strand, the sequence will be reverse complemented.

But when I do that the sequences extracted from the query are not in the good ORF.
I guess it is because I should take into account the fact that the coordinates are from nuc translated to portein but how could I do that? Should I change the coordinates manually from the bed file according to the fact that it is a negative of positif strand ?

Thank you for your help.

Blast uses a 1-based coordinate system. Bed is 0-based half open. Your code needs to convert appropriately. See for example: https://www.biostars.org/p/317112/ .


Where can I find a GWAS data set about genotypes and phenotypes for Zebrafish and rice? My goal is to work on GWAS method development as a research project. However, human genome data are usually confidential because of the identification problem, so it's very hard to get them. (application takes a long time and I may not get it at the end.)

Therefore, maybe I can try to play with the data of other species first. So I wonder where I can find some GWAS data for Zebrafish and rice. (The reason I am asking for these two is that I find researchers tend to perform genome analysis to mice, Arabidopsis Thaliana, zebrafish, and rice, and I already found where to download the GWAS data for the first two.)

I can find several data sets (like ZFIN) that have the annotated genome information, but that's only one genome. Where can I download the raw data (like genotypes of 500 zebrafishes and corresponding phenotypes) to start my own analysis?

I think it will be more productive for you if you choose one of the public human datasets and practise your GWAS methods on that. A couple of examples examples:

1000 genomes -- use "population" as your trait of association. You'll find that you get association across the whole genome, so the analysis method would change from "what bits are associated" to "what are the bits that are most consistently associated"
NCI Genomic Data Commons -- open access (and controlled access) data associated with cancer


What is the best method to estimate a phylogenetic tree from a large dataset of >1000 loci and >100 species I have a large phylogenomic alignment of >1000 loci (each locus is ~1000bp), and >100 species. I have relatively little missing data (<10%).
I want to estimate a maximum-likelihood phylogenetic tree from this data, with measures of statistical support on each node.
There are many phylogenetics programs that claim to be able to analyse datasets like this (e.g. RAxML, ExaML, IQtree, FastTree, PhyML etc). Given that I have access to a big server (512GB RAM, 56 cores), what are the pros and cons of each program.
Which programme likely to give me the most accurate estimate of the ML tree for a dataset of this size within a reasonable timeframe?

This paper claims that FastTree is almost as accurate as RAxML, while being much faster. You just have to be careful, however, that the support values output by FastTree are not bootstrap values, they are based on the Shimodaira-Hasegawa test. (Also, see this comment for the case you have very short branch lengths).   [update: However, according to the recent comparison paper mentioned below FastTree performed quite poorly in comparison to RAxML or IQ-tree.]

From what I understand, you should use ExaML only if your data is too large to be handled by RAxML in a single node. ExaML should perform like RAxML but with some parallelization overhead. For all effects I treat them as the same. I don't know of relevant advantages of phyML over RAxML (for me, it's easier to use but I am very used to phyML).

I am not familiar with IQ-tree, but its authors claim that even given the same time as RAxML or phyML, IQ-tree already finds better likelihoods more often than not (although by default it takes a bit longer to converge). A recent comparison between all these programs favoured IQ-TREE for both single-gene and concatenation analyses (with RAxML very close). It may also estimate branch support through a SH-like test only, but I'm not sure.   [update: IQ-tree offers 3 measures of support, standard bootstrap, aLRT, and ultrafast bootstrap. See OP's comment below for details.]

However, since you have few missing data, you might also want to try a single-locus tree inference followed by gene tree clustering (using treescape or treeCL) to see how spread your data is, or to see the effect of removal of outliers, or to use ideas similar to statistical binning.


Where is an up to date miRNA database and what happened to miRBase? miRBase 21 was published June 26, 2014 and was still in its growth phase. Why is it not being updated anymore or the project declared officially dead? ENSEMBL also uses miRBase as a starting point (http://www.ensembl.org/info/genome/genebuild/ncrna.html).

What do people use now as a good reference for miRNAs especially for less well annotated/sequenced species?

If you are still interested, last year miRBase generated new updates. Currently, according to ftp site the last release is 22.1. So, it is not a dead project and for more specific information you should reference the miRBase blog.
In order to get the best set of miRNA annotations, you have to define your target species or even if exists some related annotated species and then, retrieve the current annotation along different miRNA databases, but most importantly you have to be really clear about the definition of miRNA families, because it is clearly different in miRBase and also in RFAM.


Looking for a way to download all pdb files that detail a protein-drug complex I am looking for a way to download all pdb files that detail a protein-drug complex.

While I could download the entirety of the pdb and search manually (or maybe write a script to discard all proteins that do not have an HETATM entry that details a non-water, non-salt molecule), that sounds effort intensive and error prone.

It seems reasonable that this subset already exist, but I was not able to find it.

Does such a list/database exists?

Proving a negative is hard, but to my knowledge I do not know a database that has all the structure of drug bound protein. I can give some close matches.

First a few terms to be on the same page. A ligand is something that binds, which can be physiological, such as a cofactor or substrate, or not (especially man-made), such as a drug, i.e. a compound that has a physiological effect, or a hit, a small fragment from a fragment screen that will not bind with nanomolar affinity. Generally, ligands are small molecules, but may include oligopeptides and oligonucleotides. In the case of crystal structures, you also have crystallisation additives, which are not ligands, but are hetatoms in PDBs. Also, the definition of drug may vary. drugbank.ca has the following categories of drugs "approved", "experimental", "withdrawn" for pharmaceuticals, "illicit" for narcotics and "investigational" for probes and counts macromolecules as drugs.

RCSB PDB: drug mapping table is a simple searchable table with ligands that are clinical drugs only —experimental compounds (=probes) will be absent.
bindingdb.org contains many table, including a '3D' one. This table contains all the PDBs with ligands and their affinities. This contains crystallisation additives and natural ligands. Peter Curran's repo for the HotSpots algorithm (first written by Chris Radoux) has a blacklist, but you are bound to get stuff missed.
bindingmoad.org is a database with all the protein-ligand structures from PDB up until 2018, but it includes natural ligands, like cofactors and substrates as above.
https://klifs.vu-compmedchem.nl a well curated database of kinases, which has a lot of members and many ligands, but is not at all diverse. As a result it is a clean dataset often used as a testbed for novel methodologies.
I mentioned in the comments the SGC Target Enabling Packages as a small subset of proteins with many hits each. But I realise that there is not a handy PDB list and the targets are not drugs, but probes.


Differential Expression With 2 Treatments I have a expression data from a small cohort of samples taken at baseline and after 2 independent treatments. I can do differential expression contrasting T1 and T2 or I can contrast T1 vs baseline and T2 vs baseline and look at the differences. What is the difference between these two analyses? Are they both valid? If so what inferences can be drawn from each?

If you're interested in looking at the differences between two treatments then you'll end up wanting to do both a direct contrast as well as the individual comparisons to baseline.

The direct contrast will give you the genes actually differentially expressed between the two conditions. In practice, you may want to filter this a bit so you only have genes differentially expressed vs. baseline in at least one condition (e.g., to get rid of genes only slightly higher due to T1 and slightly lower in T2, but not different enough in either case to be DE). Use a pretty lax p-value threshold (e.g., 0.1 or 0.2) for this.


Classifying samples based on marker gene expression I have a few sets of marker genes that I can classify RNA-seq samples using semi-supervised clustering.  I would like to automate the process, however, I am struggling to find the ideal algorithm that could generate some kind of score for marker gene set from a given sample.

I presume that this is a standard analyses in many groups but I am not sure which method(s) are yielding good results in practice.

I would consider using gene expression signatures to classify samples (especially cancer subtypes but the same principles apply to other problems of this type) one of the classic problems of bioinformatics. Quite a lot of work has been done on methods to derive gene sets that provide good classification performance. This is slightly different from your problem since you already have a gene signature but it may still prove useful.

These methods will typically fit a model that selects a (small) number of genes from genome-wide expression data that distinguish between the cell types/conditions in question, i.e. they derive a gene signature. The resulting model then allows for the classification of new samples. I've had success using GeneRave for this purpose (but note that this was designed for microarray data, I haven't used it with RNA-seq data and don't know how well it holds up there). A more recent paper relating to this issue can be found here.

So how does that help you? One option would be to fit one of these classifiers to gene expression data for the genes you already know to obtain a model that can then be applied to new samples automatically.


What are the differences between BNGL (BioNetGen Language) and SBML (System Biology Markup Langague) formats? Per request from meta comment.

I am self-learning about whole cell modeling, specifically An introduction to whole-cell modeling and Fundamentals of Systems Biology: From Synthetic Circuits to Whole-cell Models. While my background is decades of programming, for cellular biology it is only skimming through a few books.

In the process of learning I soon discovered software/apps that use SBML (System Biology Markup language) but in referring back to "An introduction to whole-cell modeling" see that they use BioNetGen (GitHub) which has its own language BioNetGen language (BNGL)

What are the differences between BNGL (BioNetGen Language) and SBML (System Biology Markup Langague) formats?

Both formats have formal grammars for describing similar concepts, however they are different and to my knowledge one is not a proper subset of the other. While one could do a comparison of a syntactic nature, a simple Google search of just both acronyms reveals BioNetGen 2.2: Advances in Rule-Based Modeling

In here can be found:

  2.4 SBML-to-BNGL translation
  SBML is a widely-used model exchange format in systems biology (Hucka et al., 2003). Models encoded in SBML
  are ﬂat, i.e., their species do not have internal structure, which
  limits their utility for rule-based modeling. BioNetGen 2.2 includes
  an SBML-to-BNGL translator, called Atomizer (also available as a web
  tool at ratomizer.appspot.com), that can extract implicit molecular
  structure from ﬂat species (Tapia and Faeder, 2013). A full report on
  Atomizer and its application to the BioModels database (Li et al.,
  2010) is currently in preparation. However, Tapia and Faeder (2013)
  reported that an early version of the tool could recover implicit
  structure for about 60% of species in models within the database that
  contain ≥20 species. Thus, Atomizer makes available a large set of
  pre-existing models in a rule-based format, facilitating their
  visualization (Wenskovitch et al., 2014) and extension (Chylek et al.,
  2015).

While not directly related to the question this also helped in understanding more about BioNetGen

Wikipedia: Multi-state modeling of biomolecules

  Biological signaling systems often rely on complexes of biological
  macromolecules that can undergo several functionally significant
  modifications that are mutually compatible. Thus, they can exist in a
  very large number of functionally different states. Modeling such
  multi-state systems poses two problems: The problem of how to describe
  and specify a multi-state system (the "specification problem") and the
  problem of how to use a computer to simulate the progress of the
  system over time (the "computation problem"). To address the
  specification problem, modelers have in recent years moved away from
  explicit specification of all possible states, and towards rule-based
  formalisms that allow for implicit model specification, including the
  κ-calculus, BioNetGen, the Allosteric Network Compiler and others. To
  tackle the computation problem, they have turned to particle-based
  methods that have in many cases proved more computationally efficient
  than population-based methods based on ordinary differential
  equations, partial differential equations, or the Gillespie stochastic
  simulation algorithm. Given current computing technology,
  particle-based methods are sometimes the only possible option.
  Particle-based simulators further fall into two categories:
  Non-spatial simulators such as StochSim, DYNSTOC, RuleMonkey, and
  NFSim and spatial simulators, including Meredys, SRSim and MCell.
  Modelers can thus choose from a variety of tools; the best choice
  depending on the particular problem. Development of faster and more
  powerful methods is ongoing, promising the ability to simulate ever
  more complex signaling processes in the future.


How do I find split reads? How can I detect a split read in a BAM file?
Is there any sign in the CIGAR string that describes split read?

The answer to this question depends on the aligner used to generate SAM/BAM file. For BWA method that is part of the pipeline mentioned in comments, split reads will be flagged as supplementary alignment (2048). You can modify this behavior by using option -M as described in the BWA-MEM manual page.

So, in case of split alignment, you will end up with multiple aligned segments that have identical query name and they are flagged as 2048. These segments may only have part of the original read sequence, but the complete read sequence is present in at least one of the segments.

Also, there is no unique sign in CIGAR of an individual segment to identify split reads. Each aligned segment of the original read may have long deletion in its CIGAR or aligned pairs. However, without looking at all segments with the read ID it is not clear whether the segment has a novel sequence in it that doesn't align with reference or that part matches different chromosome/coordinates.


Including Picard tool in galaxy I am running the Picard tool CollectSequencingArtifactMetrics in galaxy (locally in my pc).

The output files are saved into disk and not passed to stdout and (I guess) they can't be passed as output from picad to galaxy .

In the xml file I am including:

<outputs>
    <data name='preSummary' format="tabular"  from_work_dir="pref.pre_adapter_summary_metrics"/>
</outputs>

When I click on View details, in "stderr", I see:

  [Mon Mar 11 13:43:07 EDT 2019]
  picard.analysis.artifacts.CollectSequencingArtifactMetrics done.
  Elapsed time: 0.19 minutes. Runtime.totalMemory()=3323461632

1 - It looks like it worked but shouldn't I see the highlighted information in "stdout" rather than in the "stderr"?
2- If it worked, How can I get my output files? I tracked the working directory printing pwd from the <command> but I guess all files are deleted immediately once the picard process has finished.

Any idea of what am I doing wrong?

What basic knowledge am I missing here?

Thank you!

Picard, like most tools, prints logs to stderr.
It should be in the Galaxy history. Are you using a version of picard from the toolshed (if so, which one)? You might set cleanup_job to never in the config, which will prevent galaxy from cleaning the working directory. This will enable seeing what the job actually produced.

Most likely there's more in the log file that you're not showing is us and picard itself never ran to completion. In lieu of that information, what I wrote above is usually sufficient to debug things.


multiple COSMIC id for the same mutation I would like to know why there are different entries for the same mutation in COSMIC, for example COSM6954941 and COSM12833 both refer to ERBB4 c.908C>A.

In this specific case the field Gene name in the two entries is different, COSM6954941 report ERBB4_ENST00000402597 while COSM12833 report ERBB4, but there are many cases in which the Gene name is the same, for example COSM1361 and COSM1219.

Harry Boutselakis addressed the question via the customer support of the COSMIC web site.

  Dear Ivan, Thank you for using COSMIC.

    I would like to know why there are different entries for the same mutation   in COSMIC, for example COSM6954941 and COSM12833 both refer
    to ERBB4   c.908C>A.

    In this specific case the field Gene name in the two entries is different,   COSM6954941 report ERBB4_ENST00000402597 while COSM12833
    report ERBB4,

  The mutations are annotated onto two different isoforms of ERBB4. So
  this annotation refers to two different transcripts. In these cases
  the consequence of the mutation on the mRNA and polypeptide is the
  same reflected by the syntax.

    but there are many cases in which the Gene name is the same, for example   COSM1361 and COSM1219.

  This case the consequence one mutation on the mRNA is unknown , but
  the amino acid change is the same on both cases.

  The information stored in COSMIC is curated from experimental data,
  depending on the evidence and experimental methodology used, the
  mutations are differentiated accordingly. I hope this answers your
  question.

  Kind regards,

  Harry


What are the pros and cons of the different basecallers in Oxford Nanopore Technology Sequencing? What are the pros and cons of the different basecallers in Oxford Nanopore Technology Sequencing?

I am about to start a MinION run on my laptop. What should I consider when choosing my basecaller? Can I let MinION do its sequencing and generate fast5 files and call the bases later? I think I might use Albacore. Metrichor and MinKNOW can also do basecalling. Are there reasons to choose one over another?

There are also third party free and open source basecallers that haven't been developed by Oxford Nanopore.

Of particular note is Chiron, which gave the best uncorrected assembly identity among the base callers that Ryan Wick tested. It's slower than Albacore, but appears to be more customisable, and could theoretically be used to model any DNA feature that is detectable in the electrical signal, rather than just the A/C/G/T base sequence.

There's a paper about Chiron on BiorXiv.


Should you filter GWAS hits with high standard error? I'm trying to figure out if I should be filtering out GWAS hits that have high standard error and I'm not quite sure what to do.  It seems like it might not matter, because the standard error is used to calculate the t-statistic, which is then used to calculate the p-value.  So in a way it's already built in.  But reporting SNPs that have very high standard error doesn't seem quite right.  What's the best way to handle this?

I am curious about the relationship between your p-value (or effect sizes) and standard error. I would expect the significant signals to have smaller stand error compared to the non-significant, background signals. If this is the case, there is no need to do filter.


Hypergeometic test for the overlap of two set of genomic intervals (e.g. from Chip-seq data) Suppose I have done two Chip-seq experiments. Now after the bias correction step such as MACS2, I get two set of genomic intervals (or peaks, as what is usually called) that corresponds to each of the experiment. I want to test the significance of the overlap of these two set genomic intervals that corresponds to each experiment.

The common test employed seemed to me is the Hypergeometic test. There's even a package specifically aim for this task: ChiPpeakAnno.

Recall from the Hypergeometric distrbition -- the analogy is that we have a bag of marbles with two different colors, say one blue and another red. A sample from such distribution consists of marbles drawn from this bag, and a probability measure for this distribution is the number of blue marbles in the sample.

Could someone explain what is the correspondence between the hypergeometric test and the significance overlap between two different set of genomic intervals?

There is no such thing as a hypergeometric test, at least in statistical textbooks. It's a fisher test based on hypergeometric distribution.

If it is chip-seq for the same target, i.e biological replicates, significance of overlap is not quite meaningful. You get more information by for example checking the correlation of the coverage between your replicates using deeptools or anything else.

If you are comparing between two different targets, for example, chip for two different histone marks, you can use the fisher test or chi-square this way.

Partition the genome in X bp bins, these are all your marbles
Identity those that overlap with peaks from target 1 and those that don't
Do 2. for target2
construct a contingency table using 2 and 3.

So you should end up with something like this (purely made up):

                      no_overlap_target2 overlap_target2
  no_overlap_target1               9799             112
  overlap_target1                    88               1

And if there is overlap, the bottom diagonal with be higher.. less on the off-diagonal, and you can test it with a fisher.test.


Using multidimensional scaling to visualize protein sequences by functionality I have a multiple sequence alignment (MSA) of protein sequences, with which I have performed multidimensional scaling to visualize their clustering.  Using R, I used the seqinr package to create an identity matrix (using dist.alignment).  Is there a similar method that is more "generalizable"?  Specifically, I'm interested in classifying the amino acids by some functional characteristics, such as polarity, instead of having the alignment list the amino acids by their identity (A, G, I, etc.).

Here's a very short set of fake input data, with columns 2-4 representing amino acid residue positions:

a <- c("Strain 1", "Strain 2", "Strain 3", "Strain 4")
b <- c("Polar_hydrox", "Nonpolar", "Nonpolar", "Polar_hydrox")
c <- c("Nonpolar", "Nonpolar", "Nonpolar", "Nonpolar")
d <- c("Nonpolar_cyclic", "Nonpolar_cyclic", "Polar_amide", "Polar_amide")
e <- c("Polar_amide", "Acidic", "Polar_amide", "Acidic")
x <- data.frame(a, b, c, d, e)

Is what I'm asking a feasible/correct use of MDS?  Or should I be looking into a different method entirely?  I attempted to use Gower's formula to create a dissimilarity matrix on the original dataset (i.e. using the actual amino acid identities instead of the functional groupings), but compared to what I got using dist.alignment, the MDS plots look different.

Edit: Gower's formula and the dist.alignment function are obviously calculating the distance matrices differently, but I don't know the details on the different methods.  Would using Gower's formula be a good approach for calculating distance matrices for the example dataset?

Edit: clarified functions used for calculating distance matrix.

Some of the first MSA analysis, eg of codon usage of bacteria were performed using matrix factorization (see old papers from Des Higgins). His group also published discriminative methods for finding functional residues using supervised correspondence analysis

  Wallace IM1, Higgins DG BMC Bioinformatics. 2007 Apr 23;8:135.
  Supervised multivariate analysis of sequence groups to identify
  specificity determining residues.

The method is available in the R package ade4 and Bioconductor package made4

Also have a look at the phylogenetic R package adephylo and several related packages


Finding peaks and estimating cell population sizes in multi-dimensional flow cytometry data Our research institute processes a lot of flow cytometry data, but the produced data is under-utilised due to the effort required to process it. A typical run will produce 5 million events (ideally one event per cell), with up to 14 dimensions of [ideally] log-normal fluorescence values for particular groups of cells ("populations"). Due to various systematic errors, negative values, scatter, and a non-zero "zero" value can happen, but I'm going to ignore those for the purpose of this question and assume that the data are well-distributed.

Researchers will typically probe these data using manual filters (set up specifically for each experiment) to find populations of interest. I suppose a picture might help. This one shows three identifiable cell populations, in order of size one at about (X/Ly6C:2,Y/CD86:2), one at about (4.5,3), and one small population at (2,4).

Here's another plot that has two cell populations that are close to each other, such that the population "humps" overlap substantially.

Manual filters are typically used because it can be very difficult to distinguish between a noisy data point from a large cell population and a less noisy data point from a small cell population, particularly when considering populations that make up about 0.01% of the total cells.

As an additional complication, counting these populations can be difficult when populations overlap (as in the second image). A filter/slice through the plot that separates populations could count many cells as being members of the wrong population.

Can these populations be detected in an automated way? If it is assumed that the populations are spread in a gaussian fashion at some point in each dimension, is there some method that can be used to approximate the number of cells in each population, even when populations are close by?

There are quite a few algorithms developed for the automatic classification of multidimensional flow cytometry data, you can see a (not so recent?) review here: http://www.nature.com/nmeth/journal/v10/n3/full/nmeth.2365.html. In your case, you are interested in the unsupervised methods, since you do not have data from single populations (which would provide you "labels" associated to each sample), as I understood.

In other words, if you have your data set as a matrix of 5 million samples (cells) with 14 dimensions (FC values) each, then you are looking basically for algorithms that will do a clustering of the samples into groups. By the way your description of each sample belonging to a multivariate Normal distribution leads to a Gaussian Mixture Model, which can be estimated from your data using e.g. the scikit-learn library for python: http://scikit-learn.org/stable/modules/mixture.html


blasting a refseq protein does not show the protein in the result set Can anyone explain me, why I don't find a specific protein with a blast that was took before from the NCBI refseq database?

Specifically, I was trying to blast the protein with the accession number "NP_420767" and its sequence, respectively, however that protein does not show up in the results. It not only happens when the standard options are chosen, but also, when "Reference proteins (refseq_protein)" as database in the blast options is selected.

I am really puzzled by that.. Shouldn't blast show the initial refseq entry in the results list as it is part of the refseq database and has the same sequence?

NP_420767 is represented by the non-redundant refSeq protein WP_010919826, which has the same amino acid sequence. This is not very clearly annotated, but if you scroll down to the sequence in the GenPept entry for NP_420767, you'll see the following:

CONTIG      join(WP_010919826.1:1..799)

See here for more info.


How many types of DNA methylation are there? My current understanding is:

bases A or C can become methylated due to interaction with a methyltransferase
these methylations only occur if the A or C is surrounded by a specific sequence of base pairs
different types of DNA methylations occur in different species at different rates

Is there a comprehensive list of which sequences of base pairs allow methylations to occur?

I'm currently aware of:

dcm-methylation: CCAGG or CCTGG
dam-methylation: GATC
EcoKI methylation: AACNNNNNNGTGC or GCACNNNNNNGTT
CpG methylation: CG (or CHG or CHH)

If I'm understanding you correctly, by "types of DNA methylation" you mean "nucleotide contexts where DNA methylation occurs".
This is going to be a function of the methyltransferase proteins involved in the methylation process, and this largely depends on the organism, and even cell type being studied.

Generally, what has been observed is:

The phenomenon of DNA methylation (methyl group addition to A/C nucleotides) has been observed in both eukaryotes and prokaryotes
The nucleotide context of DNA methylation differs between the organism, and even the cell type being studied, due to the methyltransferase proteins involved in the methylation process
Methylation of RNA, predominantly in N6 position As

In humans:

The DNMT1 methyltransferase is referred to as the "maintenance" methyltransferase, due to its observation of binding to hemi-methylated CpG sites (methylated C on one strand and not the other)
The DNMT3A/B methyltransferases are referred to as the "de novo" methyltransferases for their affinity of binding to naked DNA
These proteins almost exclusively methylate CpG dinucleotides in human cells
There are cases of non-negligible methylation of cytosines in CHG and CHH nucleotide contexts (brain tissues, for example)

The proteins and their effects on DNA methylation are similar for mice, as well.
Plants tend to have similar nucleotide contexts of methylation, but have much more CHG and CHH methylation than mammals.

In prokaryotes, again, the nucleotide context is dependent on the proteins.
Like you mentioned, we have that:

The Dam methyltransferase methylates the N6 position of the A in GATC
The Dcm methyltransferase methylates the C5 position of the second C in CCAGG and CCTGG
The EcoKI methyltransferase methylates As in AACNNNNNNGTGC and GCACNNNNNNGTT

NEB has a good list of its restriciton enzymes, and how DNA methylation affects their binding/cutting, which also gives some info on the type of methylation in cells.

As for a "comprehensive list", you'll likely have to find some relevant review papers in your organism(s) of interest and compile on yourself.

But also think about how you're phrasing your question.

  Is there a comprehensive list of which sequences of base pairs allow methylations to occur?

I don't think that many papers studying DNA methylation think of the nucleotide context allowing methylation to occur so much as observing the context in which methylation occurs and how it relates to proteins in those organisms.


How to predict stop codons in Illumina reads? I have Illumina MiSeq paired-end reads from 150bp amplicons mapped to my reference genome (> 1000X coverage). These reads have indels that may or may not induce frameshifts. If the indel induces a frameshift (i.e. indel is not a multiple of 3bp), I want to know if that leads to a stop codon later in the read. Is there any package/method that would allow me to do this?

Few extra details

I might be wrong but I think this is different than gene prediction as I would want the tool to scan each read in the correct frame; not try different frames to detect ORFs etc.
I expect mosaicism, so I do not want to compute a single assembly/consensus and predict stop codons in it.

Here is how I imagine such a script/tool would work:

Take the position of the first translated codon of the exon I am interested in (from Ensembl for instance?). Eg. gene abc/exon1 starts at chr1:1000.

Then for each read:

Look at which position the read is mapped to my reference. Eg. read 1 is mapped to chr1:1050.
Compute the position of the first full codon in the correct frame in my read. Eg. start codon at chr1:1000 but my read starts at chr1:1050, the first full translated codon in my read should start at chr1:1052 (I think?).
Split all my read into 3-mers starting at that position and build something like a 3-mer/position dictionnary.
Look for stop codons and report the positions.

Does it exist?

With 150bp reads, you might be able to blastx all the unique sequences you've got against the target proteins.


Is there an easy way to create a summary of a VCF file (v4.1) with structural variations? I got a bunch of vcf files (v4.1) with structural variations of bunch of non-model organisms (i.e. there are no known variants). I found there are quite a some tools to manipulate vcf files like VCFtools, R package vcfR or python library PyVCF. However none of them seems to provide a quick summary, something like (preferably categorised by size as well):

type    count
DEL     x
INS     y
INV     z
....

Is there any tool or a function I overlooked that produces summaries of this style?

I know that vcf file is just a plain text file and if I will dissect REF and ALT columns I should be able to write a script that will do the job, but I hoped that I could avoid to write my own parser.

--- edit ---

So far it seems that only tool that aims to do summaries (@gringer answer) is not working on vcf v4.1. Other tools would provide just partial solution by filtering certain variant type. Therefore I accept my own parser perl/R solutions, till there will be a working tool for stats of vcf with structural variants.

According to the bcftools man page, it is able to produce statistics using the command bcftools stats. Running this myself, the statistics look like what you're asking for:

# This file was produced by bcftools stats (1.2-187-g1a55e45+htslib-1.2.1-256-ga356746) and can be plotted using plot-vcfstats.
# The command line was: bcftools stats  OVLNormalised_STARout_KCCG_called.vcf.gz
#
# Definition of sets:
# ID    [2]id   [3]tab-separated file names
ID  0   OVLNormalised_STARout_KCCG_called.vcf.gz
# SN, Summary numbers:
# SN    [2]id   [3]key  [4]value
SN  0   number of samples:  108
SN  0   number of records:  333
SN  0   number of no-ALTs:  0
SN  0   number of SNPs: 313
SN  0   number of MNPs: 0
SN  0   number of indels:   20
SN  0   number of others:   0
SN  0   number of multiallelic sites:   0
SN  0   number of multiallelic SNP sites:   0
# TSTV, transitions/transversions:
# TSTV  [2]id   [3]ts   [4]tv   [5]ts/tv    [6]ts (1st ALT) [7]tv (1st ALT) [8]ts/tv (1st ALT)
TSTV    0   302 11  27.45   302 11  27.45
# SiS, Singleton stats:
...
# IDD, InDel distribution:
# IDD   [2]id   [3]length (deletions negative)  [4]count
IDD 0   -9  1
IDD 0   -2  4
IDD 0   -1  6
IDD 0   1   4
IDD 0   2   1
IDD 0   3   3
IDD 0   4   1
# ST, Substitution types:
# ST    [2]id   [3]type [4]count
ST  0   A>C 2
ST  0   A>G 78
ST  0   A>T 2
ST  0   C>A 5
ST  0   C>G 0
ST  0   C>T 66
ST  0   G>A 67
ST  0   G>C 0
ST  0   G>T 1
ST  0   T>A 1
ST  0   T>C 91
ST  0   T>G 0
# DP, Depth distribution
# DP    [2]id   [3]bin  [4]number of genotypes  [5]fraction of genotypes (%)    [6]number of sites  [7]fraction of sites (%)
DP  0   >500    0   0.000000    333 100.000000


How to iterate protein sequences using amino acids? I was working on certain program using python and I have been using a protein aligned sequence file in two formats, phylip (.phy) and clustal (.aln).
Example clustal file:

CLUSTAL 2.1 multiple sequence alignment

    Homo            -------------------------------------------MEMQDLTSPHSRLSGSS
    Macaca          -------------------------------------------MEMQDLTSPHSRLSGSS
    Canis           REWSSARPERSKGRRKPVDAAAVSAVQTSQTSSDVAVSSSCRSMEMQDLTSPHSRLSGSS
    Mus             -------------------------------------------MEMQDLTSPHSRLSGSS
    Rattus          -------------------------------------------MEMQDLTSPHSRLSGSS
                                                               *****************

Example phylip file:

     3   120
Apis       MPWIKCLGGG GGKARRGKPL SVSQPIHLLV KSEFEPQCHV FRTKQLRKPR
Branchiost ---------- ---------- ---------- ---------- ----------
Nematostel ---------- ---------- ---------- ---------- ----------

           PCHLCHQAVI KQASCCRVCK YICHKTCEDK GARFRDTSQQ IHKQWQADPA
           ---------- ---------- ---------- ---------- ----------
           ---------- ---------- ---------- ---------- ----------

I want to navigate through the file in a way that ;

If the columns of the file have same amino acid, it saves and prints the amino acid that exists in the column.
I need to know how can I search the columns by using the amino acids, if particular group of amino acids exist in that column or not. For example, if 'STA'or 'HY' or 'FVLIM' or 'NDEQHK' (with no '-') exists in any of the columns, then put a '@' in the file at the end of that column.
Following is the code I have been trying to manipulate the file with the following program:

Code below

alignment = open ("example.aln")
a = b = c = 0
for col in (alignment):
    num = len(set(col))
        if num == 1 and col[0] != '-':
            print (num)
        a += 1
        elif num > 1 and '-' not in col:
        print (num)
        b += 1 elif '-' in col:  # assumes 1 or more dashes
        c += 1
        print('Number of columns with the same amino acid: {}\n'
        'Number of columns with at least 2 amino acids (no gaps): {}\n'
        'Number of columns with one/more gaps: {}'
        .format(a, b, c))

These variables return the occurrences of amino acids in the columns of the file but I don't understand how to search using the amino acids.
What is the best way to iterate over the file and write into a new file with results using python? Should I use python switch statement?

In response to question 1, save the following as a python3 script (e.g., aln_parser.py):

#!/usr/bin/env python3

from Bio import AlignIO
import sys

aln = sys.argv[1]

algnmnt = AlignIO.read(open(aln, 'r'), 'clustal')

for col in range(0, algnmnt.get_alignment_length()):
    proteins = set(algnmnt[:,col])
    if len(proteins) == 1:
        print(f"Postion {col}: protein {''.join(proteins)}")

Run the script to query your alignment (clustal.aln).

python3 aln_parser.py clustal.aln

To save the output to file:

python3 aln_parser.py clustal.aln > results.txt

Requires python3 and Biopython (pip3 install biopython).


How can I extract information from .sam files? I have 10 .sam files after my bowtie2 alignment on ten single-pair sequences. I would like to build a graph based on that output data, however I don't see any other solution but just writing the alignment percentage manually after each iteration.
Is there a tool that can help me to build a .csv file so I can manipulate the output data more easily?

Use samtools flagstat with option -O tsv:

Using -O tsv selects a tab-separated values format that can easily be imported into spreadsheet software. In this format the first column contains the values for QC-passed reads, the second column has the values for QC-failed reads and the third contains the category names.

SEE ALSO:
What Does Samtools Flagstat Results Mean?
INSTALLATION:
Some old versions of samtools do not have option -O tsv. To install a recent versions of samtools, use either:
conda (macOS, Linux) (but see also: Tutorial: installing samtools with conda - overcoming two common errors):
conda create --name samtools samtools

or brew (macOS):
brew install samtools


Influence of pH on protein peptide mass I'm using the Peptides package in R to calculate the mass of peptides with the mw() method. The thing is that I want to calculate them for different pH values, and thus I wrote this little function to take the pH value into account according this image:

My code is this:

phmw = function(pept_seq , unimod_location, pH ){

  monoisotopic_hydrogen_mass = 1.00794

  if(pH < 7){

    # mol_weight = mol_weight_of_peptide + mol_weight_of_Hydrogen
    mol_weight = mw(pept_seq, monoisotopic = T) + monoisotopic_hydrogen_mass

  }else if(pH > 7){

    # mol_weight = mol_weight_of_peptide - mol_weight_of_Hydrogen
    mol_weight = mw(pept_seq, monoisotopic = T) - monoisotopic_hydrogen_mass

  }else if(pH == 7){

    # mol_weight = mol_weight_of_peptide
    mol_weight = mw(pept_seq, monoisotopic = T)
  }

  return(mol_weight)

}

Do you think that this approach is correct or not?

Your approach in general seems to be OK, but I have a few suggestions:

You are basing all your calculations on whether the pH is above or below 7. However the N and C termini have pKa values that determine whether these groups will be ionized (8 and 3, respectively). You should check whether the pH is above/below these values instead of 7.
Different amino acid side chains can also be ionized and should be accounted for when calculating the MW. Amino acid pKas can be found here
Code wise: you pass a variable unimod_location in your function but it's never used. If it's not being used then take it out.


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

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


Question: How to simulate 100k samples having 40 million SNPs in a proportion of case:control=30:70? Note: this question can also be found on Biostars

I need to perform a stress test in a GWAS tool and the duty demands a dataset (plink format) having 100 thousand samples, having 40 million SNPs in a proportion of case:control=30:70. I'm performing the command:

  plink1.9 --simulate ds1.sim --make-bed --out ds1 --simulate-ncases 30000 --simulate-ncontrols 70000

But after about 3 hours execution, I'm facing Error: File write failure., even though my environment has 256T of disk. (It's a container with a S3 (Object store) mounted through s3fs).

I was wondering if there's a clever way to do that? Perhaps breaking it into small chunks and merging up later somehow? It would prevent me to be waiting hours till an error happens and I loose either time and the sample generated so far.

Thanks in advance!

First comment: you should be using real data, rather than simulated data, as we are not very good at modelling genetic variation and recombination at a genome-wide scale. A way that this can be done with GCTA is mentioned in the PLINK manual.

I'll do some calculations to estimate the size of data you need to store:

Assume that this plink-formatted file is in tped format (one row per SNP), and the entire genome is diploid. That means that a single SNP/individual pair will take up three bytes. You want 30,000 cases and 70,000 controls, so each line from each sample will take up 300kB (plus a bit of change for the header). With 1 thousand SNPs, that's 300MB per sample; with 1 million SNPs, that's 300GB per sample; with 10 million SNPs, that's 3TB per sample; with 40 million SNPs, that's 12 TB for one sample.

Reading the plink manual, I see that the BED format is a two-bit compression for each SNP/individual pair, so one byte per 4 individuals, or about 12x compression (i.e. 1TB for one sample).

Okay, that's a lot. It's well within a 256TB storage limit, but might hit file size limits depending on the file system.

If there is a data generation issue with PLINK, it might be worth trying (as you suggest) splitting the input for simulation then merging afterwards. My guess is that doing the simulation on fewer SNPs (rather than fewer individuals) would work better. Plink has a --bmerge option for merging the compressed binary file formats.

I usually approach problems like these by working my way up from something that I know works. Try first merging two datasets of 1000 SNPs. If that works, try larger datasets. Rinse and repeat until the computer crashes, you get bored, or you reach your goal.


Getting structure and sequence related to PDB IDs I have two kinds of interactions: transient and stable. We are supposed to work on stable interaction, like interactions between two monomers in a heterodimer.

In a heterodimer there are two chains in which there are some residues between monomers which take part in the interaction of monomers and building the structure of heterodimer, and there are other residues interacting inside each monomer for building the structure of that monomer. There are some PDB IDs related to heterodimers in PDB.

Now imagine that we are trying to show which residues are physically interacting between two monomers in a heterodimer, and then we are trying two show the site of these physically interacting residues in the sequence related to that special structure. For that, we need to find the sequence related to that structure and we can download it from PDB.

But if we want to make a multiple sequence alignment of that sequence, first we need to find the refseq sequence of that and then we must blast the sequence and find homologues of that sequence.

Here there are some challenges:

we must find residues interacting in the structure of heterodimer between two monomers.
we must find the sequences of monomer chains of the heterodimer and map the interacting residues in the structure to the sequence.
we know that when a structure is solved and its monomers are sequenced, may be they are not able to completely sequence the heterodimer and residue numbers are maybe different from the related refseq sequence.

I would like to know how can we find the related refseq sequence (or sequences, in the case of heterodimer) and then how can we map the physically interacting residues in the structure of PDB ID to refseq sequence (finding the sites of these interacting residues in that refseq sequence) regarding three challenges mentioned above.

Imagine that I have 100 PDB IDs of 100 nonredundant heterodimers and I would like to find and download the structures, sequence and finally refseq sequence. What should I do?

You can download seqs and structures based on a list of PDB ids using http://www.rcsb.org/pdb/download/download.do#FASTA


Time for running ADMIXTURE analysis I am trying to run the ADMIXTURE software for the first time to analyse the structure of an in-house dataset of samples comparing to 1000 genomes project ancestral populations. I am trying to determine ancestry just from the X chromosome of those individuals.
To do that, I have downloaded the X chromosome file from 1000 genomes (ALL.chrX.phase3_shapeit2_mvncall_integrated_v1b.20130502.genotypes.vcf.gz). I have removed duplicated with bcftools and converted to plink binary with the following command in plink 1.9:
plink --bcf ALL.chrX.phase3_shapeit2_mvncall_integrated_v1b.20130502.genotypes.bcf --keep-allele-order --vcf-idspace-to _ --const-fid --allow-extra-chr 0 --split-x b37 no-fail --make-bed

Next, I modified the .bim files to change the polymorphisms ID so them where the same as their position, as a way to standarize their names. Then I filtered related individuals, also with plink and found the intersection of .bim files from 1000 genomes and my in-house dataset (a bed file with microarray genotyped variants only from X non PAR region), that was already filtered to contain only variants genotyped against the plus strand (without indels).
Then I merged both datasets with the following plink command:
plink --bfile filename1 --bmerge filename2.bed filename2.bim filename2.fam --make-bed --out filemerged

That command generated a file with all the SNPS that would have more than 3 alleles in the merged file so I removed all them from both datasets and merged them again. On that merged file I pruned SNPs by linkage disequilibrium:
plink --bfile cromXmerged2 --indep-pairwise 50 10 0.1

, assign sample sex on the .fam file and filtered individual with missingnes higher than 10%:
plink --file mergedfile --mind 0.1 --make-bed --out

I also set all heterozigous calls as missing with --set-hh-missing plink command.
With the final input I runned a supervised ADMIXTURE as following:
./admixture finalinput.bed 3 --haploid="male:23" --supervised -j4

Before running of course I created the finalinput.pop file in open office calc, by creating a list with the assigned population for each sample, with my in-house dataset left blanc. I just used african, european and east asian 1000 genomes samples because my own dataset has brazilian admixed samples. The final dataset has approximately 1500 SNPs and about 1600 individuals. My computer set up is CPU: i3-2370M with two cores, 4 threads@2.4 GHz, RAM 8 GB@1333 MHz and an SSD 120 GB Sata 3.0. OS Ubuntu 20.04 64 Bits.
The thing is that ADMIXTURE have been running for three days without showing any error or producing result yet. I would like to know if it is normal to this analysis to take so long or if it could be any problem with my dataset.
I would appreciate any suggestion or advice!!!
Thanks in advance,
Iriel

Apparently there is a bug in ADMIXTURE that prevents for converging when using at the same time the haploid mode and multithreading (flags --haploid and -j respectively). The problem is solved by avoiding using one of them.


What is Feature Barcoding technology? I heard the term Feature Barcoding a lot so I find resources to understand it. I have read docs of 10xgenomics and I see that Feature Barcoding technology is a method for adding extra channels of information to cells by running single-cell gene expression in parallel with other assays, this technology also do something with t-SNE and something with matrices.

I can't imagine the whole picture because I don't know what are those adding extra channels of information to cells, single-cell gene expression and parallel with other assays mean. And what is the input and output of this technology? Is there any paper go along with this technology?

"Feature Barcoding" refers to using antibodies linked to specific oligos. The antibodies can then bind to cell surface markers of interest and the oligos they're conjugated with turned into normal scRNA-seq reads. This effectively means you can get a FACs profile for each cell without needing to run it through a FACs sorter and individually label it.

The "feature barcoding" nomenclature comes from 10X I think, but there are a few similar protocols (maybe they're exactly the same in the end, I'm not sure about the details) such as CITE-seq and REAP-seq. You can see a talk from 10X on this here.


Question about Co-expression analysis and finding targets for lncRNAs I have a dataset with 88 tumor and 113 normal samples. Among the 50k genes after filtering there are a total of 28k genes (both mRNAs and lncRNAs)

I wanted to do co-expression analysis between lncRNAs and mRNAs to find the targets for each lncRNA. I'm thinking to use WGCNA for co-expression network analysis.

Few doubts:

1) For Co-expression network analysis (WGCNA) should I use only tumor samples (or) both tumor and normal samples?

2) If not WGCNA, if I use cor function in R with method Pearson, on what cutoff I should select the target genes?

3) Heard that with bedtools we can know which genes are 10kb upstream/ downstream of lncRNAs to find neighboring pc genes. But no idea how to do that. Could you please give me an example for this.

Thank you

You have several options to approach this with WGCNA (weighted correlation network analysis). You can run a WGCNA on the combined set, identify modules and select those lncRNAs for further follow-up that have relatively high intra-modular connectivity (i.e., are intramodular hub genes). I would suggest this approach first.

Second option is to run WGCNA on mRNA for coding genes only and screen module eigengenes and lncRNAs for association with one another. This approach is usually helpful when the two data types (mRNA and lncRNAs) are different which I suppose is not really your case (mRNA and lncRNAs were probably measured together and are overall similar).

I would start using all samples together (tumor and normal). Hopefully you'll find some interesting and plausible information; in the second pass I'd perhaps create separate (or consensus) networks for tumor and normal samples.

I am not familiar with bedtools, but you should have some sort of gene/transcript annotation with your data (it was surely aligned to a reference genome and summarized to some sort of annotated set of genes). The annotation should contain chromosome and base pair position, which would make it rather easy to find genes within a specific interval.


How to deal with duplicate genes having different expression values? I have RNA-Seq data which is FPKM. In the dataframe df first column is gene_name and the other 100 columns are samples.

Usually if it is counts data I do like following:

df2 <- aggregate(. ~ gene_name, data = df, max)

I'm not sure what do with the FPKM data if there are duplicate genes with different FPKM value for the same sample.

Lets say it looks like this:

gene_name     sample1        sample2        sample3
5S_rRNA      0.3206147    0.3327312      0.377578
5S_rRNA      0.3342000    0.0000000      0.1305166

Any suggestions please.

I assume you're familiar with the various issues surrounding FPKMs, so I'll not expound upon them.

As a general rule, you should be using gene IDs rather than gene names, since the former are unique while the latter are not. If you only have access to data quantified on gene names, then the appropriate way to merge RPKMs is with a weighted sum:

$FPKM_{gene} = \frac{FPKM_{copy1} * Length_{copy1} + FPKM_{copy2} * Length_{copy2}}{Length_{copy1} + Length_{copy2}}$

As an aside, rRNA expression levels will tend to be wrong, since they're normally excluded by either poly-A enrichment or the use of ribo-zero.


How to plot p-values in a circular barplot? Im trying to plot reactome pathway in for of circular bar plot using ggplot2

my data

    dput(df2)
structure(list(ID = 1:10, Pathway_names = structure(c(5L, 1L,
7L, 6L, 2L, 3L, 9L, 8L, 10L, 4L), .Label = c("Antigen Presentation: Folding assembly and peptide loading of class I MHC",
"Antigen processing-Cross presentation", "Class I MHC mediated antigen processing & presentation",
"Cytokine Signaling in Immune system", "Endosomal/Vacuolar pathway",
"ER-Phagosome pathway", "Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell",
"Interferon alpha/beta signaling", "Interferon gamma signaling",
"Interferon Signaling"), class = "factor"), Entities_found = structure(c(1L,
2L, 3L, 4L, 4L, 5L, 6L, 7L, 8L, 9L), .Label = c("51", "54", "55",
"57", "75", "77", "96", "117", "176"), class = "factor"), Entities_pValue = c(1.11e-16,
1.11e-16, 0.003908431, 2.54e-12, 2.77e-10, 0.0053068, 1.64e-13,
1.11e-16, 1.11e-16, 3.68e-06)), row.names = c(NA, -10L), class = "data.frame")

my code

ggplot(df2, aes(x = factor(Pathway_names), y = Entities_found)) + ######################################
  geom_bar(aes(fill = Entities_pValue), stat = "identity", position = "stack") +
  #scale_fill_manual(values = c(blank = "white", dark = "black")) +
  coord_polar() +
  theme_minimal() +
  guides(fill = FALSE)

So far im not sure how to include both p value and entities in the plot so that the pathway that are coming up can be labeled accordingly

New figure

The adjustcolor function from grDevices can be used to generate colours with transparency. What you want is colours based on the pathway and transparency to form a colour gradient for the p-value.

adjustcolor( "red", alpha.f = 0.2)

For the ggplot example, the aes function has an alpha parameter.

ggplot(df2, aes(x = factor(Pathway_names), y = Entities_found), alpha = Entities_pValue) + ######################################
geom_bar(aes(fill = Entities_found, alpha = Entities_pValue), stat = "identity", position = "stack") +
coord_polar() +
theme_minimal() +
 guides(fill = FALSE)


To find p-value of compounds (pairwise) I have next SNPs data.frame structure

> str(SNPs)
'data.frame':   1703 obs. of  4 variables:
 $ group: Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ rs1  : Factor w/ 3 levels "D/D","I/D","I/I": 1 1 2 3 3 2 1 1 1 1 ...
 $ rs2  : Factor w/ 3 levels "a/a","a/b","b/b": 3 3 2 3 3 2 2 3 3 2 ...
 $ rs3  : Factor w/ 3 levels "G/G","G/T","T/T": 2 1 2 1 1 3 1 1 2 1 ...
 ...other rs

> head(SNPs)
  group rs1 rs2 rs3 ...other rs
1     A D/D b/b G/T
2     A D/D b/b G/G
3     A I/D a/b G/T
4     A I/I b/b G/G
5     A I/I b/b G/G
6     A I/D a/b T/T

For example, I noticed that rs5 4b/4a and rs6 G/G in group A very often occur together (see below). In group B they very often occur together too. So I want to know - is it statistical regularity or not.

I can create table with all pairs in both groups

> SNPs$rs5_rs6 <- paste(SNPs$rs5, SNPs$rs6)
> tmp <- table(SNPs$rs5_rs6, SNPs$group)
> tmp

              A   B
  4a/4a G/G   1  20
  4b/4a G/G  31  83
  4b/4a G/T  14  51
  4b/4a T/T   1   0
  4b/4b G/G  37 106
  4b/4b G/T  35 119
  4b/4b T/T  11  31

So, now I need compare (find p-value) group A and group B: 4a/4a G/G in group A vs 4a/4a G/G in group B, 4b/4a G/G in group A vs 4b/4a G/G in group B, 4b/4a G/T in group A vs 4b/4a G/T in group B and etc. Only this pairs/

How can I do that?
I need p-values for all

> tmp

              A   B
  4a/4a G/G   1  20 - p-value?
  4b/4a G/G  31  83 - p-value?
  4b/4a G/T  14  51 - p-value?
  4b/4a T/T   1   0 - p-value?
  4b/4b G/G  37 106 - p-value?
  4b/4b G/T  35 119 - p-value?
  4b/4b T/T  11  31 - p-value?

You can get standardised residuals from chisq.test and convert these to p-values, but it's not recommended. The frequencies of each genotype depend on the counts of all other genotypes, so the p-values aren't independent and they aren't corrected for multiple testing, but if you must:

pnorm(-abs(chisq.test(tmp[,c('A','B')])$stdres)))


How to download RNAseq gene expression data from GTEx I would like to download RNAseq gene expression data from GTEx from the liver only. What is the best way to do this?

I have tried a few things. I downloaded several files from the GTEx website (https://www.gtexportal.org/home/datasets, under "RNA-Seq Data"), but the problems I run into are:

Not sure what is the best file: "Gene read counts" or "Gene TPMs"? What is the difference?
I would like 1 value per sample per gene. There are many duplicate genes in these 2 files. How can I collapse them to 1 value per gene?
How can I download the data for liver tissue only, instead of ALL data for ALL tissues?

Thank you for suggestions.

Gene read counts are the number of reads or fragments that are mapped to each gene. So it is an integer and the values are not directly comparable between samples due to differences in sequencing depth. TPM stands transcripts per million. It is the normalized gene expression level. Basically, you first normalize read counts by gene length and then normalize the library size to 1e6 in all libraries. As for which one should you use, it depends on what analysis you want to do and what tool you are going to use.
I have not noticed that there are duplicated genes in those files, can you show an example?
I do not think you can download expression from one tissue, but you can always subset. The whole dataset is not very large.


When performing differential expression analysis, should genes with low read counts be removed before or after normalization? I have RNA seq data which I've quantified using Kallisto.  I'd like to use tximport to transform the read count data into input for EdgeR, following the R code supplied in the tximport documentation:

cts <- txi$counts
normMat <- txi$length
normMat <- normMat/exp(rowMeans(log(normMat)))
library(edgeR)
o <- log(calcNormFactors(cts/normMat)) + log(colSums(cts/normMat))
y <- DGEList(cts)
y$offset <- t(t(log(normMat)) + o)
# y is now ready for estimate dispersion functions see edgeR User's Guide

I would ideally follow this up with EdgeR's exactTest.  However, I'd also like to remove genes with low read counts using code such as this:

filterCounts <- function(counts) {
  cpms <- cpm(counts)
  keep <- rowSums(cpms > 1) >= 1
  counts[keep,]
}

The only issue is that I can't tell whether I should remove the genes with low read counts before or after I normality the counts while importing with tximport in the first code.  I'm thinking that it would be best to do this before, editing out the abundance, counts, and length rows of the txi data frame which correspond to lowly expressed genes, but I would like this to be confirmed before I proceed.

The more genes you have the more robust the scaling factor is (among the reason why one doesn't normalize to ERCC spike-ins without a compelling reason), so I suppose in theory it's better to filter after determining the scale factor. Having said that, I'd be surprised if the results changed much either way. Unless you end up filtering out a LOT of genes the scaling factors should be fairly robust to the change.


Can you query the RCSB PDB API with a ligand and return the IDs of all entries in which it appears as a free ligand? If you're looking at the entry for a ligand, for example, adenosine, there's a link on the webpage to entries for which this entry exists as a free ligand.
However, I want to access this via the API and can't seem to find docs/info on how to do that.

I have been looking at the fields you can use to get custom reports.

I'm not sure about the RCSB, but you can get such a list from the PDBe REST API, as documented in the compounds section:

http://www.ebi.ac.uk/pdbe/api/pdb/compound/in_pdb/ADN

Alternatively, from an older service:

http://www.ebi.ac.uk/pdbe-srv/pdbechem/PDBEntry/download/ADN


How to calculate RNA copies in qPCR using R? I have these Ct data from qPCR for a series of samples A-H, I want to quantify the RNA copies of my samples using R:
Sample_ID   Replicate1_Ct Replicate2_Ct
A           22.92         24.21
B           24.29         22.16
C           24.29         22.16
D           27.91         23.76
E           21.44         22.81
F           16.50         18.31
G           17.79         16.56
H           17.34         19.09

In addition, I have measured a standard (STD) curve of copy number of the RNA in question, that I want to use as a reference:
STD_Ct      STD_copies      Log10(STD_copies)
33          4.00E+03        3
28          4.00E+04        4
24          4.00E+05        5
20          4.00E+06        6
15          4.00E+07        7
12          4.00E+08        8

It seems
This the original data form:

Now that we've straightened out the data, it seems like the simplest solution is to just linearly interpolate the copy number from your standard curve.
When I read your standard curve (std) and your sample data (ct) into R, I see that the log10(Ct) fits the copy number fairly well, whereas the linear is a bit harder:

So we can fit a model to predict copy number based on Ct on the STD data:
> fit = lm(std$Log10.STD_copies. ~ std$STD_Ct)
> fit

Call:
lm(formula = std$Log10.STD_copies. ~ std$STD_Ct)

Coefficients:
(Intercept)   std$STD_Ct
    10.6847      -0.2357

# we plug those numbers into the linear equation,
# taking the mean of the two replicates for a best estimate.

> rowMeans(ct[,2:3]) * -.2357 + 10.6847
[1] 5.130429 5.210567 5.210567 4.595390 5.469837 6.582341
[7] 6.636553 6.391424

So those are log10(copy number) estimates, of a sort. You can get much more elaborate, of course.


How to create Seurat object while RNA expression and ADT combined into one matrix I got an input matrix which RNA expression and ADT capture are combined into one file. I loaded the file into Seurat successfully, however, when I tried to create Seurat object, it threw out an error saying

Error in CreateAssayObject(counts = counts, min.cells = min.cells, min.features = min.features) :
  No cell names (colnames) names present in the input matrix

How can I create Seurat object when RNA expression and ADT data are combine together?

Given the fact that you have RNA and ADT data (probably on the same set of cells), some colnames (which represent individual cells) would be duplicated. Normally duplicated colnames are tolerated with matrices and sparse matrices, however, Seurat apparently does not do so, below is from the ?CreateAssayObject:

Details
Non-unique cell or feature names are not allowed. Please make unique before calling this function.

Possible solution to your problem: Seuart has a dedicated vignette for working with multimodal data and as you would see you will need to initiate your Seurat object with one matrix per assay: RNA and ADT. All you need to do is split your matrix into RNA and ADT, create your Seurat object with RNA data and then add the ADT data with:

seurat_obj_with_rna_only[["ADT"]] <- CreateAssayObject(counts = your_adt_matrix)

For efficiency, Seurat uses sparse matrices so don't forget to convert your data matrices to sparse.


Working with old genome builds Is working with and relying on old genome builds still valid?

For example NCBI36/hg18. Would results from papers based on old builds require LiftOver and re-analysis to be useful?

A bit of context, this is related to other post, where we have aCGH results based on old build: How do I validate a single sample ArrayCGH result?

In my opinion, it is not very reliable. LiftOver is very limited in terms of transformations it can support. The LiftOver Chain format can capture only matching regions in the same order. It means that it can account for indels, but even simple structural variations become problematic.

For instance, when a newer assembly is available, it is usually a recommended practice to remap all the reads rather than transform the existing alignments.


Cannot install chromosomer I am trying to install chromosomer but I fail. Can anybody help me, please?

$ pip install chromosomer

Collecting chromosomer
  Could not find a version that satisfies the requirement chromosomer (from versions: )
No matching distribution found for chromosomer

The output python:

$ python

Python 3.6.0 |Anaconda custom (x86_64)| (default, Dec 23 2016, 13:19:00)
[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)] on darwin
Type "help", "copyright", "credits" or "license" for more information.

Chromosomer only exists for Python 2. You should thus be able to install it via

pip2 install chromosomer

But:

Your Python installation is a bit screwed up: there are Python 2 and Python 3, which are unfortunately incompatible.

On your system python and pip seem to be aliases for Python 3, which I’d strongly recommend against (lots of tools will break). Instead, they should alias Python 2. Python 3 should be invoked via python3/pip3.


how to get access to genotypes and phenotypes used for a GWAS I'm a master's student working on genomic prediction of complex traits using deep learning.
i'm looking for a dataset of human genotypes and phenotypes that has been used for a GWAS.
The only thing i can get from the gwas i found is the summary statistics not the raw data.
Does anyone know if it's possible to get such data (for free ) for a research project?

Thanks.

How about:

dbGaP - https://www.ncbi.nlm.nih.gov/gap/ - in case you are searching for supporting datasets I leave its entry in BiŌkeanós

https://www.ddbj.nig.ac.jp/jga/ - the access might be a bit more complicated.

https://www.internationalgenome.org, not particularly large, in their FAQ they state that the individuals sequenced have not self-reported any diseases - perhaps you might want to add those with to your negative set.  Again, its entry in BiŌkeanós.

If you are looking for more comprehensive information, you could try biobanks - e.g. https://www.ukbiobank.ac.uk. These are usually open (or reasonably priced) for research but the approval process can take time.
A side note - some of GWAS catalogues will have a link to the study where the association originated from, e.g. here. You can sometimes find the datasource there.


COSMIC Genotypes and Phenotypes I'm trying to find files containing genotypes and phenotypes for cell lines in COSMIC.  On their downloads page I found links for mutation data, which I guess could be converted to genotypes, but I haven't been able to find any phenotype information.  I'm specifically interested in drug sensitivity information.  Does anyone know if it's possible to get more user friendly genotype files and whether phenotype data exists publicly?

It seems like I was wrong and COSMIC does store information about cell lines.

You need to download the files from Genotypes:

  Genotypes

  Files listing the SNP calls for each cell line identified by PICNIC
  analysis of Affymetrix SNP6.0 array data. Both a simple genotype (AA,
  BB - homozygous or AB - heterozygous) and a complex interpretation of
  the genotype are given (for example, in a triploid region of the
  genome the genotype maybe AAB). 'genotypes/.'

You need to register (and get a license) to be able to download the data. Presumably the cell lines are described in the heading of the file. Then you can sort them out by cell line.

Then once you know the germ line identifier look for the paper describing it to know more about their phenotype. (But maybe there is a database with phenotype of known cell lines)