text:
how to extract a different substring for each row of a dataframe in R?
I have a dataframe with putative alignments of thousands of probe sequences. I would like to pull out the adjacent nucleotide from the genome for each alignment. My dataframe includes seq.name for the target chromosome and I've added an adjacent_pos variable with the location on the target chromosome that I want. Using read.fasta from the phylotools package gives me a genome in the format of a dataframe with seq.name and seq.text for each chromosome and unplaced scaffold.
Working on a test data set with a single chromosome, I successfully used something like this:
adj <- cbind(alignments$adjacent_pos, alignments$adjacent_pos)
alignments$adjacent_base <- str_sub(genome$seq.text, adj)

Now back to the whole genome, I can't wrap my head around how to direct str_sub to the correct seq.name for each seq.position. I'm sure there's an elegant solution, but I haven't been able to find it!

please can you provide a small reproducible example for your alignments and genome objects? the merge() function from the data.table package will surely help you.

Say your file looks like this:
cat test.fa
>chr01
ACCGGTT
>chr02
ATGCATGC

Something like what you have done, to read in the fasta and extract the sequence. Key part is creating a GenomicRanges object:
library(GenomicRanges)
library(BSgenome)
library(Biostrings)

adj <- data.frame(chr=c("chr01","chr02"),start=c(2,3),end=c(2,3))

    chr start end
1 chr01     2   2
2 chr02     3   3

adj = makeGRangesFromDataFrame(adj)

fasta <- readDNAStringSet("test.fa")
SEQ = getSeq(fasta,makeGRangesFromDataFrame(adj))

as.character(SEQ)
[1] "C" "G"

or if the fasta file is huge, you can use Rsamtools, indexing the fasta file first:
library(Rsamtools)
fasta = FaFile("test.fa")
indexFa('test.fa')

as.character(getSeq(FaFile("test.fa"),adj))
chr01 chr02
  "C"   "G"

but it's not just one location from each chromosome. I have about 160k alignments, each needing a different location pulled out from a 3G genome that has 21 chromosomes plus over 7500 unplaced scaffolds. Still perhaps the GenomicRanges object is the way to go?

yes.. the code above works for as many locations as you have. You just need to construct the genomic ranges object like I did above, and I suggest using the Rsamtools solution if your file is huge

adj should be GenomicRanges object...if it was clear from the solution above

thank you I will try

brilliant, worked like a charm! thank you

you're welcome :) glad it worked for you


text:
how to easily make complex Euler/Venn diagrams in R?
From my microarray data I have selected three subsets of probes. I would like to display these subsets as a Euler diagram to visualize the overlaps, but am struggling with inputting the combinations into Eulerr. All three subsets need to be nested within the full set of probes; one subset is fully contained within another, and the third overlaps with the other two. Is there a package into which I could feed the vectors of probe names and let it determine the overlaps for me? Many thanks in advance!

Since you are in R, consider using the UpsetR library.  It doesn't make Venn diagrams but it helps to visualize overlaps between any number of groups.
http://gehlenborglab.org/research/projects/upsetr/#:~:text=UpSetR%20is%20an%20R%20package,based%20on%20groupings%20and%20queries.

You might consider different ways of using the venneuler package. I'm not super thrilled with its usability (or with the accuracy of the overlaps), but you can get it to make custom Venn diagrams.
Here is a thread with some suggestions along these lines.
Below is some sample R code that I slightly modified from an old project that I used to make a 3-category Venn diagram.
I didn't run this to make sure that it works because of an irritating Java dependency issue on my OSX laptop that I don't feel like fixing right now, but as long as you have a JDK on your machine and can install the packages it should work.
# note irritating rJava dependency of venneuler
library(plotrix)
library(venneuler)

#A-B sharing: 200414
#B-C sharing: 17561
#A-C sharing: 16764
#All share 14201

# get coordinates for circles- replace as needed
# this function is supposed to plot stuff, but doesn't actually
ven = venneuler(c(C=48147-17561-16764+14201,
                  B=347181-17561-200414+14201,
                  A=325971-200414-16764+14201,
                  "C&B"=17561, "B&A"=200414, "C&A"=16764, "A&B&C"=14201))

cats = c('A','B','C')

# note- this following only shows up in the .pdf printer, not in the X11 utility plotter
pdf('venn.pdf')
plot(c(-.5,1.5), c(-.5,1.5), type='n', xaxt='n', yaxt='n', bty='n', xlab='', ylab='')
draw.circle(ven$center[1,1], ven$center[1,2], ven$diameters[1],
        col = hsv(1,1,1,.5))
draw.circle(ven$center[2,1], ven$center[2,2], ven$diameters[2],
            col = hsv(.5,1,1,.5))
draw.circle(ven$center[3,1], ven$center[3,2], ven$diameters[3],
            col = hsv(.1,1,1,.5))

text(ven$center[1,1]+.065, ven$center[1,2], labels=cats[3], cex=1.3)
text(ven$center[2,1], ven$center[2,2]-.2, labels=cats[2], cex=1.3)
text(ven$center[3,1], ven$center[3,2]+.2, labels=cats[1], cex=1.3)

dev.off()

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
How to calculate Gene Ontology terms in python
I am testing different Protein-Protein Interaction networks alignment tools. I observed that all alignment tools show some Gene Ontology terms similarity. How can I calculate Functional similarity / average functional similarity for Gene ontology's terms Biological processes, Cellular components, Molecular functionality of any two proteins networks

Welcome.  What tools have you been using?

hubalign ......

Do you need it in python? I know one very nice package in R that does that.

@Llopis it would be nice in python but if any R package you know, kindly mention that

You can calculate gene ontology similarity with the GOSemSim R package (paper)(disclaimer I'm a contributor to this package). You have several similarity scores implemented, some of them are across the three subontologies.
In python you have the GOATOOLS (paper). I think there was another Python package but now I don't remember its name.

thanks @Llopis. i looked into goatools already but i don't think that they are have methods to calculate GO CC, MF or BP

Please, mention in the questions what you have done, because now I feel that I wasted my time trying to answer you when your problem is something else. Did you look at the notebook ?

sorry if you feel that way. but did you saw any GO calculated in goatools ?

I linked a notebook explaining how they do it..

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How to find the position-specific weight matrix (PSWM) for a list of genes having different sizes?
Wiki has the algorithm for finding position-specific weight matrix (PSWM) for a list of nucleotide sequences having an equal number of nucleotide. How do we find the PSSM for a list of nucleotide sequences having different sizes?

If your sequences are different lengths, this implies a gapped alignment with insertions/deletions. This is not possible with a regular PSSM, so you have to create a separate PSSM for each group of homologous sequences of the same sequence length.
The alternative is to use a different approach that allows gaps like profile HMMs (as implemented in HMMER3 or HH-Suite3).


text:
Colouring amino acid substitutions by scoring matrix
I have a list of site-specific amino acid substitutions from a mutagenesis experiment which I am visualising by colouring the structure of the protein in PyMOL. I would like a way to colour the substitutions by some approximation of the biochemical (dis)similarity between them.
So far, I have tried to create a colour gradient of scores from an alignment scoring matrix (e.g. BLOSUM62, therefore a 16-shade white to red gradient; a substitution of A to R for example gets a -1). However, I am fairly certain this is not the best interpretation of of the log-odds scores in the BLOSUM matrix, and I wondered if anyone had any thoughts on how I might better accomplish this?
Note that other colour schemes such as hydrophobicity, Tailor etc. aren't suitable for my purpose.

If you wanted to plot the diversity of an MSA Cosurf (https://consurf.tau.ac.il/) would the tool to go for (consurfDB is the precalculated version). For a set of point mutations and using PyMOL only (ouch!), I would suggest checking this PyMOL snippet, which gives a rough ∆∆G, which could be set as a b factor (alter resi xx, b=1.23) as played with that way.

*) I mention ∆∆G because that would be a metric of change in stability. Whereas a substitution matrix is a metric of how uncommon such a change is from a selection point of view. If your experiment is a Lenski style experiment (long term evolution experiment), that makes sense, but not if it's directed evolution following in vitro mutagenesis (mutazyme, magnanese, wobbly analogues etc) even if with in vivo selection (e.g. phage display, biosensor etc etc.) as the mutational spectrum is a product of the in vitro step not biological selection of a population.

Hi @MatteoFerla many thanks for your thoughts. I think the Gibbs energy approach is maybe not quite what I am looking for because we know that both the start and the end point fold and are active - we are particularly interested in IDing important substitutions for changes in activity between the two isoforms (this being a somewhat unguided experiment, based on ancestral sequence reconstruction)

Ah. That makes sense if it's an ASR —plus you probably have some rad epistatic mutations, which the ∆∆G of each in isolation would have been dismal. A substitution matrix is poor metric for dissimilarity because it would just light up the surface while in reality the slow evolving core mutations are the interesting ones. Consurf would reflect this, but is not specific for your changes. Is surface/core the reason why hydrophobicity, volume etc etc. are not suitable?

Indeed. So ideally we would like some way to winnow this (obviously other than exhaustive combinatorial mutagenesis, which we have in the pipeline).


text:
Why does BWA MEM orientation contradict my library prep method
I have some RNA-seq data from a stranded paired end library prep, with dUTP and UDG preparation, so the orientation should be RF (confirmed with sequencing provider). I assembled the reads with Trinity, specifying RF for the strandedness, and I am mapping the reads to the transcripts with bwa-mem (I am aware bowtie2 or STAR might be better for this, but bwa is what I presently have installed, so just using it to check).
In the bwa output, I see lines # candidate unique pairs for (FF, FR, RF, RR): (0, 5372, 2, 0), which seems to contradict the expected orientation. Is there a reason that this is expected, or does it indicate a problem?

Not familiar with the protocol, but I guess RF means different things here. If I am right, with your RF protocol, the first read comes from the reverse strand of the transcript and the second read from the forward strand, but the pair orientation is still FR:
------------------->  transcript
  R2 ---->   <---- R1

The RF pair orientation should be
<----   ---->

You can confirm this by checking the mapping strand of read1.


text:
How can I get a list of papers where a particular mutation has been mentioned?
I have recently tried a few tools (MutationFinder, tmVar, PubTator etc.) that extract mutations from a given text. However, for my work I need something that will take a list of mutations as an input and return a list of pubmed id in which the mutations has been found.
Example input:
mut_list = [Y702F, L483Y]
Example output:
pub_list = [12215417, 18537569]
Can anyone give any idea?
Thanks!

Welcome to the site. Could you elaborate more your question? What did you try to solve this problem?

You could use the pubmed API to query it directly.  Typically I use Biopython for my pubmed queries.  There are some good examples in the cookbook and on other SO posts.
This example from the cookbook link above has the guts of what you need.  You will just have to iterate over your list.  Depending on the number of queries you want to do you might have issues getting throttled by NCBI
In this example, we will query PubMed for all articles having to do with orchids (see section 2.3 for our motivation). We first check how many of such articles there are:

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.egquery(term="orchid")
>>> record = Entrez.read(handle)
>>> for row in record["eGQueryResult"]: ...
        if row["DbName"]=="pubmed": ...
            print(row["Count"]) 463

Now we use the Bio.Entrez.efetch function to download the PubMed IDs of these 463 articles:

>>> from Bio import Entrez
>>> handle = Entrez.esearch(db="pubmed", term="orchid", retmax=463)
>>> record = Entrez.read(handle)
>>> handle.close()
>>> idlist = record["IdList"]

This returns a Python list containing all of the PubMed IDs of articles related to orchids

A follow-up question: What will be the value of the db parameter if I want to search all databases? I read through Biopython's Entrez.esearch documentation as well as source code but wasn't able to figure it out.


text:
Tumor sample and heterozygous SNVs: what is the genotype of the normal cells?
In the following picture, you can see that in normal cells (yellow) we have two copies of each gene, while in tumor cells (violet), we just have one copy due to a clonal mono-allelic deletion. I do not really care about the genotype of the two genes, I just need to know that they are heterozygous.

On the other hand, suppose we are talking about somatic heterozygous SNVs and therefore we do not consider deletions anymore. In this case I would like to know the genotype (AA, AB or BB) of tumor and normal cells. Based on my understanding, I expect the AF of the SNVs in the tumor cells to be roughly 0.5 (since I ask for heterozygous SNVs). But what about normal cells? What is their genotype? If I ask for heterozygous SNVs, does that imply that normal cells are homozygous for that alleles?

Normal cells also have SNV and mutations, in fact recognizing which mutations and SNV are the ones that cause the tumor is an active research area. But it is impossible to say which is their genotype, as you would need to sequence them to know it.

Your initial information/image and actual question seem pretty unrelated to me, but I'll do what I can.
If you have a heterozygous SNV in a malignant population of cells, it really depends on when the mutation occurred for you to know the genotype of wildtype cells. It could have been there since birth (a germline mutation passed down from the parents) or have been acquired throughout the person's life and clonally expanded (somatic mutation).
Unfortunately, you have no way to differentiate between these two types of SNVs just by sequencing the tumor cells. If you want to know true somatic mutations that are found only in the tumor cells, you need to sequence a matched normal control from the same patient (preferably from the same cell type). This is why most large cancer sequencing studies include not only tumor samples, but normal cells from the same patient - see TCGA or ICGC for examples.


text:
Use of heterozygous SNPs in cancer research: why?
When reading about allelic fraction (AF) and SNPs in cancer research, they always mention the fact that they're using heterozygous SNPs (informative SNPs). Why is this? Why can't we use homozygous SNPs?
In this case, what is the reference base (from the reference genome)?

I’m no longer working in tumour sequencing so I’m by no means an expert. But in a nutshell, the reason is that, as indicated, homozygous SNPs aren’t informative: if your allelic fraction is 100%, meaning you only observe a single nucleotide at a given position, we don’t know whether we’re dealing with a tumour related mutation. In fact, we probably don’t because the sample likely contains healthy as well as tumour cells.
By contrast, tumour-specific variants always have an AF < 100%, in other words they exhibit some degree of heterozygosity.
In fact, in healthy tissues you’d expect to see one of three cases: your variant allele (= sequence base ≠ reference) has a frequency of either 0% or 100% (homozygous), or 50% (heterozygous), ignoring rare somatic mutations. With tumour samples, by contrast, the fraction will vary depending on the penetrance of the variant in the tumour, and on the fraction of healthy to tumour cells in the sample.

I'm not sure what I'm talking about, but I guess one could re-phrase Konrad's answer in terms of information theory.
The information content of a message (a packet of bits) is inversely proportional to the frequency of those bits. For example, if you randomly extract a word (the message) from an English book and the word is "the", you don't get much information about the book because "the" is very frequent and you see it everywhere. If, instead, the word happens to be the much rarer "bioinformatics" you start having a pretty good guess what the book is about.
About SNPs, homozygote reference genotypes have very little information because most of the genomic positions in most of the people are homozygote reference. Heterozygote and homozygote non-reference are rarer and thus more informative.
So, I don't think there is anything wrong with using homozygote SNPs. They are ignored because they are so common that their information content is very small.


text:
How many reads do I need to cover the entire genome?
Suppose my genome is 3 million bases and that my reads are 100 nucleotide long. I need to know how many reads I need to cover the entire genome.
I start from using the equation $C = \frac{N \cdot L}{G}$ where C is the coverage, N the number of reads, L the length of a read and G the length of the haploid genome. I also know that the average coverage follows a Poisson distribution, so I need this for C.
On my slides, in fact, I have the following:
$$
C = \frac{N \cdot L}{G} \approx ln \left( \frac{G}{L \cdot \epsilon}\right),
$$
if we assume that this is the average coverage needed to cover the entire genome with probability $1-\epsilon$.
The problem is that I do not understand how to arrive at the expression on the right. I thought that if we want to be sure to have a coverage of at least 1 on average, we would have had to compute $1 - P(X=0) = 1 - e^{-\frac{N \cdot L}{G}} = 0.99$, if we take $\epsilon=0.01$. But I do not get the same result as with the equation above.

Let's see what your slides are claiming $\epsilon$ to be, using:
$$ C = \frac{NL}{G} \approx \ln \left( \frac{G}{L\epsilon} \right) $$
We can rearrange to get
$$ \epsilon \approx \frac{G}{L} \times e^{-C} $$
Now, what does it mean for an entire genome to be covered? Having a clear understanding of this definition is crucial, and your misunderstanding seems to partially stem from this.
I claim it's natural to consider a genome as consisting of $G/L$ bins. In this setting, we have $N$ balls that each randomly goes into a bin. Then, the entire genome is covered when all bins contain at least $1$ ball/read.
From the Poisson distribution, we know that P(a bin has 0 balls) = $e^{-C}$.
Then, we define $\epsilon$ to be the probability that at least one of the $G/L$ bins has 0 balls. Given the way we've modeled the problem, this is exactly equal to
$$ \epsilon = 1 - (1 - e^{-C})^{G/L}. $$
From above, your slides are claiming $ \epsilon \approx \frac{G}{L} \times e^{-C}$. This might be strange, as it appears to be equal to the expected number of bins with 0 balls, instead of the probability that there is at least 1 bin with 0 balls. Furthermore, it's strange because it contains only one exponential, compared to our exact derivation which has a double exponential.
The mismatch in the number of exponential is a hint that another mathematical approximation is at play. We can use a first-order approximation of the binomial series to obtain
$$ (1-e^{-C})^{G/L} \approx 1 - (G/L)e^{-C} $$ which holds particularly when $|e^{-C}| \times \frac{G}{L} << 1$.
In your case, $\frac{G}{L} = 30,000$, so the approximation only really starts getting accurate when $C > 15$ or so.
Following through with our approximation, we get
$$ \epsilon \approx 1 - (1 - (G/L) e^{-C}) $$
$$ \epsilon \approx \frac{G}{L} \times e^{-C}. $$


text:
De novo motif discovery in protein sequences
I am trying to build a features matrix to be used for Random Forest based classification. I'd like to add, as features, short motifs which are common to all the protein sequences belonging to a specific gene family.
I tried to use MEME but I'd like to know if there are also other good tools (online or command line). Keep in mind I have more than 10.000 sequences and that I prefer short motifs (3/4 aa long).

I think MEME is a good tool for your purpose, but there are others as well. I can think of InterProScan for example (although I am not sure if it is really de novo). Here a summary of available tools.


text:
Combining several different affymetrix expression microarray experiments for further MDS analysis
Assume that I have three different Affymetrix microarray experiments, 2x U133 (from 2015 and 2017) and 1xU133 Plus 2.0 (2019). Each array contains several different samples, some of them are common for 2 or 3 arrays. I'd like to combine these thee results to perform further MDS analysis using whole sets.
I'm aware about some batch effect that I'd like to avoid.
When I do MDS using separately normalized data from three experiments, I can see that expression pattern of sample present in two of 3 matrices is the same but the normalized signal values are raised by a certain value. It visibly affects further MDS analysis.
My questions are:

Is there some bioconductor package or strategy to do such analysis avoiding bias and batch effect?
Can I just normalize all 3 experiments together using RMA? And then MDS analysis will be unbiased.

Thanks in advance!


text:
How to change sequence format in my alignment file?
I have fasta file with alingned several sequences (from MUSCLE), when I open it (e.g. notepad++) they look like:

And I want dot format of identities like  and then save it in txt file.
I failed using MEGA, BioEdit, Jalview and CodonCode, how can I do this?

Hi ... this is definitely in MEGA. Could you check again? Try updating to MEGA X?

@Michael yes I tried. Export to xls and xlsx - crashed. Csv - fine, but only first 1000nt. Unfortunately I have about 40 sequences ~25kb length each. MEGA can't handle this amount of data.

okay, it must be done because it is very easy to code e.g. via AlignIO in Biopython ... you'll need to adapt it, but all you are doing is splicing through an alignment.


text:
Problem with FTP server regarding data upload via Webin-CLI
I'd like to upload my data (only two sequences) into my project using Webin-CLI. I don't know why i so complicated process...
I read deeply into tutorial and what I did?

I created account and study (Project).
I downloaded latest release of Webin-CLI.
I downloaded and filled spreadsheet template for annotated sequences (https://www.ebi.ac.uk/ena/submit/webin).
I prepared the command for the files validation:

java -jar webin-cli-8.0.0.jar  -username my.email@my.email -password ********* -context sequence -manifest manifest.txt -validate
The validation was succesfull. So I decided to submit data:
java -jar webin-cli-8.0.0.jar  -username my.email@my.email -password ********* -context sequence -manifest manifest.txt -submit
And I recieved this error:
INFO : Your application version is 8.0.0
INFO : Please upgrade to Webin-CLI version 7.0.1 or later if you see the following error: Failed to initialise validator. Could not retrieve BioSample.
INFO : Connecting to FTP server : webin2.ebi.ac.uk
INFO : Creating report file: /my/path/./webin-cli.report
WARN : Retrying connecting to FTP server.
WARN : Retrying connecting to FTP server.
WARN : Retrying connecting to FTP server.
WARN : Retrying connecting to FTP server.
WARN : Retrying connecting to FTP server.
WARN : Retrying connecting to FTP server.
ERROR: Failed to connect to FTP server. Failed to upload files to server because of a system error.
ERROR: Some or all submissions failed. Please see application logs.

I read that it's good idea to use Aspera. So I installed it and I added -ascp to the command line and this is the result:
java -jar webin-cli-8.0.0.jar -ascp -username my.email@my.email -password ********* -context sequence -manifest manifest.txt -submit
INFO : Your application version is 8.0.0
INFO : Please upgrade to Webin-CLI version 7.0.1 or later if you see the following error: Failed to initialise validator. Could not retrieve BioSample.
INFO : Uploading files to : webin.ebi.ac.uk
INFO : Creating report file: /my/path/./webin-cli.report
WARN : Aspera upload failed. Client exit code : 1
WARN : Client error :

I checked the connection with default Aspera server (ping) amd it worked. I do not have any firewall. I'm struggling with this problem for long time. Of course I wrote to ENA support however I haven't recieved the answer so far... Maybe you know how to fix it?
Session Stop  (Error: Client unable to connect to server (check UDP port and firewall))

WARN : Retrying file upload.
WARN : Aspera upload failed. Client exit code : 1
WARN : Client error :


text:
Affymetrix tags (same ID's) present in different places of the genome
could you explain me how to interpret such situation?
Same tag aligned in several distint places, near different genes...
Is the annotation able to distuinguish such 'same' tags, during DGE analysis?
Why tags arent unique?

Is the annotation able to distinguish such 'same' tags, during DGE analysis?

No, you'll end up discarding such probes (assuming you have a reason to actually use microarrays still).

Why aren't tags unique?

Platforms such as HG-U133 date back to 2002 (meaning they were designed a few years before that), when the human reference genome (and the reference transcriptome) were not nearly as "finished" as they are today. Probes that may have made sense given what we knew almost 20 years ago may no longer makes sense now given what we know. Further, microarrays are rather conservative, in that they tend to primarily gain probes over time and not have them removed (lest companies lose customers who are relying on given probes).

Thanks, for the answer! Ok, so I have the normalized expression signal from this tag. How can I distuinguish to what gene is this signal corresponds to? For example (very abstract) I have some ambiguous tag which is located in three different genes. So the signal is cumulative value of thesee 3? As a result I am not able to define the expression of these genes separately?

As I said, you'll end up excluding this probe most of the time. In modern arrays genes will have multiple probes associated with them (sometimes multiple per transcript). If you really wanted to incorporate this probe, you'd distribute it in some sort of weighted fashion to all genes/transcripts for which it's informative (in a similar to what salmon/kallisto/rsem do for RNA-seq reads). Depending on how elaborate you wanted to get that could be a lot of work, so typically you'd just exclude such probes instead.

I always assumed that these tags, such as 239963_at you show, mean that the sequence is not specific of a genome region. This means it cannot be assigned to any gene for DGE.
It could be due to changes on the reference genome or due to a faulty design of the tag sequence. However once in the microarray template they cannot be changed without introducing more problems, that's why doing RNA-seq prevent this.
Sorry I don't have any reference for this (but if someone wants to add this edit the answer).


text:
GWAS, MWAS, EWAS: what are the (in)dependent variables?
I started reading some papers on X-wide association studies, where X can be metabolome, epigenome, etc... The authors usually describe which are the dependent and which are the independent variables of their models. Taking as an example a metabolome-wide association study (MWAS), the authors use as model the following: metabolite ~ condition. That is, each metabolite is considered the dependent variable and the condition (e.g., smoker vs. non-smoker) is considered the independent variable. I was a little bit confused, so I looked into an example study using GWAS: in this case the model can be written as condition ~ SNPs, which is what I would expect (knowing nothing about GWAS but something about Data Science) and quite the opposite from the MWAS described above.
Does the choice of the dependent variable depend on the questions asked in the study, or does it exist some sort of agreement? Moreover, regarding the GWAS model described above, I read that usually researchers fit a model for each SNP: what's wrong with fitting a model which consists of several independent variables (i.e., y ~ SNP_1 + SNP_2 + ... + SNP_N)?

Fitting a model that consists of many independent variables is done, but serves a different purpose. There is nothing inherently wrong with it. A GWAS answers the question "given phenotype and genotype data, which genotypes are associated with the phenotype?" Your proposal answers a question more like "given phenotype and genotype data, how can we best predict phenotype with genotypes?"

See also: https://biology.stackexchange.com/questions/111538/why-would-greml-in-gcta-treat-snp-effects-as-random-effects/111539#111539

To answer the first part of your question, the dependent and independent variable of X-WAS is kind of arbitrary and dependent on the question you asked. But it gradually becomes a convention in the field after the initial name and concept are accepted by the community. For example, GWAS from the very beginning is written as condition ~ SNP, and there is no doubt about it. Later, people did something like condition ~ gene expression and called it TWAS, which is also widely accepted now. In your case, MWAS is indeed metabolite ~ condition, and the only important thing is to know the convention in the research community so people are on the same page when they talk to each other.
As for why running GWAS using single SNP. In a GWAS, it is typical to use 500K to 7M SNPs. So this is a large number. It is obvious that we cannot include all the SNP in one single model, then the question becomes how many SNPs to include in each model, and how to choose those SNPs. This is very hard to decide. Also, due to link disequilibrium, many SNPs are correlated, which means co-linearity, a thing you want to avoid in linear models. It is not clear how to correctly interpret the effect sizes in a model with many SNPs as well. Therefore, due to feasibility, simplicity, and interpretability, it might be best to perform regression on every single SNP.
Now since we are running GWAS in an SNP-by-SNP fashion, there are a lot of work that needs to be done after GWAS in order to find the true causal SNP. Many fine-mapping methods have been developed to tackle this problem.
Interestingly, there are some multi-variate GWASs methods, but they do not come in the flavor of simple linear models. Instead, sometimes complex Bayesian statistics are used, for instance.


text:
Methylation data: beta values are normally distributed
I am trying to incorporate Methylation data (Illumina 450K) into my pipeline. I was provided with two versions of the same data: one (dataset $a$) normalized with ComBat and the other (dataset $b$) normalized using the residuals method (i.e., taking the residuals of the models fitted on the covariates). Looking at the distribution of the beta values, I noticed that:

The values for dataset $a$ have a typical bimodal distribution with peaks at around $0$ and $1$
The values for dataset $b$ are roughly normally distributed with the mean at around $0.5$.

I am not that familiar with Methylation data, so maybe this is expected. Is it? I did not find any paper showing similar results.

$\beta$ values are the ratio of methylated to unmethylated probe intensities.
This means that $\beta \in [0, 1]$, so it can't be normally distributed, by definition.
Its distribution may have sections that look like bell curves, but it isn't normally distributed.
You can take a look at this paper by Du et al. [1] for details about the statistics of $\beta$ values.
That said, the distribution of all $\beta$ values from an array or bisulfite sequencing (or some other DNA methylation assay) should be bimodal with modes at 0 and 1.
Nearly all CpGs in the human genome are methylated or unmethylated.
There may be some hemi-methylated sites, but for the 450K array you mentioned, there should be few of these.
It sounds like you might be looking at some residual value, or some other value from the analysis in the paper for dataset $b$ other than $\beta$.
Could you be looking at $M$ values, or something else?
[1]: Du et al., BMC Bioinformatics, 2010

It's not clear from your question what the statistical design of the experiment is, so I can't judge whether or not a distribution around 0.5 would be expected from the residuals after fitting.
Given that dataset a looks like what you expect, I'd be inclined to use that (and perhaps ask the person who gave you this data why dataset b looks so different, and why it ended up with a 0.5 mean).


text:
Get genomic coordinates using GenomicFeatures by HGNC gene names
I want to get coordinates of human genes from my list (consisting of hgnc genes id) using GenomicFeatures and TxDb.Hsapiens.UCSC.hg19.knownGene R packages from Bioconductor.
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb=(TxDb.Hsapiens.UCSC.hg19.knownGene)

my_genes = c("INO80","NASP","INO80D","SMARCA1")

select(txdb, keys = my_genes,
   columns=c("TXCHROM","TXSTART","TXEND","TXSTRAND"),
   keytype="GENEID")

However, it doesn't' work because txdb doesn't take hgnc identifiers; how can it be solved? I couldn't find any appropriate keytype that will support hgnc and not sure how to match hgnc id I have and GENEID from txdb.

Answer from @devon-ryan, converted from comment:
Try a UCSC gene ID instead. For those genes they'd be: uc001zni.3 (INO80), uc002vaz.4 (INO80D), uc004eup.4 or uc004eun.4 (SMARCA1) and uc001coj.2, uc010olq.2, uc010olr.2, uc001coi.2 or uc001col.2 (NASP). You might have to omit the versions (the .4 part, for example).


text:
How to map selected genes to Metabolic pathway Maps
I have a selected Arabidopsis Genome Initiative (AGI) list for RNA seq and proteomics data, how can I map them to metabolic pathway maps to vilualize(e.g. TCA cycle / FA / Photosynthesis) in KEGG or AraCyc or any freely available pathway maps?
Could someone share an example if possible?

KEGG has a webtool for this, KEGG Mapper.
You can just copy-paste your gene ids in there and see how your pathways are colored.

Thank you Benn :)


text:
How to visualise /add Arabidopsis T-DNA Lines As Tracks On IGV
We have done some RNA-Seq on Arabidopsis T-DNA lines and mapped the reads on to Arabidopsis genome. Now we wanted to confirm the T-DNA deletion of those lines by loading them on IGV. As a control we also wanted to load T-DNA lines as tracks on IGV to make sure that the T-DNA lines we chose were correct. I am wondering where can i get bam files for those confirmed T-DNA lines?

Sorry I might be missing something. You mentioned mapping the reads onto the genome. The output is a bam file right? You can just create an index for it and load it into IGV

Are you asking how to see the location and names of all the T-DNA lines, like shown here?


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

Hello Akhil, this is not the greatest place for open questions. You could try one of the forums or just google it.

As others have said, there is no short cut and you will need to spend some time with a good book.
I would say there are 3 broad areas of understanding you have to have as a structural bioinformatician:

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

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

Thank you. If you could also suggest me some resources(your top picks) for the first two areas, it'll be really helpful.

Of the top of my head "Introduction to Protein Structure" by Branden and Tooze is good. In terms of databases have a look at PDB, UniProt, Pfam, CATH and InterPro.

Thanks @jgreener

A good starting point could be this book by Gu and Bourne:
Structural Bioinformatics, 2nd Edition

Thanks, I'll check it out. But I was hoping for a good blog or some tutorial(web/video) kind of resource rather than a book

@AkhilVerma — Books are far more informative. They contain more information, undergo quality control if published by reputable publishers and are generally better written, more comprehensive and more reliable. You cannot learn anything worthwhile without putting in some effort. Certainly not structural biology. How’s your chemistry, incidentally? You’re going to need that too


text:
Normalize by counting CNV within a pool of VCF files
I’m (trying) using the GATK4 germline CNV calling pipeline. I successfully got 57 VCFs from my sample batch, called with segments (obtained by merging the contiguous intervals), like in a classic VCF :
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  2046745451-1006_S4
M       3288    CNV_M_3288_15907        N       <DEL>,<DUP>     .       .       END=15907       GT:CN:NP:QA:QS:QSE:QSS  2:5:9:17:6:21:21
1       69071   CNV_1_69071_70028       N       <DEL>,<DUP>     .       .       END=70028       GT:CN:NP:QA:QS:QSE:QSS  1:0:1:204:204:204:204

But I got way too much of those intervals, more than 10k. I would like to know if there is an existing tool which count the different segments (variants / intervals common by +/= 75% of their length) in one VCF and gives me the count of the different segments overlapped by segments in other sample in my batch. By counting the most redundant segment, I could determine which are background noise and maybe decrease the number of variants in my VCF by filtering.
Thank you.


text:
Get results of keyword search on Pfam via python script
I'm interested in all proteins that are in any way associated with Danio rerio. I decided to look them up at Pfam data base and when I just make a keyword search, I get a a nice list which looks like this http://pfam.xfam.org/search/keyword?query=danio+rerio&submit=Submit.
Important thing to note is that most of the associations of these protein families can be found in seq_info column.
My question is: how to extract that list from a python script? Moreover, I want to make a script that when run, will fetch the latest release, filter out the proteins that cannot be found in D. rerio and do some additional processing.
I've been searching for ways to directly download files from Pfam ftp server and mine for the proteins that are associated with D. rerio, however, I didn't find a single way to do it. I simply cannot see where do they store the information about species in which the protein can be found.
I would like to retrieve everything about protein families of interest, because the next step I will need to do is to find the DNA binding proteins and for finding a way to filter out the ones of no interest, I would liketo have as much data as possible.

If I understand your question correctly, then you can use the Pfam proteom search to get the data you are looking for.
Go to the

Pfam website -> Browse -> select the letter D from proteoms category
  -> D. rerio

There you can find a download link under the section Domain composition:

You can download a list of all regions for this proteome from our FTP
  site.

The download link seems to be bound to the species ID and if you want to automate your data retrieval you may be able to simply download the file from there in your python script.

Exactly what I needed. Although, the species ID is not that informative but I think I can manage that easier in the future. Thank you!

One thing to note though is that seq_info at Pfam is retrieved from UniProt database so in order to retrieve that data, one needs to make a mapping from list of proteoms and uniprot.gz file at the Pfam ftp server. The tricky part comes if the two databases are not synced. But that's not an issue that should be discussed here.


text:
A basic question about phylogenetic inference with Bayes-MCMC
I am new to the world world of Bayesian-MCMC inference and have been working my way through various texts. I am having some difficulty understanding how MCMC works. I understand that it is a way of sampling from the posterior distribution. What is puzzling me, however, is how the probabilities of the posterior distribution have been calculated. That is, if the Markov Chain is sampling from the posterior distribution, that posterior distribution must have probability values on some axis. How were those calculated?

In short, the posterior is a probability value between 0 and 1 which is calculated from Bayes rule, which In phylogenetic usage is given as the product of the prior probability and the likelihood divided by a normalising constant. The posterior probability allows you to state the probability of some hypothesis given the data at hand (as opposed to the likelihood parlance, which is probability of data given hypothesis).
In a phylogenetic setting, the sampler works (roughly) like this: propose some value of a parameter (topology, branch lengths, etc.) and use the specified prior and calculated likelihood for that parameter value given the data to calculate the unnormalised posterior. Next, propose a new value, and repeat the same step. Compare the two values, and if the posterior probability of the second is greater, accept this parameter value and repeat the cycle again. In this way the Markov Chain moves through posterior space until it reaches a reasonably stable optimal approximation of the posterior distribution. Here, we want to get as many samples as possible to maximise the accuracy of our posterior. The maximum a posteriori (MAP) value of the parameter and its credible intervals, etc, will come from the distribution of these samples.

Thank you so much! This is very clear and helpful. There is just one detail that I am hung up on, namely this step: "use the specified prior and calculated likelihood for that parameter value given the data to calculate the posterior." How can this be calculated without knowing the marginal probability/normalizing constant?

@Namenlos this is actually related to a property of the Metropolis-Hastings algorithm: because each step of the chain explores ratios of a value which is proportional to the posterior (that is, the likelihood multiplied by the prior; the numerator of Bayes' rule), the denominator of the posterior cancels out. Samples from the chain therefore give you an estimate of the posterior probability. I've edited my answer to reflect this.

Extra details can be found in this excellent introduction to Bayesian phylogenetics

Ah, I see. Thank you again. This is starting to make sense now.


text:
Why does `pratchet()` in `phangorn` not find the most parsimonious trees?
I called pratchet() on a dataset with binary values for the characters and it returned a parsimony value of 35. I then calculated the parsimony value for the dataset with a tree that I myself constructed and it was 7. How is this possible at all? I thought that pratchet() was supposed to return the most parsimonious tree on the basis of the data set, but that clearly did not happen in this case. Is there other software that is guaranteed to find the most parsimonious tree?

Please provide a minimal working example of code you have tried. That way a concrete and useful answer can be provided.

Perhaps the tree which you have calculated is not being caught in the search heuristic? Have you tried (as a trial) setting the tree which you calculated as the starting topology, just to check that it is indeed being calculated on by pratchet()?
Alternatively, you could try playing around with the search heuristic settings, like varying from ratchet to NNI, or you could also try using an external program, like PAUP*.

Thank you, I tried these. When I set the start tree to the tree that I created manually, I got the following result:                                                                    Error in edgeMatrix[, k + 1] <- c(p1, e12, e34, p2, 1L) : number of items to replace is not a multiple of replacement length                                                         What does this mean?

I am unsure. Have you verified that you provided a bifurcating newick tree file and that the data matrix has the same number of taxa as the tree?

The number of taxa is the same, but the tree is not fully resolved. So I presume that is then the problem.

For dendrograms with less than about 10 OTUs/taxa/leaves it is possible to evaluate all possible dendrograms using an exhaustive search algorithm. However, for more than about 10 OTUs, heuristics need to be employed to try and "home-in" on and score a smaller set of "better" trees. However, these heuristics are not guaranteed to find the best dendrogram.
Without some specific code it is hard to do more than speculate as to the problem you are encountering, especially where you get an error using data you manually create.


text:
Do any computational phylogenetic methods enable the specification of ancestral states?
Various phylogenetic algorithms estimate ancestral states of a phylogenetic dataset. Is there a way in either maximum parsimony, distance-based methods, or Bayesian inference to indicate what the ancestral states of characters were?

In what respect? Do you mean as in looking at the evolution of a character on a tree? If so, there are algorithms for constrained Brownian motion. If you mean in terms of tree estimation, then knowing the ancestral states implies that you know something about the relationships already - if so, you can use constraint trees for some relationships in ML and Bayesian inference.

Otherwise if you are just looking for ancestral character state estimation, look at the ace() function in the ape package, or make.simmap() in phytools.

Thank you. I had in mind specification of ancestral states in terms of tree estimation, so I will investigate constraint trees in ML and Bayesian inference.

For maximum parsimony, Fitch's method seems simple and intuitive, and can be done in the R phangorn package: Fitch, Walter M. (1971). "Toward Defining the Course of Evolution: Minimum Change for a Specific Tree Topology". Systematic Zoology. 20 (4): 406. doi:10.2307/2412116. ISSN 0039-7989. JSTOR 2412116
There are various methodological criticisms of parsimony that you can find in the literature, but its appeal is you can easily understand why the algorithm came to the conclusions it did.
RAxML can do maximum likelihood reconstructions based on a tree and tip values, and I think Mr. Bayes can do a Bayesian version.

Curious, it seems like Systematic Biology has the same paper here

The 70s sure were a crazy time!

Yes this approach is done in parsimony, i.e. you can indicate the outgroup and the ancestral state is calculated from there. I assume this is also performed using 'Beast' because it is a non-reversible matrix. Regular ML and Bayes approaches use reversible mutation matrices therefore the tree is calculated first, you can then specify an outgroup to calculate the ancestral state. However, in a reversible matrix the ancestral state is not used in the calcuation. I hope that makes sense.


text:
Using cellranger mkgtf w/ NCBI GFF3 files
I am using several cellranger tools/pipelines to analyze some data, specifically, mkgtf, mkref, count, and aggr. I am using mkgtf and mkref to include a PI-synthesized protein to my reference. I have finished the analysis but am not seeing the protein in my output, so I am going back to the start to verify the process.
Unlike the documentation (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references), my GTF file was not downloaded from ENSEMBL or UCSC, but NCBI (https://www.ncbi.nlm.nih.gov/nuccore/AY678269.1?report=genbank). This GFF file does not have the same key nomenclature as seen in ENSEMBL or USCS.
I did use the mkgtf tool, but the created file is identical to the input.
cellranger mkgtf custom_ref/AY678269_genes.gtf reference_data/AY678269_mkgtf_genes.gtf

I went ahead and used cat to add the fasta and gtf file from the synthesized protein to the respective reference fasta and gene files, then ran mkref.
The s.protein is in the fasta file:
>lcl|AY678269.1_cds_AAV52169.1_1 [protein=tandem-dimer red fluorescent protein] [protein_id=AAV52169.1] [location=1..1431] [gbkey=CDS]

as well as the gene list
##sequence-region AY678269.1 1 1431
##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=32630
AY678269.1      Genbank region  1       1431    .       +       .
AY678269.1      Genbank CDS     1       1431    .       +       0

However this protein is not found in the output...surely I must have made an error. Does anyone have any ideas about what I could be doing incorrectly?
Do I need to change the entire nomenclature of the GFF keys to match what I see from ENSEMBL (e.g., gene_id ""; gene_version ""; transcript_id, etc...)?
Any assistance is greatly appreciated.

Regarding your comment:

This GFF file does not have the same key nomenclature as seen in
  ENSEMBL or USCS

The link you have provided corresponds to an entry in Genbank format and not to GFF format. I believe you need a file in GTF format in particular.
Nevertheless, given the fact that you have manually added the corresponding entries to the reference fasta and gtf files prior to mkref, you should have seen your gene of interest in the output. I would double check one of the original fasta headers and original gtf entries in the reference file and match these for the new entry.

Thank you for the response. I figured out what the problem was, I was given an incorrect accession number! However, after finding & reading through this: https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/tutorial_mr, I can see how my approach may have been problematic anyway because of the NCBI keys. Fixed it & I'm creating the custom reference. FYI, the link does go to a genbank file, but when downloading the information, you have the option of choosing the format (fasta, genbank, gff), but you are absolutely correct about the GTF vs GFF. Danke!


text:
Is there a way to calculate the number of parsimony informative characters in a morphological dataset in R?
I have a dataset of binary morphological characters and would like to know how many of them are parsimony informative. Is there a way to calculate this with an R package?

The ape library (Analyses of Phylogenetics and Evolution) using the MPR method (Most Parsimonious Reconstruction).
Information on ape is here.
The format
install.packages(“ape”)
library(ape)
MPR(x, phy, outgroup)

x :   a vector of integers.

phy : an object of class "phylo"; the tree must be unrooted and fully dichotomous.

outgroup : an integer or a character string giving the tip of phy used as outgroup.

Rooting the tree is recommended. In the example rooting isn't used but an outgroup is used as taxa t1, that will be fine.
Generating random data,
vectx <- rpois(50, 1) # make random data
tree <- rtree(50, rooted = FALSE)  # also rtree(50) make random tree
MPR(vectx, tree, "t1")

Thats it, there's nothing more. Data, tree = steps.
Clarification the question is number of parsimonious characters. As a strict definition the number of informative characters is the same as the number of steps and thats what MPR will describe ... thats easy.
There are other possible meanings

The number of parsimonious characters/steps to make a tree is calculated via MPR as above.
if its the actual changes at a given position in the alignment (or vector) for each node in the tree - thats the answer below - but thats complicated.
If its simply is this site (not character) is parsimoniously informative and this one isn't then the "alignment" or vector could be used one character position at a time and subject to MPR and the step length would determine how parsimoniously informative each site was ... thus point 3 is easily doable in R

Note Calculating the actual sites on the alignment and associated changes per node is easily done outside R, but within R I don't know if thats doable.
Wait, I'm wrong it's doable but I would envisage its difficult to implement. Calculating parsimonious characters for each given node in a tree - yes I do that - doing that in R ... errhhh thats not me. There are Pythonic ways in, but I've not implemented them.


text:
How is the instantaneous rate of an instantaneous rate matrix determined?
In molecular evolution, I know that the instantaneous rate of the rate matrix is the limit of the rate as time approaches 0, but what I don't understand is how this rate is established or what I can do with this number.
For instance, if a discrete character is thought to have changed from 0 to 1 five times on a tree with a root age of 40 million years, is the proportion 5 / 40 million the same as the instantaneous rate of change? If not, how exactly are they different?
From a different perspective, let's assume an (absolute) instantaneous rate of 0.0001 and let's say I have a tree with a root age of 2 million years. Can I multiply 0.0001 by two million to get the expected number of transitions of that character across the tree?

related: https://bioinformatics.stackexchange.com/questions/13605/what-does-instantaneous-rate-matrix-mean

@MaximilianPress I saw this post. Although both questions concern rate matrices, they are nevertheless quite distinct.

I agree; note that I didn't mark it as duplicate. However, they have a substantial overlap and it is worth annotating as such, for anyone coming across this question.


text:
add task to nextflow and resume
in test.nf (dsl2) I have
   ...
   workflow{
       a = Channel...
       processA(a)
       processB(processA.out.emitA)
       processC(processB.out.emitB)
    }

The pipeline completed successfully.
Now I want to add processC_2 that takes the output of processB.
How do I run the workflow to bypass processA execution?
I added processC_2 to test.nf and then did:
nextflow run test.nf -resume run-name

But it doesn't run the newly add process.
How do I achieve my goal?
Thank you!

It is actually very hard to answer this using only the information you have provided. It'd be necessary to actually see your test.nf file and read the processes. The execution seems right and I'd guess that the processB may be not emiting any output as shown in https://www.nextflow.io/docs/latest/process.html#outputs.
It is also possible that when using the -resume together with a specific run-name you may be recovering the state of that specific run as shown here https://www.nextflow.io/blog/2019/demystifying-nextflow-resume.html. You can try to execute it without the run-name or without the -resume flag so you can check whether the problem is with the emition and channelling of outputs between processes, or with the -resume flag.


text:
Use conditional in workflow in Nextflow dsl2
How should I do:
...
if ($(params.aligner) == "bowtie") {
    align_bowtie2(get_fq_info.out.fq_info)
}
else {
    align_bwa(get_fq_info.out.fq_info)
}
...

You can move that into your workflow definition. There's an example in the docs of how best to use workflow parameters:

A workflow component can access any variable and parameter defined in
the outer scope:
params.data = '/some/data/file'

workflow my_pipeline {
    if( params.data )
        bar(params.data)
    else
        bar(foo()) }

Your example workflow would look like:
workflow {

    fq_info = get_fq_info.out.fq_info

    if( params.aligner == "bowtie" ) {
        align_bowtie2( fq_info )
    }
    else {
        align_bwa( fq_info )
    }

    ...
}
```

Is this inside the main.nf file or inside the process.nf file? Anyway, the param.anything variable should already be available in the nextflow context when executed, thus, you could simply do:
if (params.aligner == "bowtie") {
    align_bowtie2(get_fq_info.out.fq_info)
} else {
    align_bwa(get_fq_info.out.fq_info)
}

Check out this main.nf of a nextflow pipeline of mine, it contains some working if/else snippets.

Here it is another StackOverflow question on conditional nextflow pipelines that might help.

Hope it helps.


text:
Can bwameth accept unaligned Bam?
Does bwameth accept bam (unaligned) type as input or I will have to convert unaligned bam to fastq (samtools fastq bam) to use it as input?

No, it's a bwa mem wrapper for bisulfite data essentially and bwa mem accepts fastq only.

You can pipe uBAM to bwamem with samtools fastq in.bam | bwa mem -p - > out.sam, provided that paired-end reads are grouped by name.


text:
samtools mpileup skipping read
I run the command:
samtools mpileup -O -s -q20-B -Q20 -f hg19.fa -r chr1:569929-569931 myFile.bam

and get:
chr1    569929  G   7   ...,,., EEEEEEE NTTVTTU 53,48,42,60,30,29,27
chr1    569930  C   6   ...,.,  EEAEEE  NTTTTU  54,49,43,31,30,28
chr1    569931  G   6   ...,.,  EEEEEE  NTTTTU  55,50,44,32,31,29

I am interested in position 569930.
I inspect the bam with:
samtools view myFile.bam chr1:569929-569931

and I see:
NS500355:NS500355:HHYN5BGXB:2:21311:2303:9255   99  chr1    569869  0   76M =   569900  107 CTCCATAACGCTCCTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGA    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE<EEEEEEEAEEEEEEE    NM:i:2  MD:Z:5C3C66 MC:Z:76M    AS:i:66 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9321,76M,1;
NS500355:NS500355:HHYN5BGXB:4:22603:25857:2799  99  chr1    569877  45  76M =   569896  95  CGCTCCTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCAC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE    NM:i:1  MD:Z:1C74   MC:Z:76M    AS:i:74 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9329,76M,1;
NS500355:NS500355:HHYN5BGXB:1:12105:11868:12011 99  chr1    569882  51  76M =   569900  94  CTCATACTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACC    /AAAAEEEEEAEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9334,76M,1;
NS500355:NS500355:HHYN5BGXB:2:13301:2681:6114   99  chr1    569888  9   76M =   569952  140 CTAGGCATACTAACCAACACACTAACAATATACCAATGATGGAGCGATGTAACACGAGAAAGCACATACCAACGCC    AAAA/E/EE//AAEEEEE/EAAEAA//EEAEEE/E/EEEE/////A/EE/EE/EE/E//EE/EEEE/E/AEA####    NM:i:4  MD:Z:6C19C15C29G3   MC:Z:76M    AS:i:57 XS:i:52 RG:Z:../H2O XA:Z:chrM,+9340,76M,5;
NS500355:NS500355:HHYN5BGXB:2:23110:10125:9334  99  chr1    569888  51  76M =   569904  92  CTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCC    AAA/AEE6EEEE6E///E/A/AAE/E/AEEEEEEE/EE/EEEAEEEEAEE//EEEEEEAEEEAEEAE/AEE/AEE<    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9340,76M,1;
NS500355:NS500355:HHYN5BGXB:1:13209:10692:11045 83  chr1    569888  53  18S58M  =   569888  -58 CGTGTGCTCTTCCGATCTCTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAG    EE<A/EE//AE/AA///AAE///EEEEE//6///E//EE/6EE/EEEEEEEA/A/A/E/E//AEAEEEEE6A6AAA    NM:i:0  MD:Z:58 MC:Z:59M17S AS:i:58 XS:i:53 RG:Z:../H2O XA:Z:chrM,-9340,18S58M,1;
NS500355:NS500355:HHYN5BGXB:1:21304:26480:15673 83  chr1    569890  11  76M =   569796  -170    AGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCAC    AEEEEEEEAAEEEEEEEEEEEE/EEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,-9342,76M,1;
NS500355:NS500355:HHYN5BGXB:4:11602:21237:3028  83  chr1    569900  51  76M =   569900  -76 ACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCAC    EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,-9352,76M,1;
NS500355:NS500355:HHYN5BGXB:3:13504:4369:8812   99  chr1    569901  51  76M =   569901  76  CCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9353,76M,1;
NS500355:NS500355:HHYN5BGXB:2:11203:19774:18704 99  chr1    569904  11  76M =   569952  124 ACACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGT    6AAAAEEEEEEEEEEEEEAEEEEEEEEEEEEEEEE6EEEEAEEEAEEEEEEEEEEEEEEEEEEAE/E6EEAEEEEA    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9356,76M,1;
NS500355:NS500355:HHYN5BGXB:2:11309:13954:17567 99  chr1    569905  11  76M =   569972  143 CACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTC    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9357,76M,1;
NS500355:NS500355:HHYN5BGXB:1:21202:13252:14049 99  chr1    569905  11  76M =   569927  98  CACACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTC    AAAAAEE6EEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEEEE/EEEE<EEEEEEEEEEE    NM:i:0  MD:Z:76 MC:Z:76M    AS:i:76 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9357,76M,1;
NS500355:NS500355:HHYN5BGXB:1:22311:21964:12072 83  chr1    569908  52  5S71M   =   569908  -71 GATCTACTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTG    <AEEEA6EEE6EEEEEEEEEEEEEE/EEEEEEEEEEEEEEEEE/EEEEEEEEEEEEEEEEEEAAEEEAEEEAA6AA    NM:i:0  MD:Z:71 MC:Z:71M5S  AS:i:71 XS:i:66 RG:Z:../H2O XA:Z:chrM,-9360,5S71M,1;
NS500355:NS500355:HHYN5BGXB:4:12602:7858:19418  99  chr1    569909  15  76M =   569917  84  CTAACCATATACCAATGATGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTCCAAA    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEE/AAEEEEEEEEEEEEEEE    NM:i:1  MD:Z:74G1   MC:Z:76M    AS:i:74 XS:i:71 RG:Z:../H2O XA:Z:chrM,+9361,76M,1;
NS500355:NS500355:HHYN5BGXB:2:12205:19456:7598  83  chr1    569927  0   76M =   569907  -96 TGGCGCGATGTAACACGAGAAAGCACATACCAAGGCCACCACACACCACCTGTCCAAAAAGGCCTTCGATACGGGA    EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAA/EEEEEEEEEAAAAA    NM:i:1  MD:Z:56G19  MC:Z:76M    AS:i:71 XS:i:76 RG:Z:../H2O XA:Z:chrM,-9379,76M,0;

Based on mapq and positions (respect to 569930) I see that the read with cigar=18S58M is skipped by mpileup.
WHY? The soft clip is way far from position chr1:569930

Examining the mpileup output for a single read on its own is a good way to figure out why bases are not appearing as expected.
So if you suspect that the cigar=18S58M read is the missing one (you are correct), prepare a SAM file containing just that read:
@SQ SN:chr1 LN:248956422
NS500355:NS500355:HHYN5BGXB:1:13209:10692:11045 83  chr1    569888  53  18S58M  =   569888  -58 CGTGTGCTCTTCCGATCTCTAGGCCTACTAACCAACACACTAACCATATACCAATGATGGCGCGATGTAACACGAG    EE<A/EE//AE/AA///AAE///EEEEE//6///E//EE/6EE/EEEEEEEA/A/A/E/E//AEAEEEEE6A6AAA    NM:i:0  MD:Z:58 MC:Z:59M17S AS:i:58 XS:i:53RG:Z:../H2O  XA:Z:chrM,-9340,18S58M,1;

and run it through mpileup:
$ samtools mpileup -Q0 -B oneread.sam

This will show you that this read's base quality dips below 20 at the position you are interested in. Thus, as finswimmer suspected, your -Q20 is removing it at this position.

Why do you believe it's the read with cigar=18S58M? The mapping qualities given in sam files are different to those in the mpileup output. I guess samtools is doing some recalculation.
You are using the parameters -q20 -Q20. So you skip all reads with mapping quality <20. This will result in 7 reads got to mpileup. Furthermore you are skipping bases with base quality <20. Have a look at the base qualities for the read starting at position 569904. If I count it correct you end up with a quality of / at position 569930 which means 14. So I guess this read is dropped.
BTW: samtools mpileup is deprecated (only when creating bcf/vcf)  and should be replaced by bcftools mpileup.

Thank you, I missed the '/' quality for a 'E'. Now I understand.

The mpileup subcommand has been split between the two commands. Textual “mpileup output” stays in samtools (and always will), but VCF/BCF output has moved to bcftools. This groups the latter alongside the other VCF-related facilities and in particular means that both halves of the calling pipeline bcftools mpileup … | bcftools call … are in the same program, which simplifies updates and improvements to this pair of commands that work together.


text:
nextflow `fromPath` not passing all files with *
main.nf is
#!/usr/bin/env nextflow

params.input = ''
tchan = Channel.fromFile(params.input)

tchan.println()

I have a.bam, b.bam and c.bam
I run:
nextflow run main.nf --input *bam

But it only prints a.bam
According to the documentation, wildcard character it is interpreted as a glob path matcher
What am I doing wrong?
Thank you

Shell expands wildcards, try ls * vs  ls "*" .

got it!
I should run:
nextflow run main.nf --input "*bam"

instead of
nextflow run main.nf --input *bam


text:
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.

Picard was running (exit status was 0 when I set to 1 if there was an error) The log was the same as when I run the command in the terminal. I just excluded the prefix "./" in naming the output and it worked (I guess galaxy add a path to the filename with is not "./" ?).
Thanks!


text:
Installing DESeq2 in Ubuntu
I am trying to install DESeq2 in my Ubuntu with R version 3.5.1. According to the package repository in Bioconductor the version should be 3.5.
> R.version

platform       x86_64-pc-linux-gnu
arch           x86_64
os             linux-gnu
system         x86_64, linux-gnu
status
major          3
minor          5.1
year           2018
month          07
day            02
svn rev        74947
language       R
version.string R version 3.5.1 (2018-07-02)
nickname       Feather Spray

then I attempt to install DESeq2, following the instructions in Bioconductor:
> if (!requireNamespace("BiocManager", quietly = TRUE))
       install.packages("BiocManager")
   BiocManager::install("DESeq2", version = "3.8"

and at the end I get:
Warning messages:
1: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘RCurl’ had non-zero exit status
2: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘GenomeInfoDb’ had non-zero exit status
3: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘annotate’ had non-zero exit status
4: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘GenomicRanges’ had non-zero exit status
5: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘genefilter’ had non-zero exit status
6: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘geneplotter’ had non-zero exit status
7: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘SummarizedExperiment’ had non-zero exit status
8: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘DESeq2’ had non-zero exit status

I have visited $2*10^{9999}$ discussions on this issue and none of them helped me solving the problem. In addition, I didn't have any problem installing other packages like limma and edgeR, so I assume that this might be specific to DESeq2 or one of it's dependencies.
Thank you for the help!

First I installed this libraries:
sudo apt-get install libcurl4-openssl-dev libxml2-dev

I then installed RCurl:
BiocManager::install("RCurl")

and lastly:
BiocManager::install("DESeq2", version = "3.8")

Have you tried just typing,
install.packages("BiocManager")

Without all the proceeding if statement? If this doesn't install correctly nothing else will.  Personally I would be tempted install each of the named libraries independently. If that works okay, e.g. its not a problem with the mirrors then I'd contact the developer michaelisaiahlove at gmail.com

BiocManager is successfully installed. The problem is with DESeq2...

Thanks it worked! I will post the answer for others with the same problem

Cool, happy to hear! Systems admin is the bane of every bioinformatician.

I noticed that some packages sometimes failed to compile depending on the version of gcc present on the system. So far, most of my problems were solved using a minimal R 3.5 installation, gcc 4.9 and trying to install everything using Bioconductor.


text:
Add new tool to galaxy
I am trying to include a new tool to the galaxy main menu, following these instructions.
tool_conf.xml doesn't exists only tool_conf.xml.main and tool_conf.xml.sample and none of them seem to be working. I make changes and nothing happens when I run galaxy.
What am I missing?
Thank you.

The .main and .sample file files usually come with default installations and may have been removed by whoever setup this Galaxy instance. You can look at the originals online and just base your own file on those: https://github.com/galaxyproject/galaxy/blob/dev/lib/galaxy/config/sample/tool_conf.xml.sample

Thanks @Devon, I removed the suffix .sample and it worked. I believe that galaxy folks could improve/update their manual.


text:
imgtrefseq vs imgtblast
In IMGT ftp-download (flat files) page, there are 2 optional fasta files (imgtblast and imgtrefseq) with their corresponding blast databases to download. Is the only difference the syntax and additional annotation in the refseq files?
Thank you!


text:
use same output in two processes in nextflow dsl2
This is my workflow:
pre_align()
pre_align.out.single_fastqs.view()
get_fq_info(pre_align.out.single_fastqs)

align_bwa(get_fq_info.out.fq_info)

align_bowtie2(get_fq_info.out.fq_info)

where I want to use the same output from get_fq_info as input for either align_bwa and align_bowtie2`.
As it's shown here only align_bwa is executed. If I comment align_bwa THEN align_bowtie2 gets executed.
How can I fix this?
Thank you!

I have solved this problem in the past by duplicating the output channel of the desired process in order to feed the same output into multiple downstream processes.

As far as I know, DSL2 should be able to handle this. Have you specified an output qualifier for your process get_fq_info (e.g. path)? I ran into a similar issue when I forgot to specify it.
For example, this will result in the error process align_bowtie2 declares 1 input channel but 0 were specified:
process get_fq_info {
   input:
   path reads

   output:
   "output.txt"

   shell:
   """
   myprocess.sh ${reads} > output.txt
   """
}

And is probably solved by this:
process get_fq_info {
   input:
   path reads

   output:
   path "output.txt"

   shell:
   """
   myprocess.sh ${reads} > output.txt
   """

```


text:
snakemake get environment variable for rule name and log path
When a rule fails in snakemake workflow, it returns a descriptive message in log file like following:
Error in rule step1:
    jobid: 99
    output: results/foo/step1/sample5/outdir/output1.tsv, results/foo/step1/sample5/outdir/output2.tsv
    log: logs/step1/sample5_out.log (check log file(s) for error message)
    shell:
    ...
        (exited with non-zero exit code)

Since I run workflow over many samples, I have defined a postrun.sh bash script in jobscript.sh directive that runs after snakemake {exec_job} section and notifies me via userbot in slack channel when a step failed for one or more samples.
jobscript.sh
        ########################### START SNAKEMAKE CMD ############################
    {exec_job}
        ############################ END SNAKEMAKE CMD #############################
    exitstat=$?
    export exitstat

    echo "END at $(date)"
. "${{HOME}}"/.smk_confs/postrun.sh

postrun.sh
## execute after running each of snakemake job
echo 'Running postrun.sh'
ERRMSG="MYJOB ID: $PBS_JOBID failed at $(pwd) on $(hostname) for ${USER} with exit status: ${exitstat}."

## ping to slack
ping_slack -m "$ERRMSG"
sleep 5

I am looking for a snakemake exported environment variable that I can use in postrun.sh to know which rule ID failed for one or more samples. For now, I can think of using job scheduler variable matching job name, e.g, PBS_JOBNAME, but better if I have a variable from snakemake API for rule ID or preferably, path to log file.
Thanks

You could add an exit trap to the beginning of the shell scripts inside your Snakefile, which is executed after each script is done and sends a message when there was an error:
 shell:
    """
    trap 'if [ $? -neq 0 ]; then ping_slack -m "Error in rule xxx with input {input}"; exit 1' EXIT
    ...

I believe this needs to be added to all of workflow rules. I am hoping to set a trap similar to existing onerror handler of snakemake but prefer to execute it at completion of each rule rather an entire workflow.

yes, that would be needed in each rule.


text:
High Phred Quality score VCF file
I am working on ION Torrent S5 XL and I am facing to an issue with VCF output generated by the variant caller of thermofisher. It stating on VCF 4.1 that Phred Quality Score is on a range of 0 to 255. In my VCF I have values over than 255. Do you know if there is a way to transform these values?

Let me rephrase to see if I understood: your output exceeds the max value explained by its manual? Could you post an example? My only guess is that either the structure is wrong or there is a bug, or the program is not well documented.

Hi Llopis, here's an example: chr17 7577104 . AGGA AAAA,AAGA,AGA,AGAA,AGGG,AGTA 35144.8 PASS. 35144.8 is QUAL Score into the VCF

Please can you [edit] the question to include it in the right formatting? Also could you post the command you used to get this VCF ?

See this post on the GATK forums.
It is normal to see values above 255.  The QUAL score is the posterior probability that all samples in your callset are homozygous reference. I quote Sheila's response from the thread posted above:

The QUAL score is the Phred-scaled posterior of AC = 0. We use the AC priors and the PLs to get the likelihood of the data given each AC, then use those to get the posterior probability for each AC. From there, the calculation is 1 - Pr{AC > 0}.

Please see the GATK post in Methods and Algorithms for the nitty gritty details. Beware, Bayesian statistics abound.

GATK's HaplotypeCaller and freebayes are two more variant caller that output QUAL values greater than 255. The reason for that is, that each read have its part on calculation that value. This results in the higher the read depth at this position the higher the QUAL value.
To get something useful out of this value, the usual thing that is done, is to calculate the quality per depth (QD) by just divide the QUAL value by the read depth.
Note whenever you try a different variant caller, you have to find out whether the QUAL value is read depth dependend or not.


text:
Negative scale-free topology
I am trying to build gene co-expression networks with WGCNA by combining in-house and publicly available RNA-seq datasets. I am interested in identifying gene networks associated with different environmental conditions. I used SVA to adjust for batch effects and "protected" different environmental conditions in the model as I don't want to lose information about them. However, when I used the "batch-adjusted" data to create networks, all the scale-free topology indices were negative (max = -0.05).
I tried generating the networks with softpower = 16 as it is suggested for sample sizes 20-30. However, only five networks were generated which is far less than I expected.
What might be the reason for getting negative scale-free toplogy indices?
Can combining studies from different environmental conditions cause the networks to be not biologically meaningful?

Welcome to the site. I have some questions before being able to provide an answer: Could you show (the code) how did you modify your data to remove the batch effects with sva? How many samples and genes do you have? Are all the RNA-seqs built against the same genome version? How did you call the network construction, full automatic or manual?

Whenever you have batch or other nuisance variables driving variation you're going to get odd networks. You tried using SVA, which is a good method to use, but it seems in this case that it's able to optimally handle this dataset. I also strongly suggest that you use pickSoftThreshold() to choose a softpower, since that may produce better results (though since you already have negative values I'm not holding by breath). In general, you're going to need to play with either SVA or filtering to remove the batch-related genes that are causing problems.

Thank you for the response, Devon Ryan. Removing the low expressed genes (log2(FPKM) < 1) before using SVA has partly resolved the issue. I still don't get a high scale-free topology, but now I have positive numbers (max = 0.6). The resulting modules look trustworthy.

@lence any update on the question? (you seem the same user, but your user id is different, do you need help to merge both accounts ?)


text:
Is there a way to import tables in PubMed Central to MySQL or other Relational databases?
I need to use these PMC tables from the medical journals for the purpose of integrating the data with a decision-making software. Though, I can see that one can access full-text for an article from PMC here (which includes the article table), but is there a way to extract the table from this XML document and export it to a local MySQL database or any other relational database?

Is connecting to a remote database an option? Also what data do you need? All of it? Only the paper? Only the abstract?

I need only the tables from the Published articles (mainly to access individual patient characteristics like Age, Gender, Marital Status, Disease History, etc.) Connecting to a remote database is surely an option.

If you can use R, consider using web scraping.
See here : http://bradleyboehmke.github.io/2015/12/scraping-html-tables.html

The tables are already present in the xml files that PubMed allows to download, but there is no way to restructure those xml components into the tables. For example, this study - (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5970438/) that extracts tables from xml and restructures them into MySQL tables using text-mining algorithm.


text:
Any input as to what to study in order to work into genetic engineering from math background?
Any input as to what to study in order to work into genetic engineering from math background?
I have troubles understanding what the "interface" would be between mathematicians and biologists/biotechnologists.
Wikipedia links the first model of an organism having built in 2011:
https://en.wikipedia.org/wiki/Computational_genomics#First_Computer_Model_of_an_Organism
And it based on statistically analyzing biology/chemistry papers. This reminds of a typical "stats research". So there's not necessarily direct interaction between the empirical and the formal. However, for new applications, it's given that there must be some interaction, in order to synthesize knowledge?

I don't understand the question. Do you seek advice to how to proceed in your career or do you want to know how does an "empirical" and a a "formal" science interact?


text:
Identifying Indels from Chromatograms
I have around 100 chromatograms (.ab1 files) from Sanger sequencing a genome at loci believed to have an indel.
I'm new to interpreting this kind of data in general, but I've read a bit on the general idea—mostly in guides like this. What I'm not able to figure out from the resources I've checked is how to identify insertions and deletions in a chromatogram.
I'm most concerned with  heterozygous indels, and these seem to be simpler. (And since my samples are (supposed to be) larger indels that are more rare and more likely to be het, I think.) A heterozygous indel, according to my understanding, would like something like consecutive heterozygous SNPs—overlapping peaks, maybe not with the same height, but with the same horizontal placement—until the end of the sample.
In other words, something like the top track in the open window in this image.

(Source: CodonCode.)
(I'm trying to do this manually, and not with packages like this for several reasons, including: I want to see how this is actually done, I want to do it as thoroughly as possible—I will probably check my results with an automated package—and I was asked to do it like this.)
How can evidence of (het) indels like this be distinguished from noise in the chromatogram? My guess is it includes things like

the overlapping peaks last until the end of the sample (and don't start at the beginning, but somewhere in the middle); and
the peaks are still at the same horizontal positions; i.e., one is not shifted right or left relative to the other?

Is this correct?

With my own chromatograms, I have a lot of doubt about distinguishing noise from indel. Here's an example, visualized with 4Peaks. (Full chromatogram here.)

On the left, the chromatogram appears high quality. There are sharp, regularly spaced, single peaks. On the right, there appear to be two sets of peaks lined up with each other. Does this constitute evidence for an indel? Does the homopolymer region—I've noticed in many of these samples, the putative indel seems to start in a homopolymer region—cast doubt on this?

I've also noticed that in some situations like this, there appear to be not two but three peaks.

Is one of these just noise and can be disregarded? Or does this cast doubt on the validity of the data—which is supposed to be coming from the diploid human genome?

If you want to call heterozygous SNPs and INDELs, you may try polyphred. I believe there are a couple of other similar tools, but I have forgotten their names.

How can evidence of (het) indels like this be distinguished from noise in the chromatogram?

Your two points are correct. Basically, at some position, you start to see two streaks of overlapping peaks. Once identified, the two streaks can be aligned to each other if you shift one of them. This is a strong constraint, but it is hard to tell by eye because it is not obvious how to assign peaks to streaks without running an automatic tool. Doing this is even more difficult given noisy traces like your last figure. I think you'd better run polyphred or similar.

the putative indel seems to start in a homopolymer region—cast doubt on this?

The INDEL mutation rate in long homopolymer runs is orders of magnitude higher than elsewhere in the genome. This has two implications: 1) you expect to see more homopolymer indels; 2) PCR polymerase is more likely to introduce errors in homopolymers––a higher fraction of homopolymer indels tend to be wrong.

Thanks! I'll look into polyphred and other tools. To be clear, what's your take on the indel in the homopolymer region—that the chromatogram is inconclusive itself? What can I do to verify it?

@Randoms Overlapping peaks alone are not strong evidence. You need to check if the overlapping peaks can be separated into two sets and the two sets can be aligned with a shift. I can't do this in my head. Even if there is an indel there, the chance of seeing a PCR error is still high. Run polyphred first. PS: finswimmer is right that overlapping peaks should be twice as "low". This makes the second figure less likely to be true.

I think you're right that running a tool is really the way to go. Thanks! I'm going to try Poly Peak Parser: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4525701/

to understand how to interpret sanger sequencing chromatograms I'd like to give a brief introduction to how this data are generated (sorry, if you already know this)
The target region is amplified by PCR. The following sequencing reaction is quite similar to a PCR, despite two important differences:

Only one strand is used as a template, by adding just 1 primer instead of two.
In addition two normal dNTPs, we have modified ddNTPS. Those are fluorescence labeled and the 3'-OH is missing, so there is no chance to add more dNTPS.

The result of this reaction are single stranded DNA fragments of different length. The last base is fluorescence labeled. So we have fragments that looks like this:
A
AT
ATC
ATCG
ATCGA

Of course there a many copies of these fragments.
Now these fragments are separated by capillary electrophoresis. Small fragments are detected first. We can differ the several bases by its fluorescence color.
Let's imagine a DNA strand like this:
ATCGA

What happens if we had on the second strand a SNP?
ATCGA
ATAGA

Now the sequence would detect two signal of equal intensity on base 3 - one for C and one for A. This have two effects in the chromatogram.

Two peaks of equal heights with different colors
The relative peak height is just half the height, compared to those sequences that have just one base at this position.

Let's see what happen in case of a heterozygous deletion. An alignment would look like this
ATCGA
AT-GA

But our capillary electrophoresis doesn't know about alignment. It only knows about fragment length. So this is what it will detect
ATCGA
ATGA

With the beginning of the deletion it will detect to different signals. Just like in case of the SNP, the relativ signal intensity is now just about the half.
With this knowledge let's go back to your initial questions:

How can evidence of (het) indels like this be distinguished from noise
  in the chromatogram? My guess is it includes things like

the overlapping peaks last until the end of the sample (and don't start at the beginning, but somewhere in the middle); and

They don't start at the beginning, if your indel is not right after the primer, correct. In case of an insertion or deletion you have strands of different size. So if your target region is quite small you may have luck, that you can see the end of the overlap at the end of the sequence, because the shorter strand reaches its end.

the peaks are still at the same horizontal positions; i.e., one is not shifted right or left relative to the other?

Yes. And one of the best things to distinguish between noise and true heterozygous variant is a relativ drop off of the signal intensities by approximately 50% compared to homozygous samples.
Concerning your second chromatogram: This is not an indel. As the other already said, this is a mistake made by the polymerase. Invitro polymerases doesn't like homopolymer stretches. They begin to add extra bases or "forget" one. As we do PCR this error accumulates. In my experience up to 9 time the same base is not a big problem, 10 only if you have luck, from 11 upwards your sequence are terrible. If you are interested in the sequence part behind this stretch, you have to sequence from the other direction.
fin swimmer

Thanks for your thoroughness! Interesting that you bring up the relative intensity of heterozygous peaks relative to a het peak at the same spot--I'm not quite sure I'll be able to use this, though, because I've only sequenced one sample, so I don't have multiple sequencing results that cover the same locus.
I'm interested in to know how you're confident that the chromatogram with the homopolymer does not show an indel—surely it's possible that there is an indel in a homopolymer region, right?

You see no significant drop-off of the average peak height and the additional base peaks are to small. That's why it's unlikely that there is a heterozygous indel. But as I wrote, the homopolymer stretch is to long to be sure. I wouldn't give any definitively answer here.

Downvote without comment? Not a good style :(

Honestly I would not try to call heterozygous indels with Sanger sequencing.  It's just not the right technology.
If you were trying to confirm a known indel that didn't flank a homopolymer stretch, and you had sequence from the other side that got messy at exactly the same point, maybe you could safely conclude there was something going on there, and possibly you could eyeball the exact right sequence of mixed peaks, but I don't see that what's what you have in your examples at all.  It just looks like noise.

Thanks for the input. By the other side do you mean the other strand—like forward and reverse Sanger sequencing—or the same strand, just on the other side of the variant?

For me sanger sequencing is still the gold standard for calling indels :) Compared to short-read sequencing is it quite easy to discover indels of any length.

But mixed indels?  Unambiguously no matter what the surrounding sequence?  If there is an indel in that poly-A region, It would be exceedingly difficult to prove.

Sanger sequencing produces a bit of messy signal at the start and a lot at the end of the sequence, but should have a region of clean signal in the middle. If the signal is getting mixed up in the middle, it suggests that there are sequences that are out of phase with each other. That could be a heterozygous INDEL, but it could also indicate that multiple template sequences have been amplified by the same primer.
Sanger sequencing is frequently carried out in both directions (forward and reverse complement), and if similar things are happening in both directions at the same place, it makes it more likely that there is a heterozygous INDEL and not something else funny going on. It's a very good idea to have homozygous templates for all expected variants, if possible. If the secondary peaks in the signal match what is observed in the homozygous variants, that can also be used as evidence. In the absence of control sequences, it is also possible to eyeball the peaks to find self-similarity.
Where triple peaks are seen, that could point to a region where there are two INDEL variants in close proximity, or multiple regions getting amplified at the same time, or even a single INDEL variant with more than two variant states (e.g. a homopolymer with 8 As that is sometimes 7 As, and sometimes 9 As).

Thanks for your response! Can you expand on what you mean by "positive and negative homozygous variant templates"?

Sorry, I was thinking of positive and negative controls; I've changed it to clarify as "homozygous templates for all expected variants"

There are specialized tools to deconvolute (separate) heterozygous InDels in a chromatogram trace using a reference genome. For instance, Indigo is a web application for this, the command-line application is called Tracy.

One important point that @gringer briefly touched on is that there is often off-target amplification when performing PCR on whole-genome DNA.
So, the noise in the chromatograms very well could be from another part of the genome and not related to indels at all.
The only way to ensure that this is not a problem is to clone your PCR products into a vector, transform into something like E. coli, grow up a single colony, cut the product with restriction enzymes, size-select, then sequence the correct-sized insert. Or you could just re-sequence your PCR product with Illumina or PacBio to accurately sequence the complete population of molecules -- indels and all. Check out:
Hebert, Paul DN, et al. "A Sequel to Sanger: amplicon sequencing that scales." BMC genomics 19.1 (2018): 219.
The safe thing to do is to not read too much into the noise of Sanger chromatograms unless it is a simple and very clean SNP.

Before making such a big stunt with cloning, NGS, PacBio etc. one should check the  specificity of the primer, make changes to the PCR protocol, make sure that the clean up steps work, ... If you are working on region that have no homologues, off-target amplification isn't normal.