text:
Smallest group size for differential expression in limma (bulk RNA-Seq)
I am reading Smyth et al. (ref. 1). I want to run differential expression analysis on a bulk RNA-Seq dataset in which each group is composed by 2 samples. In the paper previously cited it is written that:

Genes must be expressed in at least one group (or in at least three samples across the entire experiment, where three was chosen as this is the smallest group size) to be kept for downstream analysis.

Is it possible to use limma DE analysis also with groups composed by only 2 samples? NB. It is possible that in this particular dataset the smallest group size is 3.
If not, which alternative should I use?
Update
I have 17 samples. My idea is to test one group vs the rest of groups (all together). So at least the reference group would be composed by > 10 samples. In this case, what can be a feasible analysis for DE?

Smyth, G. K., Law, C. W., Alhamdoosh, M., Su, S. & Ritchie, M. E. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. F1000Research 5, 1408 (2016).

The problem with most of the methods is that use the gene's variance for each group, which can't be calculated (reliable) when the sample is <= 2. Also, statistically, it would have extremely low power, so the conclusions couldn't be trusted much.
You could also calculate the "raw" fold change by yourself (ie without the variance estimation and adjusting). If you had one more sample per group you could use DESeq2 which was thought for such cases. But from the paper:

However, if there are two or fewer replicates for a condition, these samples do not contribute to outlier detection, as there are insufficient replicates to determine outlier status.

Which makes me unsure if it will work well for just two samples as it won't help to detect outliers genes.

Thanks. I updated my question. In particular, I was thinking about using your suggestion and just take the raw fold-change. Is it correct to take genes that are considered outliers (e.g. > 3 StD) when comparing the group to the rest of samples?

The problem with such a few number of samples is that it is hard to know what is an outlier and what it is "normal". If you compare 2 vs 15 you still have two problems. 1) you still need to calculate the variance for the group of 2, which has the same problem as initially and 2) the comparison will be meaningful? If you considered them to be different groups is with a reason, if you mix them it will help you to understand what are the differences between them?

You could show the samples and the relevant data about these samples (which groups, do they belong, if they are biological or technical replicates or the differences between the conditions, ...) to help you better.

The dataset I'm using is this one on GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE109125). There's no citation yet and I didn't find in the metadata if they are biological or technical replicates. However, I think they are biological replicates. The group is the first part of sample name (before the '#' sign). The full normalized matrix can be downloaded here: ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE109nnn/GSE109125/suppl/GSE109125%5FGene%5Fcount%5Ftable%2Ecsv%2Egz

This may sound obvious, but can you provide some context in the answer on why the variance cannot be calculated reliable when the sample is <= 2?

@gc5 With less than 2 you don't have variance, and with just 2 it is not good enough, as it is defined by the difference between the points and the mean, and if you add another point you add 33% more data! But that would be a good question for stats.SE, who could give a more reasonable and detailed answer than me :\

While I share Llopis' concern about estimating variance from 2 samples, the statement you quoted is about avoiding false positives from genes that are only expressed in a few samples. It's fairly common to exclude genes that expressed in fewer samples than the smallest group even if the number of samples per group is much higher than 3.

Good point, thanks. I also think the sentence meant that in this specific case.


text:
How to split bedfile into non-overlapping regions and compute aggregation function on duplicate segments
I have a bed-like file in which I have several overlapping regions with associated features. I want to split the overlapping regions into disjoint regions in order to process the file (e.g. with pandas), group by regions and compute an aggregate function.
The input data format is the following (overlapping are marked with *):
chr9    8855925    8986425    -1  *
chr9    8855965    8986325     0  *
chr9    8893857    9022105    -1

I want to split into disjoint intervals, then I want to groupby per-interval and compute an aggregation function (e.g. the mean of the last field):
chr9    8855925    8855965    -1
chr9    8855965    8893857  -0.5
chr9    8893857    8986325 -0.67
chr9    8986325    8986425    -1
chr9    8986425    9022105    -1

I was thinking of using bedtools and pandas to aggregate. However, I was unable to find the right tool to do this job in the bedtools/pybedtools documentation.

To make disjoint intervals, you could use BEDOPS bedops --partition, piping to bedmap --mean to get the mean signal over disjoint regions.
Starting with the input bedgraph file, convert it to five-column BED with GNU awk, putting the signal in the fifth column per UCSC convention:
$ awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }' /tmp/in.bedgraph | sort-bed - > /tmp/in.bed
$ cat in.bed
chr9    8855925 8986425 .   -1
chr9    8855965 8986325 .   0
chr9    8893857 9022105 .   -1

The disjoint set from these intervals would look like this:
$ bedops --partition /tmp/in.bed
chr9    8855925 8855965
chr9    8855965 8893857
chr9    8893857 8986325
chr9    8986325 8986425
chr9    8986425 9022105

We can pipe these intervals to bedmap --mean with the original five-column BED file as a "map" file, calculating the mean signal over the disjoint intervals:
$ bedops --partition /tmp/in.bed | bedmap --echo --mean --delim '\t' - /tmp/in.bed > /tmp/answer.bedgraph
$ cat /tmp/answer.bedgraph
chr9    8855925 8855965 -1.000000
chr9    8855965 8893857 -0.500000
chr9    8893857 8986325 -0.666667
chr9    8986325 8986425 -1.000000
chr9    8986425 9022105 -1.000000

The --mean option calculates the arithmetic mean. Other aggregation functions are available to process the signal or score (fifth) column in map BED files, like --stdev, --median, --mad, etc.
When calculating statistics, the default precision is six digits. The --prec <n> option can be used to specify desired precision in the score result.
See bedmap --help or the online documentation for a full description.
To avoid making and cleaning up intermediate files, bash process substitutions could be used to generate a one-liner:
$ bedops --partition <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') | bedmap --echo --mean --delim '\t' - <(sort-bed in.bedgraph | awk -vOFS="\t" '{ print $1, $2, $3, ".", $4 }') > answer.bedgraph

This can be useful for cleaner pipelines or with processing very large, whole-genome scale inputs on systems with limited disk space.

Thanks for the detailed answer. I got confused with the numbers, I should have used very simple segment start and end. Your answer has the expected solution. I'll update my question.

On macOS, use gawk instead of awk since the -v option is invalid. With homebrew brew install gawk.

source: Stack Exchange

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

sort: votes


text:
Scaling by linear regression against the number of reads
I am trying to build the preprocessing pipeline presented in The Tabula Muris Consortium et al. (pp).
It is a pipeline to preprocess single-cell sequencing data. There is one step that is not clear:

Counts were log-normalized (log(1 + counts per N)), then scaled by linear regression against the number of reads (or UMIs), the percent of reads mapping to Rn45s, and the percent of reads to ribosomal genes.

I understand the first part (I assume that log in this context is log2), but I need help on understanding how to scale by linear regression against the number of reads, the percent of reads mapping to Rn45s, and the percent of reads to ribosomal genes.

Have you contacted the paper authors? They would be able to give you a better answer, and your comments/questions would help them improve the paper when it is properly published.

@gringer You are right, I posted here because I thought there was some standard procedure.

The linear regression is straight-forward enough, but regressing against only three values (if I'm interpreting that correctly) is a recipe for disaster (it sounds vastly worse than even ERCC spike-ins...and those aren't exactly ideal).

@DevonRyan ignoring by now the small number of regression values, how do you scale using the linear regression? I am unable to see how the linear regression is used in this case. Can you elaborate more (maybe with an answer)?

Did you figure how did they preprocess the tabula muris data?

@yuqi_yuqi see https://bioinformatics.stackexchange.com/a/3225/1771

I don't know if this question has been solved already, but what they try to do is equalize the depth of sequencing for each cell. Therefore, they scale for the total number of reads. If you regress out (via linear or negative binomial regression) the differences in the number of reads per cell, you end up with cells that have been sequenced with the same depth.
In my opinion they use the ribosomal genes in the same way. They are considered to be a kind of housekeeping genes, which you can use to equalize the sequencing depth.
And yes, the log in this context is log2, used to get fold-change values instead of counts.

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

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


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

While it is obvious that scRNA-seq data contain lots of zeroes, I couldn't find any detailed explanation of why they occur except for short notices along the lines of "substantial technical and biological noise". For the following text, let's assume we are looking at a single gene that is not differentially expressed across cells.
If zeroes were caused solely by low capture efficiency and sequencing depth, all observed zeroes should be explained by low mean expression across cells. This however does not seem to be the case as the distribution of gene counts across cells often has more zeroes than would be expected from a negative binomial model. For Example the ZIFA paper explicitly uses a zero-inflated negative binomial distribution to model scRNA-seq data. Modelling scRNA-seq as zero-inflated negative binomial seems widespread throughout the literature.
However assuming negative binomial distribution for the original counts (as measured in bulk RNA-seq) and assuming that every RNA fragment of the same gene from every cell has approximately the same (low) chance of being captured and sequenced, the distribution across single cells should still be negative binomial (see this question  for related math).
So the only remaining possible cause is that inflated zero counts are caused by PCR. Only non-zero counts (after capture) are amplified and then sequenced, shifting the mean of the observed gene counts away from zero while the pre-PCR zero counts stay zero. Indeed some quick simulations show that such a procedure could occasionally generate zero-inflated negative binomial distributions. This would suggest that excessive zeroes should not be present when UMIs are used - I checked one scRNA-seq dataset with UMIs and it seems to be fit well by plain negative binomial.
Is my reasoning correct? Thanks for any pointers.
The question
How can we distinguish between true zero and dropout-zero counts in single-cell RNA-seq? is related, but provides no clues to my present inquiry.

It may be necessary to distinguish between methods that use unique molecular identifiers (UMIs), such as 10X's Chromium, Drop-seq, etc, and non-UMI methods, such as SMRT-seq.  At least for UMI-based methods, the alternative perspective, that there is no significant zero-inflation in scRNA-seq, is also advocated in the single-cell research community.  The argument is straight-forward: the empirical mean expression vs. dropout rate curve matches the theoretically predicted one, given the current levels of capture efficiency.
Examples
Svensson Blog
A couple of blog posts from Valentine Svensson argue this point rather pedagogically, and include citations from across the literature:
Droplet scRNA-seq is not zero-inflated
Count-depth variation makes Poisson scRNA-seq data negative binomial
bayNorm
There is a more extensive preprint by Tang, Shahrezaei, et al. (BioRxiv, 2018) that claims to show a binomial model is sufficient to account for the observed dropout noise.  Here is a snippet of a relevant conclusion:

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

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

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

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

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

I highly recommend giving this a read, especially because it nicely situates other common generative models (e.g., binomial, Poisson) as valid simplifying assumptions of the multinomial model.
It should be noted that this same group previously published a work (Hicks, Irizarry, et al. 2018), mostly focused on non-UMI-based datasets (SMRT-seq), where they showed evidence that, relative to bulk RNA-seq, there was significant zero-inflation.

One more preprint arguing that zero inflation is not a good model: https://www.biorxiv.org/content/10.1101/477794v1.abstract

The Biostars thread turned out helpful. The most interesting possible cause, not mentioned in the Ian Subery's answer, is that due to bursty nature of transcription, the true distribution of transcript counts across cells can be bimodal with a peak at zero even assuming a simple model of transcription such as the random telegraph model. See for example Dattani & Barahona 2017 for a more detailed discussion.
However, the random telegraph model predicts that the even with peak at zero, there should usually be non-negligible probability of having small, but nonzero counts, which is not always the case in zero-inflated scRNA-seq data: see for example the Isl1 gene in Retina dataset (click Explore -> search for "Isl1")

I know of no references for this, but in general, I would say that your reasoning is sound. I would just add that in contrast to what I suspect you have simulated, not all transcripts are equally likely to be captured and amplified. We don't really understand what the determinants of this are, but for example, GC content is definitely related.

Thanks for the note. I understand there are betweeen-gene differences due to GC content etc. but I am only looking at one gene at a time so it shoud not matter. Do you believe that the sequencing results could be noticeably influenced by within-gene differences - e.g. because a specific fragment of a long transcript is much better captured/amplified than others letting us observe large non-zero counts only in cells where this fragment was captured and zeroes otherwise?

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


text:
TruSeq Exome Kit bed file old version
I have to run an analysis using Whole Exome Sequencing (WXS) data prepared with the TruSeq Exome Kit (2012).
I need the .bed file for the kit and I found it here.
However, this file was posted on Nov. 16, 2015, while the kit used in my case is from 2012. I searched but I didn't find the bedfile for that specific version online.
Will I run into problems if I use this version, and in this case, where can I retrieve the bedfile for 2012's version?

Care to explain why the downvote?

I haven't down voted you, but you don't show what have you done to understand what problems can you have when using an older file, nor the places you have searched for the bed files of 2012. Or maybe it was from biology.se, who knows ?

I usually search before writing a question, so also in this case I searched online to understand what could be the difference. Since online I couldn't find anything I wrote this question. Also, the fact that I put the reference to the Illumina kit page (the only page in which I could find reference to the bedfile of the kit) disproves that I didn't put the places I have searched. So I don't see the point of downvoting, since I: clearly stated the problem, provided details about what I've done, provided links to what I've found.

Well, those are imaginary internet points, don't mind them too much!

Actually I don't mind them :) but since I put some effort to search and write I couldn't understand why the downvote.

Usually it is more problematic to work with older versions of data with newer protocols than the other way around.
In this case, I think that the features annotated in the (newer) bed file should be more accurate than the older one. Unless you specifically want to compare your experiment with a previous experiment performed with a bed file from 2012 you'll be fine using the newer one.
However make sure that you are correctly using the right file and that you report (or record) each version of protocol and data used in case some doubts arise.


text:
Estimating BAM file from compressed fastq file size
Is there a way to estimate the size of a BAM file will have after mapping with BWA?
The input file are two mates fastq files, compressed with gzip, each one about 70G.

A BAM file should have a roughly similar size to a compressed FASTQ file, because they're both compressed files that contain similar information. BAM will be slightly larger, because it contains information about mapping and differences from the reference.
If you are concerned about space and have a reference sequence available, you can store mapped reads as CRAM format, which will save a lot more space.

Could also be smaller if you remove duplicate reads. I doubt it will ever be significantly different in size from the original fastq file, but a +/- a few hundred megabytes (more for larger files) or so is to be expected.

A sorted BAM will be (potentially substantially) smaller than an equivalent FASTQ file, even when containing mapping information.


text:
Genotype quality in VCF for haploid calls
I cannot find in the VCF specs what is the default value of Genotype Quality GQ for haploid calls. In such case, is the default the missing value .?
From the specs:

If any of the fields is missing, it is replaced with the MISSING value. For example if the FORMAT is GT:GQ:DP:HQ
then 0 | 0 : . : 23 : 23, 34 indicates that GQ is missing. If a field contains a list of missing values, it can be represented
either as a single MISSING value (‘.’) or as a list of missing values (e.g. ‘.,.,.’ if the field was Number=3). Trailing
fields can be dropped, with the exception of the GT field, which should always be present if specified in the FORMAT
field.

From this I assume that there is no default value, GQ might be dropped from the FORMAT field, correct?

Can you give an example? Why would there be a default quality?

Thanks @terdon, I probably misread the specs. I updated the post with a part of the specs that indicates the GQ field can be dropped from the FORMAT field, if I understood correctly.

Thanks, but I'm still a bit confused. You are already quoting the relevant section that explains that yes, . is what you put when one of the samples or alleles in the file is missing a GQ value.  Are you asking if GQ is a required field for a valid VCF file or if you can have a VCF without GQ?

Sure. My doubt derives from the fact that I haven't understood if, for haploid calls, the GQ equals the QUAL field of the variant, or if it's dropped because we only have one allele.


text:
Assess ploidy estimation through SNV and CNV results
I want to assess ploidy estimation (i.e., diploid vs tetraploid) in a set of samples for which I have both results from the variant caller (i.e., SNVs) and from the copy number caller.
What are the possible approaches that can be used in this case? I know that in some in cases it might be possible to use VAF to assess ploidy estimation, but I haven't found exactly how this (possibly simple) approach is defined.
Note: I already have a caller that I am testing. My question is more related to exploratory data analysis or manual inspection of the caller's results.
VAF distribution examples might be similar to what is done in CNAqc (image below extracted from their online tutorial):

What is the underlying data quality? Would you be able to go in and visually confirm alleles in a single sample in IGV or whatever? Approaches such as here use allele frequency at heterozygous SNPs in single samples to try to address this question, so I think that's the relevant question. I have to say that I don't understand the figure or the tutorial at first glance, it appears to be looking for de novo mutations (?) and seeing their frequency.

Thank you. I have updated the question regarding the tutorial - I'm only interested in plotting the VAF distribution for the different genotypes, that is the first example I could find online to explain one of the tools for the analysis (sorry it that confused you). I will read the paper you linked, thank you.


text:
When can we have duplicated variants in a VCF file, excluding artifacts?
I read these two questions on duplicated VCF variants: 1, 2.
However, it is not entirely clear to me when duplicated variants can arise (excluding merging artifacts). I want to understand if duplicated variants should be always considered artifacts, or if there are specific use cases where they are needed.

Do you have a specific problem in mind? There are multiple reasons why duplicate variants may appear (especially when merging datasets), and Stack Exchange doesn't really work well when there are multiple "correct" answers for the same question.

I don't have a specific problem in mind. I wanted to understand when to expect duplicate variants, since I noticed that some algorithms filter them. Although I understand merging artifacts, I wonder if there are additional corner cases where these duplicate variants appear. The accepted answer can be a community wiki where multiple "correct" answers can be collected, I believe.

Okay, probably better asking this question on Reddit r/bioinformatics then.


text:
Handling of sex chromosomes in male when computing CNV copy ratio
Reading from CNVkit documentation, Mathworks documentation, Bionano documentation, for CNV calling, it seems that the common practice for array is to compute the linear copy ratio by normalizing against the expected depth for the reference ploidy, i.e.,:
$$CR = \frac{D_{OBS}}{P_{REF}}$$
where $D_{OBS}$ is the observed depth, and $P_{REF}$ is expected depth for the reference ploidy for a specific segment (i.e., expected depth for 2 copies for females, and males outside sex chromosomes or in PAR regions).
Is that correct? Are there alternative approaches used in array copy number calling?


text:
Drawbacks of upper quartile normalization for scRNA-seq data
I would like to use Upper Quartile normalization for scRNA-seq data defined as:

The upperquartile (UQ) was proposed by (Bullard et al. 2010). Here each column is divided by the 75% quantile of the counts for each library. Often the calculated quantile is scaled by the median across cells to keep the absolute level of expression relatively consistent. A drawback to this method is that for low-depth scRNASeq experiments the large number of undetected genes may result in the 75% quantile being zero (or close to it). This limitation can be overcome by generalizing the idea and using a higher quantile (eg. the 99% quantile is the default in scater) or by excluding zeros prior to calculating the 75% quantile.

In particular, I'd use the 99% quantile to avoid normalizing by small number or zero. Is there any other drawback for using it on scRNA-seq data?

Bullard, James H, Elizabeth Purdom, Kasper D Hansen, and Sandrine Dudoit. 2010. “Evaluation of Statistical Methods for Normalization and Differential Expression in mRNA-Seq Experiments.” BMC Bioinformatics 11 (1). Springer Nature: 94. doi:10.1186/1471-2105-11-94.

My primary concern with using the top ~1% or so in upper quantile normalization is that it's going to be prone to the same robustness issues that RPKM/FPKMs run in to. That is, if for whatever technical reason you have to have a fair bit of variability in a couple very highly expressed genes (typically rRNAs, but one can imagine other genes) and the set of genes going into the UQ normalization is composed primarily of those then the normalization results will simply reflect irrelevant technical variability.
Of course, this all assumes that it's even appropriate to normalize scRNAseq data with a standard method like this. In some cases it is (e.g., when comparing control vs. treatment with a largely homogeneous cell population), but if you're dealing with very different cell types then I worry that the standard method is going to wash away things like differences in bulk RNA content between cell types. Presuming there are spike-ins, they should have the same normalization applied and then checked to ensure there's nothing greatly amiss. Having said that, if you're going to be doing tSNE and want to look for cell types, then odds are good you'll still be able to see most of them (after all, it's unlikely that RNA quantity is the only difference between them, though not removing that difference would make things easier).
The problem of global shifts in expression
The concept of global shifts in expression is not unique to scRNAseq, in fact the most common example of this is Myc up/down regulation in relation to cancer. To explain the issue, let's suppose that you have the following (very small) matrix of counts:
       sample1 sample2
gene1        1       2
gene2       10      20
gene3       15      32
gene4      100     198

As you can see, sample2 has roughly twice as many reads for each gene as sample1, so any normalization is going to cut all of its counts in half (or the reverse for sample1, or perhaps a mixture of the two...depending on the method). That's usually fine, since the assumption with all common normalization methods is that differences like this are simply due to differences in sequencing depth between the samples. But what happens if sample2 actually has twice as much RNA as sample1? Well, then you've just normalized away all of the differences between the samples (I encourage readers to go through the literature on Myc, since that's a very nice example of this, where a mostly-global shift in expression led to apparent down-regulation of genes with constant expression). Such a situation doesn't occur very commonly in bulk RNAseq, but when it does you need to use something like spike-ins (dosed to match cell number or DNA quantity) to maintain the global difference.
If this is only rarely an issue with bulk RNAseq, then why is it more likely to be an issue in scRNAseq? The answer is that bulk RNAseq benefits from the averaging out of many different cell-types, such that ensemble averages/medians don't tend to change much due to treatments...but this isn't going to be the case at the single-cell level, since cells can vary vastly in size (compare a sperm and an oocyte...every RNA should be differentially expressed between them). This difference in base RNA content between different cell types will be washed away by common normalization methods, which can lead to either losing differences or creating differences in genes where there shouldn't be (again, see examples of this in the Myc literature). Whether this really becomes an issue will depend on what your goals are with your data and how much of a global expression difference you actually have and the type of analysis you want to perform.
At the end of the day, just remember that common RNAseq normalization methods assume that (A) there's no global expression change between samples and that (B) if there are very many DE genes that the direction of change is not heavily skewed in one direction. If either of these assumptions are violated then you can run into issues.

Thanks. Could the primary concern go away by using a lower quantile (such as 90%)? Does removing rRNAs prior to UQ normalization mitigate this issue? Regarding spike-ins, I removed them after filtering per Spike-in/Endogenous RNA ratio (retaining only cells with low ratios). In theory, the dataset I'm using should be composed by a largely homogeneous cell population, can you elaborate more on how this method is going to wash away things like differences in bulk RNA content between cell types?

The more genes you use the more robust it'll be, so yes that's helpful. I'll update my answer to make mention of the issues with normalizing when there are global changes.

Thanks, now it's clear. I think I will normalize using spike-ins.


text:
TCGA bam file duplicates (aligned reads) and related last part of TCGA barcode
I am running an analysis on TCGA bam files (GCRh37/hg19) of whole exome sequencing data.
I usually find one bam file for each sample. However, at least one sample has two bam files, in which the last part of the name (before the extension) has the suffix .1 or .2 or .3...
For example, C500.TCGA-BL-A3JM-01A-12D-A21A-08.1.bam and C500.TCGA-BL-A3JM-01A-12D-A21A-08.2.bam.
I didn't find an explanation for this last suffix. Are they duplicates, or.. ? Is there some wiki about the last part of the filename?

If you check the bam list of TCGA, you'll see that they use different numbers for the bam files (maybe different settings). What the difference is, I don't know, but if you want the file that they used for their paper, it is C500.TCGA-BL-A3JM-01A-12D-A21A-08.2.bam, the other one is not in the list.


text:
ComBat for batch effects removal on scRNA-seq data
Is it possible to use sva's ComBat for batch effects removal on scRNA-seq data?
Is there any difference between RNA-seq and scRNA-seq data which doesn't allow to use ComBat on single-cell data? (I am thinking about extreme sparsity, > 80%). If so, are there alternatives?

see here for a simiar [unanswered] question.

@gringer Indeed, that was another question of mine :)

Look at this recent paper that uses ComBat on scRNA-seq data for batch effect removal and states that it "successfully does so".
I also suggest that you check out this publication on Distribution Matching Residual-Nets. Authors evaluated their method also on scRNA-seq data and thus it may be something you are looking for. I personally played a bit with their code and found it straightforward to use.

There are competing claims regarding how well ComBat works on scRNA-seq datasets. A recent paper from John Marioni's introduces a mutual nearest neighbors method that seems to outperform ComBat in at least some relevant scenarios. In general, it's best to look at tSNE or other diagnostic plots after batch correction to see if the results seem reasonable (this advice holds generally).


text:
Detect transcript isoform abundance for a specific gene in scRNA-seq
I want to detect the count of isoform transcripts for a specific gene in scRNA-seq data. Data is coming from cells of Mus Musculus.
For transcript isoforms I mean the different alternatives provided by known annotation. For each one I need the total transcript count in each cell.
By looking at the count matrix, alternative transcripts are not defined (see for instance this dataset).
Should I come back to the fastq? Which tool may I use to produce a count matrix with counts per each alternative transcript?

Transcript quantification is a difficult enough problem as it is, when you add the extra difficulty of going from the low read numbers available in scRNAseq if gets even more difficult.
Added to this, many scRNAseq techniques (including 10X, which is what the linked dataset is) only capture the 3' end of the gene, so you can distinguish transcript with different 3' end, but the data simply doesn't exist for doing much more with standard tools.
Lior Pachter does claims that there is benefit in looking at "transcript compatibility counts" rather than just gene counts, and that you can use his Kallisto tool for this, and the these 3' end tagging tools do give you some data for the rest of the transcript.
You could also isolate transcripts with differing three 3' ends and map reads to these using the something like the UMI-Tools pipeline.
Alternatively you could look at things on a junction level - you won't capture all junctions, but you may be able to look at differential junction usage for those you do capture using something like SUPPA or MISO. For this I would use UMI-Tools to get de-duplicated BAM files, split them on cell identity and then pass them into your tool of choice. I don't really how much you will capture.
The last thing to say is that this might be a space worth watching for future developments (cough).
Data from some other single cell methods, like SMART-seq2 or C1 will sequence the whole transcript and you may be able to apply standard tools such as Salmon, Kalisto or RSEM to get transcript usage for them, although as I said, transcript level quantification is noisey, single cell data is noisey, single-cell level transcript data will be very noisy.
Whatever, you are going to have to go back to the raw sequencing data.
Disclosure: I am the author of UMI-Tools.

Looking for alternative splicing in scRNA-seq is difficult, but there is at least one tool that is specifically designed for it:
DISCO. The site does not have too much info, but there is a presentation about it. It uses MISO output as the input, which in turn uses SAM.


text:
Are mitochondrial genes to exclude in scRNA-seq such as ribosomal genes?
In this answer, it is stated that ribosomal genes should be excluded prior to normalization in scRNA-seq as contaminants.
Do mitochondrial genes have to be excluded as well? I plotted the top 50 expressed genes for a specific dataset and they tend to appear often (for example MT-ATP6). My assumption is that, given that they work for mitochondrial function and may be highly expressed, they can dilute the signal of gene differential across cell types but expressed at lower levels. Is this biologically sound?
Additionally, in this course, mitochondrial genes are used to filter cells when they contribute above a certain ratio of total RNA of single cells. However, I could not find an explanation for this filtering. What is the underlying process this procedure is used to account for?

It's worth distinguishing between the ribosomal RNA gene (Rn45s) and the many genes that code for ribosomal protein (which start with RPS, RPL, MRPS, or MRPL). mRNA's for ribosomal protein subunits are an intermediate stage used for protein production. You might catch a few lucky ribosomal RNA (Rn45s) copies that were transcribed recently and still haven't been processed, but depending on your RNA-seq protocol, there's a chance they are acting as a functional part of the ribosome rather than an intermediate stage.

According to Ilicic et al. (2016), on upregulation of mtRNA in broken cells:

There is an extensive literature on the relationship between mtDNA, mitochondrially localized proteins, and cell death [34, 35]. However, upregulation of RNA levels of mtDNA in broken cells suggests losses in cytoplasmic content. In a situation where cell membrane is broken, cytoplasmic RNA will be lost, but RNAs enclosed in the mitochondria will be retained [...]

According to this, I think that mtRNA levels relative to endogenous RNA can be used as a control for low quality (broken) cells, in a similar way to ERCC spike-ins.

Ilicic T, Kim JK, Kolodziejczyk AA, Bagger FO, McCarthy DJ, Marioni JC, Teichmann SA. Classification of low quality cells from single-cell RNA-seq data. Genome biology. 2016 Feb 17;17(1):29.

This is correct, though I'd add that while the abundance of mtRNA is indicative of broken cells, analysts should ensure that filtering done based on "percent mtRNA" should be guided by the dataset at hand, not a hard threshold.  I have seen several analysts exclude all cells with more than 5% mtRNA just because that's what is done in a Seurat tutorial.   In practice, the threshold is (a) arbitrary, (b) can vary, and (c) should be informed by visual inspection.

Mitochondrial genes do not always have to be removed, though you may need to do so for a variety of technical reasons.
Mitochondrial reads are innately different from rRNA due to rRNA typically being excluded during normal library preparation (so any of it getting carried over is essentially a contaminant). Some scRNAseq programs (e.g., RaceID) can use sampling at various stages, which can break if you have a lot of your signal going to mitochondrial transcripts. Then it's not that they should be excluded because they're from mitochondria, but rather that the results are heavily affected by the presence of VERY highly expressed genes and your results are suspect if you don't remove these RNAs. When your methods aren't affected by that (e.g., you use more robust scaling methods) or you don't have super high expression of mtRNA then that's not explicitly needed.
Further, sometimes you very much want to keep mtRNAs, since you might be interested in metabolic changes. Since I work in an institute partially dedicated to immunology, we often look at things like metabolic changes involved in immune activation and removing mtRNAs would hinder that. But if you're a priori not interested in metabolism then excluding genes you won't care about isn't unreasonable.


text:
Why is ribosomal RNA difficult to remove even with Poly(A) selection?
In this answer (actually in a comment), it is stated that:

As you've noticed from your own analysis, the ribosomal genes have quite variable expression across cells. They're expressed everywhere, and quite difficult to completely remove from a sample (even with polyA selection), with success of removal depending on things like the amount of mechanical damage that the RNA has been exposed to.

Which are the most relevant issues that make Poly(A) selection fail? Looking online, Poly(A) is suggested often but I didn't find an explanation for the different ways in which it can fail.

We've found ribosomal RNA to be less of a problem with sequencing that depends on polyA, which suggests the issue might be in the library preparation, rather than the selection.
Many polyA RNA library preparation methods involve amplification, rather than selection, which means that existing transcripts that are present in very high abundance (such as rRNA) will still show up in the mixture.
But polyA-tailed RNA isn't the only RNA present in a sample. For a comprehensive look at RNA, the whole-RNA sample preparation methods can't use any form of selection to identify target molecules. In this case, it's typically recommended to use a form of ribosomal RNA removal, where a probe to ribosomal RNA is attached to magnetic beads (or a column, or similar), and used to fish out out ribosomal RNA, leaving the remainder behind. Unfortunately those probes don't work perfectly, either due to the ribosomal RNA being damaged at the probe site, or the probe not matching the ribosomal RNA molecule, or the RNA just not getting into the right physical state to bind with the probe. All these things mean incomplete removal of ribosomal RNA, ending up with a situation where the highly-abundant ribosomal RNA still pokes through in the reads.


text:
Regressing out unwanted sources of variation in single cell RNA-seq data
I have a dataset of single cell count data and I want to regress out the variation caused by the number of UMI's and the percentage of mitochondrial genes.
I know that count data is discrete data, and commonly follows a negative binomial distribution. To regress out these two confounders, should I use a linear model or GLM (poisson/negative binomial)? And how can I determine the optimal choice?

What do you mean by single cell count data? If you mean mRNA counts, then your data are compositional and should be treated as such.

yes mRNA counts, as I mentioned in the title. What regression technique is applicable for compositional data?

If I understand you correctly, "I want to regress out the variation caused by the number of UMI's and the percentage of mitochondrial genes" means that you want to mitigate the constrained geometry arising from multinomial sampling (aka the unit-sum problem). There are multiple transforms that can alleviate certain effects of compositionality (alr, clr, iqlr) or project the data into an unconstrained Euclidean space (ilr). The choice depends on subsequent analyses and their properties. ALDEx2, a widely used RNA-seq package, uses clr and iqlr. Take note, that they all assume certain properties.

Cleaning data before doing analysis can be more important than the method of data analysis. For UMIs, there's an obvious cleaning step that can be done: it would be better to filter out duplicates (prior to generating counts) than to try to incorporate that information in the count analysis.
I'm not familiar with single cell analysis, but I am aware that DESeq2 can be adapted to work with single cells, and incorporating additional variation is just a matter of adding an additional covariate into the DESeq2 formula model.

Thank you for your answer. In the meantime i found out that a linear regression is commonly used, but a negative binomial would be the favorable way to go.

Seurat has as part of its protocol a step where you filter based on UMI counts and percent mitochondrial
http://satijalab.org/seurat/pbmc3k_tutorial.html

Yes, but it uses a linear regression as a default. My question is specifically why a linear or a negative binomial regression.

If you just want to filter UMI related problems and cells with high mitochondrial to nuclear ratio, maybe there is no need to use modeling.

For UMIs: You can filter UMIs that are not supported by at least x reads (for instance 3 reads) and that probably arise from sequencing errors.
For high mitochondrial to nuclear ratio: You can compute (number of molecules of mitochondrial genes / number of molecules of nuclear genes) and plot this ratio for all the cells and try to put a threshold value for cells that have high values for this ratio. Thus, filter these cells that behave strange compared to the rest.

Thank you for your reply. I know these are standard threshold protocols, mostly dealing with badly sequenced cells or dying cells, however, this is not a correction for amplifcation bias. So I don't think this answer is what i'm looking for.

@DCZ I was following the Seurat tutorial for clustering cells (http://satijalab.org/seurat/pbmc3k_tutorial.html) when I found something very similar to your question. Literally, in the removing unwanted source of variation step: "In this simple example here for post-mitotic blood cells, we regress on the number of detected molecules per cell as well as the percentage mitochondrial gene content." It seems that they used a linear model, extracted from Buettner et al., 2015 (https://www.nature.com/articles/nbt.3102). Is this what you are looking for?. They even have an implemented function.

That's correct. The Seurat tutorial indeed does this, but the scaleData function uses a linear regression by default. As count data follows a negative binomial distribution, this seemed a peculiar choice to me.


text:
Database of copy number alterations in different cancer types
Is there a database of somatic copy number alterations (CNAs) in different cancer types? That is, annotations of which CNAs are more commonly seen in different cancer types.

I don't know if there's a database that does exactly what you want, but there are some places that might help you figure this out, especially if you already have a list of CNAs/genes/regions in mind.
ICGC is the first that comes to mind, as it has samples from many types of cancer and copy number data for most of them. Depending on what you want though, you might have to apply for access to download the data you want. It's rather annoying, and their data download tools are atrocious to use, but they have a lot of data.
DGV might be another option, though their data browser is rather dated feeling. They publish a lot of gold-standard germline CNVs and such.
COSMIC is a large database for cancer-associated mutations, including CNVs, so it might be worth a look for you.


text:
Least present k-mers in the human genome
What are the least present k-mers in the human genome at different sizes?
Starting with k=4 and going up in size until k=10, what are the k-mers least seen (or not at all) in the human genome? I am only interested in the reference human genome, so I am discounting SNPs/Indels in the population.
If this is not pre-calculated somewhere, what tool is recommended to use that starts with the GRCh38 reference as input?

If you are interested in least present k-mers, you might be interested also in the concept of minimum absent words (these two things are related). There are many papers on this topic, see, e.g.: http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0029344. Another phrase to google is "word avoidance".

You can use the Jellyfish software to calculate the k-mer profiles up to length 31.
From the instructions in the user guide:
The basic command to count all k-mers is as follows:
jellyfish count -m 21 -s 100M -t 10 -C reads.fasta

To compute the histogram of the k-mer occurrences, use the histo subcommand (see section 3.1):
jellyfish histo mer_counts.jf

To query the counts of a particular k-mer, use the query subcommand (see section 3.3):
jellyfish query mer_counts.jf AACGTTG

To output all the counts for all the k-mers in the le, use the dump subcommand (see section 3.2):
jellyfish dump mer_counts.jf > mer_counts_dumps.fa

I am trying jellyfish now.

I tried jellyfish, and got the list of kmers. But I couldn't get the list of kmers of count=0.

You can also use R. I give you an example of only chr1 and only kmer=4.
library(BSgenome.Hsapiens.UCSC.hg38)
library(Biostrings)

genome <- BSgenome.Hsapiens.UCSC.hg38
kmers <- oligonucleotideFrequency(genome$chr1, 4)
kmers

m <- as.matrix(kmers)

m[order(m),]

Grabbing chromosomes for hg38:
$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
$ for fn in `ls *.fa.gz`; do gunzip $fn; done

Via kmer-counter and Python, here's how to search for kmers of length 7 from chromosome chrY:
#!/usr/bin/env python

import sys
import subprocess
import itertools

k = 7
chr = 'chrY'
fastaFile = '%s.fa' % (chr)
kmerCmd = 'kmer-counter --fasta --no-rc --k=%d %s' % (k, fastaFile)

try:
    output = subprocess.check_output(kmerCmd, shell=True)
    result = {}
    for line in output.splitlines():
        (header, counts) = line.strip().split('\t')
        header = header[1:]
        kmers = dict((key,int(val)) for (key,val) in [d.split(':') for d in counts.split(' ')])
        result[header] = kmers
except subprocess.CalledProcessError as error:
    sys.stderr.write("%s\n" % (str(error)))

kmers = result[chr]
comp = {'A':'T', 'C':'G', 'G':'C', 'T':'A'}
for kmerList in itertools.product('ACGT', repeat=k):
    kmerKey = ''.join(kmerList)
    kmerCompKey = ''.join(reversed([comp.get(b,b) for b in kmerList]))
    if kmerKey not in kmers and kmerCompKey not in kmers:
        kmers[kmerKey] = 0

for key, val in sorted(kmers.iteritems(), key=lambda (key,val):(val,key)):
    sys.stdout.write("%s\t%s\n" % (key, val))

This script will print a two-column tab-delimited text file to standard output, where the first column is the 7mer and the second column is the count of that 7mer and its reverse complement over chrY (including zero-counts):
CGACGCG 20
CGTCGCG 20
CGCGATA 23
TACGCGC 25
...
AAATAAA 33521
TTTCTTT 34014
GAATGGA 35361
AATGGAA 36906
TTTTTTT 103093

With a couple tweaks, this script could run over a range of k and chromosomes for your reference genome.

authors: Hášíl Páůďýál


text:
Does the MAPQ=0 fraction of a BAM file depend on the insert sizes?
When doing Illumina 2x150bp sequencing of genomic DNA, and after aligning the reads to GRCh38, does the percentage of the non-N fraction of the human genome as MAPQ=0 depend on the insert sizes of the genomic fragments?
This is, for two identical samples with identical final coverage, sample A having an average insert size of 250bp and sample B having an average of 450bp, would the fraction of MAPQ=0 change between the two?

Yes, as a general rule of thumb mappability increases with insert size (up to a limit) and read length. Whether this will actually occur in a given case will depend more on how randomly the sequencing samples from the genome to begin with (i.e., if the library prep happens to select for/against high mappability regions then the insert size won't much matter). The basic reason behind this is that if one end of a read aligns to a repetitive element that the likelihood of being able to use the other end as an anchor will increase a bit if it's further away. Now having said that, Illumina sequencers in particular like a certain fragment length for efficient sequencing, so you can't go too crazy on increasing the insert sizes and still get decent output (unless you modify your library prep).

This is, for two identical samples with identical final coverage, sample A having an average insert size of 250bp and sample B having an average of 450bp, would the fraction of MAPQ=0 change between the two?

For human, I expect sample B will have fewer MAPQ=0 reads (I don't have concrete data to prove). A large part of human genome consists of ALU repeats. The average size of an ALU is ~320bp. If your insert size is ~450bp, you can map reads into most of isolated ALUs in principle. With ~250bp insert, fragments in the middle of ALUs won't be mappable uniquely.
Note that to see the difference, you need to use a mapper that sufficiently takes the advantage of paired-end information.

Does bwamem take advantage of it?


text:
Filter BAM file for read pairs where one or both of the reads starts with a given sequence pattern
What's the most straightforward way to filter a BAM file for read pairs where one or both of the reads starts with a given sequence pattern? Would a tool that deals with UMIs work for this? Which one would be best?

Are you interested in READs that start with a particular sequence or ALIGNMENTS that start with a particular sequence, because they are not the same thing for 2 reasons:

The sam format records sequences in genome orientation. Thus if the sequence of the read starts AAAA, but the read maps to the '-' strand the BAM seq field will not start AAAA, but end TTTT

Mappers may clip or otherwise only map part of a read. When they do this, they might only record part of the sequence of the read in the SEQ field. Thus a line with the CIGAR 5H90M will not contain the first 5nt of the read seq

First option
My solution uses the following steps:

use picard sortsam to sort the records on query-name (not samtools sort because the order is not the same between java and C )
use jjs (java scripting engine) and the htsjdk library to build a bufferof reads having the same name. If any read starts with a pattern, print the whole buffer. The htsjdk should be already contained in the picard.jar.
sort back to coordinate order using samtools.

All in one:
java -jar /path/to/picard.jar SortSam I=in.bam O=/dev/stdout  SO=queryname |\
jjs -cp /path/to/picard.jar  script.js |\
samtools sort -T TMP -o out.bam -

with script.js
var out=java.lang.System.out;
var File = Java.type("java.io.File");
var SamReaderFactory = Java.type("htsjdk.samtools.SamReaderFactory");
var SamInputResource = Java.type("htsjdk.samtools.SamInputResource");
var SAMFileWriterFactory = Java.type("htsjdk.samtools.SAMFileWriterFactory");
var samReader =  SamReaderFactory.makeDefault().open(SamInputResource.of(java.lang.System.in));
var samWriter =  new SAMFileWriterFactory().makeBAMWriter(samReader.getFileHeader(),true,java.lang.System.out);
var iter=samReader.iterator();
var buffer=[];
for(;;)
        {
        var rec=null;
        if(iter.hasNext()) rec=iter.next();
        if(rec==null || (buffer.length>0 && !rec.getReadName().equals(buffer[0].getReadName())))
            {
            var i=0;
            for(i=0;i< buffer.length;++i)
                {
                if(buffer[i].getReadString().match(/^AAAA/)) break;
                }
            if(i!=buffer.length)
                {
                for(i=0;i< buffer.length ;++i)
                    {
                    samWriter.addAlignment(buffer[i]);
                    }
                }
            if(rec==null) break;
            buffer=[];
            }
        buffer.push(rec);
        }

samReader.close();
samWriter.close();

A second option is to use samtools view and awk to export the read names:
 samtools view  input.bam | awk -F '\t' '{if($10 ~ /^AAAA/) print $1;}' | sort | uniq > names.txt

and to use Picard FilterSamReads and READ_LIST_FILE to filter the reads.
UPDATE: 3rd solution:
I've updated my tool samjdk: there is now a new option --pair to work with query-name-sorted BAM files.
The command would be:
java -jar picard.jar SortSam I=S1.bam O=/dev/stdout  SO=queryname |\
java -jar dist/samjdk.jar --samoutputformat BAM --pair \
     -e 'return records.stream().anyMatch(R->R.getReadString().startsWith("AAAAAAG"));' |\
samtools sort -T TMP -o out.bam -

Does this account for the fact that BAM stores the sequence of a read in genome orientation, rather than read orientation? So if the sequence of a read mapped to the plus strand starts with AAAA it will be recorded in the BAM as ending with TTTT.

@ia no but it would be easy to test this: (R.getReadNegativeStrandFlag()?R.getReadString().endsWith("TTTT"):R.getReadString().startsWith("AAAA"))

I got it to work now with samjdk as explained in 3rd solution. Thanks a lot!

For at least one in a pair matching the sequence, if the read names of pairs are identical, some grep magic could do the job:

Print the sequence at the beginning of the line and read the name in the second column
Find all lines starting with the sequence
Extract only read names

And then just grep all found reads and save into the BAM file:
SEQ=<my_sequence>
samtools view unfiltered.bam | \
   awk '{print $10 "\t" $1}' | \
   grep "^$SEQ" | \
   cut -f 2 > list_of_matching_names
samtools view unfiltered.bam | \
   grep -f list_of_matching_names | \
   samtools view -b > filtered.bam

Two small notes:

You might want to consider also reverse complementary sequences.
You also might want to keep header of BAM file (you can extract it at the beginning and add it in the end).

I am one of the authors of the UMI-tools pacakge, and although you probably could hack together a way to do this with our package, its not what its designed for and the solutions above are probably better.
The best way to do this would be to filter the fastq before mapping if thats a possibility, if it has to be the BAM, then here is a python implementation similar to Pierres java solution above. This one fiddles about with strands and cigar operations a bit more.
First sort the reads by name with samtools sort -n
then
import pysam
inbam = pysam.AlignmentFile("mybam.bam")
outbam = pysam.AlignmentFile("mybam.filtered.bam", "wb", template=inbam)
seq_to_find = "AAGG"
rc_seq_to_find = "CCTT"

def test_and_out_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find):
    '''This function takes a list of reads and check if any of them
    have the prespecific sequences at their start. If the sequence
    is detected the whole buffer is output. The start of the read
    is the end of the alignment if the read is on the '-' strand.
    Also checks that the read doesn't start with hard-clipped sequence'''

    do_output=False

    for segment in read_buffer:

        if segment.is_reverse:
            # the last operation is stored in the first part of the
            #last of the cigar tuples. Operation 5 is Hard-clip
            output_condition = (not segment.cigartuples[-1][0] == 5) and \
              segment.query_sequence.endswith(rc_seq_to_find)
        else:
            output_condition = (not segment.cigartuples[0][0] == 5) and \
              segment.query_sequence.startswith(seq_to_find)

        if output_condition:
            do_output=True
            break

    if do_output:
       for segment in read_buffer:
           outbam.write(segment)

last_name=""
read_buffer=[]

# Loop over reads in the file
for read in inbam.fetch(until_eof=True):

    # if the name hasn't changed, the name is part of the
    # same "pair"
    if read.query_name==last_name:
        buffer.append(read)
        continue

    # if we got here then the a new pair is being examined
    # test the old pair and output if neccesary.
    test_and_output_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find)

    # make a new buffer with just the single new read.
    last_name=read.query_name
    read_buffer=[read]

# the last set of reds will not have been output.
test_and_output_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find)


text:
Index a BAM file using pysam
(How) can you index a BAM file using pysam?
When I tried the intuitive pysam.index I got:
import pysam
my_bam = pysam.AlignmentFile("regular_bwamem_mapping.bam", "rb")
pysam.index(my_bam)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-5-a0b3a04ecb2c> in <module>
----> 1 pysam.index(L_bam)

~/miniconda3/lib/python3.7/site-packages/pysam/utils.py in __call__(self, *args, **kwargs)
     59             args,
     60             catch_stdout=kwargs.get("catch_stdout", True),
---> 61             save_stdout=kwargs.get("save_stdout", None))
     62
     63         if kwargs.get("split_lines", False):

pysam/libcutils.pyx in pysam.libcutils._pysam_dispatch()

TypeError: object of type 'pysam.libcalignmentfile.AlignmentFile' has no len()

disclaimer: when I post a question and answer, it's because I figure the answer while I write the question.

I have tried same command...
error
[E::hts_idx_push] Unsorted positions on sequence #8: 3411049 followed by 3378665
[E::sam_index] Read 'A00404:40:HFGYNDSXX:4:1258:24623:27633' with ref_name='chr8', ref_length=145138636, flags=1107, pos=3378665 cannot be indexed
Traceback (most recent call last): File "Insertx.py", line 283, in  main() File "Insertx.py", line 263, in main pysam.index(str(''.join( data))) File "/usr/local/lib/python3.8/dist-packages/pysam/utils.py", line 69, in call raise SamtoolsError(
pysam.utils.SamtoolsError: 'samtools returned with error 1: stdout=, stderr=samt

@anju there must have been something that gone wrong with sorting the file. Otherwise it would not complain about it.

Oh you silly sausage, pysam.index takes a bam file name, not a python object.
import pysam
pysam.index("regular_bwamem_mapping.bam")

will index your .bam file.

I should print "Oh you silly sausage" and tape it above my monitor.


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

My code is this:
phmw = function(pept_seq , unimod_location, pH ){

  monoisotopic_hydrogen_mass = 1.00794

  if(pH < 7){

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

  }else if(pH > 7){

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

  }else if(pH == 7){

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

  return(mol_weight)

}

Do you think that this approach is correct or not?

The code is correct as far as I can tell but there’s a lot to be said about style. The code can be made a lot shorter and more readable. Furthermore, don’t use T instead of TRUE, and don’t use return` (which is unnecessary here). Generally I suggest reading the Tidyverse style guide, even if you decide not to follow it.

I'm curious as to the practical application? Is your interest specifically in peptides in solution at different pH values? Just wondering because most peptide MW prediction is with reference to mass spectrometry, where pH is not really relevant.

If you want completely accurate monoisotopic mass you should use the mass of the proton rather than hydrogen. (The difference is the mass of electron).
I second the question of practical applicability since each amino acid will have individual pKa, moreover, in the peptide neighboring amino acids will interact. Thus the protonation and deprotonation will happen at the different pH.

Actually I'm trying to mimic the mass spec for some peptides and see if I can characterize them as detectable by their m/z.
@neilfws Didn't understand exactly what you said. The algorithm behind the mw() assigns the masses of the a.a's from mass spec knowledge ? But one can run a mass spec experiment with different pH values ( Something that changes the mass of the total peptide). Which of these masses does the mw() algorithm uses ?

I'm just wondering whether protons gained by protonation due to low pH are relevant to the ions generated during MS. For example in electrospray MS "it became evident that there is no correlation between solution charging and electrospray charging" - https://www.hindawi.com/journals/ijac/2012/282574/. I would think solution protons would be lost during fragmentation and ion generation.

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

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

Amino-acid content (eg Lys, Asp) can have drastic effect on observance in a particular MS-mode.
You should consider including Peptides::pI() or Peptides::charge() to determine the likelihood of (de-)protonation in your calculations.


text:
How do I find identical sequences in a FASTA file?
I want to create a database for a proteomics study. Therefore, the mapping from a given sequence to a protein ID has to be unique. I am wondering whether there is already a built-in function in Biopython for that, but I could not find any. The duplicated sequences should be merged into a single entry and there respective IDs should be just concatenated (see expected outcome below).
Using the following example input
>prot_id_1
MAWIGLDISKLFVENSRDAAA
>prot_id_2
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_5
MAWIGLDISKLFVENSRDAAA
>prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_7
MAWIGLDISKLFVENSRDAAA

I currently do it as follows:
from Bio import SeqIO
import numpy as np
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC

# Read entries into dictionary, gi: sequence
record_dict = SeqIO.to_dict(SeqIO.parse("test.fasta", "fasta"))

# Extract sequence as string
record_dict = {gi: str(record_dict[gi].seq) for gi in record_dict}

# Check for duplicated sequences
seq_id_dict = {}
for k, v in record_dict.items():
    seq_id_dict.setdefault(v, []).append(k)

# Filter sequences that appear more than once
filter_dict = {seqi: li for seqi, li in seq_id_dict.items() if len(li) > 1}

# Merge identifiers to a single one
merged_ids = {"_".join(sorted(li)): seqi for seqi, li in filter_dict.items()}

# Remove all entries that are associated with duplicate IDs
for gi in set(np.concatenate(filter_dict.values())):
    record_dict.pop(gi, None)

# Add the merged IDs again
record_dict.update(merged_ids)

# Create list of Seq objects
final_seq = [SeqRecord(Seq(seqi, IUPAC.protein), id=gi, name='', description='') for gi, seqi in record_dict.items()]

# Write file
SeqIO.write(final_seq, 'test_output.fasta', 'fasta')

which gives me the desired output:
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
>prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS

Is there a smarter way of accomplishing this?

The trick would be to swap the key in the dictionary to be the sequence itself.
Also I would recommend using a different separator that "_" since that is what the current ids have so that you can easily separate the individual ids from the concatenated id. I used a pipe "|" in this example. Also I just manually wrote the FASTA output instead of using Biopython to do it.
from Bio import SeqIO
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)
with open("Output.fasta", 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write(">{}\n".format('|'.join(ids)))
        output.write(seq + "\n")

Yes, that works fine and is more condensed than my solution (your dedup_records is the same as seq_id_dict). But true, all the additional filtering is not needed; Friday afternoon... :)

(I wait with acceptance for a while just to see if there is already some kind of built-in for that)

There is no built in function that I know of that does the id concatenation.  If it did exist it would likely be doing the same logic.

Seems like a good approach if fitting the file into memory is ok, but why not use biopython to write also? This asymmetry seems odd, plus that way you won’t dump all the sequence on a single line

Mainly because I find the BioPython interface cumbersome to write out sequences in a fasta format when the data is not a Seq object already.  @cleb can use Biopython if he wants the main changes are in the use of the defaultdict with the swapped key.

I see, well still consider using record.description not record.id to ensure you get the whole headers. Also you could save memory by storing the sequences as a hash rather than the full sequences

OP referred to the protein ID so thats what I used.  I agree that hashes could be useful to handle large input amounts.

@Chris_Rands: Mind adding the hash solution!? In my case I only have the IDs, no further description but if you could add a more detailed answer I am more than happy to upvote it.

@Cleb I meant to say, the reason I didn't add a hash solution is because it would be tricky to resolve hash collisions (although they'd be very rare)

@Chris_Rands: No worries, the solutions posted are currently sufficient; would have just been a 'nice-to-see' type of thing...

Just as a small variation to @Bioathlete's answer in case you want to write the fasta using Biopython (e.g. to add names and description):
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)

# this creates a generator; if you need a physical list, replace the outer "(", ")" by "[" and "]", respectively
final_seq = (SeqRecord(Seq(seqi, IUPAC.protein), id="|".join(gi), name='', description='') for seqi, gi in dedup_records.items())

# write file
SeqIO.write(final_seq, 'test_output_final.fasta', 'fasta')

This will give the desired output:
>prot_id_1|prot_id_5|prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_2|prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH

I think you mean it creates a generator not an iterator.

@Bioathlete: Of course, fixed, thanks for spotting it.

using datamash ans some awk
$ cat input.fa | paste - -  |\
cut -c2- | sort -t $'\t' -k2,2 |\
datamash  collapse 1 -g 2  | tr "," "_" |\
awk '{printf(">%s\n%s\n",$2,$1);}'

>prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH

This can also be done quite easily with awk:
awk '{if(NR%2){printf("%s\t",$1)}else{print}}' filename | sort -k2,2| awk '{if($2==SEQ){gsub(">","",$1);ID=ID"|"$1} else{if(SEQ!=""){printf("%s\n%s\n", ID,SEQ);}SEQ=$2;ID=$1;} }END{printf("%s\n%s\n", ID,SEQ)}'
The first awk converts the fasta file to a tab separated file with format ID\tSequence, which is then sorted by sequence by sort. The last awk goes through the sorted file looking at the sequences: if the sequence in the current line is the same as that in the previous line, it appends its id to the ID variable. Otherwise (i.e. we have encountered a new sequence), it prints the previous concatenated ids and the corresponding sequence.

Here is my free program on Github  Sequence database curator
(https://github.com/Eslam-Samir-Ragab/Sequence-database-curator)
It is a very fast program and it can deal with:

Nucleotide sequences
Protein sequences

It can work under Operating systems:

Windows
Mac
Linux

It also works for:

Fasta format
Fastq format

Best Regards

Hello Eslam, thanks for your contribution. It would be nice if you could also include how to use your program to achieve what the original poster has asked for... Here we try to be as concrete as possible both when asking and answering :-).

I agree with @KamilSJaron: Please show how to use your package to solve the issue from above (ideally with the same input). That will help others to get started. Thanks!

If you have a very big fasta file(60 GB) why use SeqIO.parse and write the output file using biopython?
Why not use SimpleFastaParser instead (you get a tuple(name, seq)) and you dont need to wait so long as using the parser and using Seqrecords.
What do you think about it guys?
Thanks
Paulo
Kind of:
from Bio.SeqIO.FastaIO import SimpleFastaParser
from collections import defaultdict

dedup_records = defaultdict(list)
for (name, seq) in SimpleFastaParser(filename):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[seq)].append(name)
with open(filename, 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write(">{}\n".format('|'.join(ids)))
        output.write(seq + "\n")


text:
Convert population allele count into population allele frequency by using Perl or Python (translate available Bash into Perl or Python)
I already know the Bash AWK solution and also R code for the question that I am asking. But, my file is so big and with R it takes very long time and I am afraid of AWK's mistake for this computation, so I wonder if there is a Perl or Python solution for my problem.
I have a population allele count data, and it looks like this:
1   0   0   0   0   0   0   0   0   0   1   2   1   0   0   0   0
0   2   0   0   0   0   0   0   0   0   0   4   0   2   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
2   2   0   0   2   1   0   0   0   0   2   4   0   0   0   2   0

Columns are populations and rows are SNPs. I have two rows per SNP (one row for the number of copies of allele "A" in each population and one row for the number of copies of allele "a").
In the example above, the first and second rows are number of alleles for the SNP1, the third and fourth rows are alleles for SNP 2, and like this until SNP20000000.
I want to calculate population allele frequency per SNP for all populations: frequency of allele A at SNP1 in population 1 = number of copies of allele A in population/total number of A/a gene copies in population and frequency of allele a at SNP1 in population 1 = number of copies of allele a in population/total number of A/a gene copies in the population. This means for each SNP number of copies of each allele should be divided to sum of the number of copies of allele "A" and "a" for each population. This is my desired output:
1   0   0   0   0   0   0   0   0   0   1   0.333333    1   0   0   0   0
0   1   0   0   0   0   0   0   0   0   0   0.666667    0   1   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0           0   0   0   0   0
1   1   0   0   1   1   0   0   0   0   1   1           0   0   0   1   0

As I said, I have an R and Bash solution, but is there a way to do this estimation in Perl or Python?
Here is the Bash solution, but I do not know how to translate it into Perl or Python? I am learning Perl...
awk '{for(i = 1; i <= NF; i++)
    tally[i] += $i};
    (NR%2) == 1
    {for(i=1; i<=NF; i++)
        allele1[i] = $i};
        (NR%2) == 0
        {for(i = 1; i <= NF; i++)
        allele2[i] = $i;
            for(i = 1; i <= NF; i++)
                if(tally[i] == 0)
                    tally[i] = 1;
                for(i = 1; i <= NF; i++)
                    printf allele1[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i <= NF; i++)
                    printf allele2[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i <= NF; i++)
                    tally[i]=0}' MyData |
sed 's/\t$//g'

It is strange that with R it takes so much time. Usually with matrices it is quite fast. Maybe you can post the code in R to translate it to python or perl, or improve it. perl is usually better for text editing, as with numbers it is a bit more difficult but it is quite fast usually. What have you tried in those languages? Did you translate your R code? Did you start from zero?

File is 8GB. I included bash solution but do not know exactly how to translate it into perl or python!

If you have a solution that works in on language then I dont' think this is the right forum.  It is no longer a bioinformatics question it is a programming question that happens to use biological data.

It seems like you do 7 loops in your awk code. I am not an expert of awk, but this calls for an optimization, I think it could be reduced to 4 or so. How is your code in R? (I am more confident of my R skills, and perhaps I could help you) And as Bioathlete says this might suit better to code Reviw. Where they could help you speed and improve your code

Could you also share the R solution for for those of us that find awk one liners to be completely impenetrable?

@ LIopis, OK. what is code review? is that what you meant? is that a forum?

@Anna, code review is when someone looks your code to improve it. See this site in the stackexchange network.

Your "bash" code is an awk solution

I don't know how fast you need it to be, but this R solution runs on 20 million SNPs in under 2 min on my laptop:
allele_counts <-matrix(rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0,
                     0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 0, 0,
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     2, 2, 0, 0, 2, 1, 0, 0, 0, 0, 2, 4, 0, 0, 0, 2, 0), 10000000),
                   ncol=17, byrow=TRUE)
allele1<-allele_counts[ c(TRUE,FALSE), ]
allele2<-allele_counts[ c(FALSE,TRUE), ]
total_counts <- allele1+allele2
allele_counts[c(TRUE,FALSE), ] <-allele1/total_counts
allele_counts[c(FALSE,TRUE), ] <-allele2/total_counts

Edit
I decided to a head-to-head test of this solution vs. bli's pandas solution and your awk solution
I used readr for file IO:
library(readr)
allele_counts<-as.matrix(read_tsv("allele_counts.txt", col_names=F))
...
write_tsv(as.data.frame(allele_counts), "allele_freq_R.txt", col_names=F)

Here are the timings:
cat Benchmarks/R.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
103.7556        0:01:43 2211.16 2518.59 2209.11 2209.15 0.00    27.58   73.62

cat Benchmarks/Python.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
132.6425        0:02:12 4872.45 6042.28 4870.89 4870.92 2.67    46.59   99.92
cat Benchmarks/awk.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
126.4455        0:02:06 4.22    413.55  0.93    0.97    0.00    130.68  91.28

The R solution is slightly faster than the pandas solution and uses half as much memory. Your awk one-liner is about the same speed as the pandas script, but requires a tiny fraction of the memory footprint.

Intelligent use of vector reusing! I would recommend scale  though

Thanks. I nicked the idea from here

Thanks for the benchmarks. Maybe you can answer this question of mine, about the meaning of the benchmarking values given by psutil (via snakemake): https://stackoverflow.com/q/46813371/1878788

I'm new new to benchmarking with Snakmake and also couldn't find any useful documentation for it. I did my best to answer with the info I've found, but I guess it's a case of the blind leading the blind

Here is an implementation in pure python. Probably not the fastest, but you should not encounter any memory issues with this.
Note: this is python3, use as follows: cat snp_counts.tsv | python allele_frequencies.py
#!/usr/bin/env python
"""
allele_frequencies.py
"""
import sys

is_odd = True
for line in sys.stdin:
    parts = line.strip().split()
    if not parts:
        continue
    allele_counts = [ int(p) for p in parts ]
    if is_odd:
        allele1 = allele_counts
        is_odd = False
    else:
        allele2 = allele_counts
        allele_sum = [ sum([a1,a2]) for a1,a2 in zip(allele1,allele2) ]
        allele1_freq = [ a / tot if tot else 0.0 for a,tot in zip(allele1, allele_sum) ]
        allele2_freq = [ a / tot if tot else 0.0 for a,tot in zip(allele2, allele_sum) ]
        print(*allele1_freq, sep = '\t')
        print(*allele2_freq, sep = '\t')
        is_odd = True

thanks for the code. I tried your code but I get an error message:  File "./allele_frequencies.py", line 21
print(*allele1_freq, sep = '\t')
^
SyntaxError: invalid syntax.

@Anna1364 Are you sure you are using python3?

In python, if your files can fit in RAM, you might try using pandas and cytoolz.
import pandas as pd
from cytoolz import interleave

def odd(num):
    return num % 2

def even(num):
    return 1 - (num % 2)

# Read in two separate dataframes
alleles_1 = pd.read_table("snps.txt", skiprows=odd, sep="\s+", header=None)
alleles_2 = pd.read_table("snps.txt", skiprows=even, sep="\s+", header=None)
sums = alleles_1 + alleles_2
# https://stackoverflow.com/a/45565745/1878788
pd.DataFrame(list(interleave([
    alleles_1.div(sums).fillna(0).values,
    alleles_2.div(sums).fillna(0).values]))).to_csv(
        "freqs.txt", sep="\t", index=False, header=False)

I don't know how efficient this is. On my workstation it runs in 1m27 on data obtained by taking 1 million times the example given in the question. See heathobrien's answer for benchmarks on larger data.

Interesting, though I actually used the example times 10 million. I added some proper benchmarking to my anwser

@heathobrien Indeed, I had misread your code. With 20 million snps, on my workstation where Firefox is eating one third of the 16GB RAM, my script fails on a MemoryError after 4 minutes.

Here's my quick Perl attempt after misreading the question, which divides each row elements by the row total:
cat popCountData.txt | perl -lane '
  my $total=0;
  grep {$total += $_} @F; ## create per-row total
  foreach (@F) {printf("%0.2f ", ($total) ? $_ / $total : 0)} @F; ## divide elements by total
  print ""; ## print a line break (note: line will end with a space)
'
0.20 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.20 0.40 0.20 0.00 0.00 0.00 0.00
0.00 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.00 0.25 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.13 0.13 0.00 0.00 0.13 0.07 0.00 0.00 0.00 0.00 0.13 0.27 0.00 0.00 0.00 0.13 0.00

The challenge here is that the processing is in chunks of two lines, which means a couple of things:

The previous line (at least for even-numbered lines) needs to be stored
The sum of counts (or whatever statistic is necessary) will have different processing on each line

Unfortunately, I can't quite see how to get from your example input to your example output. For example, how does a count of 2 and 4 turn into a frequency of $1/3$ and $2/3$?

Because in that colum there is a 2+4 = 6 So 2/6 = 1/3 and 4/6 = 2/3, The scale is by each two rows and column wise.

@gringer, thanks but it does not give the output that I want. I added bash solution.Yes, LIopis is right


text:
How do I save selected sequences in seqkit to a file?
Consider:

seqkit grep -nrp ";s__nosocomialis;"  ribogrove_11.217_sequences.fasta

Using this command lets me select all sequences related to the organism "A. nosocomialis" in the RiboGrove database.
How do I save these selected sequences to a new ".fasta" file?

Based on the docs, seqkit grep has a -o "output file" option, so if you get the correct output with:
seqkit grep -nrp ";s__nosocomialis;" ribogrove_11.217_sequences.fasta

Then this command will (in theory) give you your desired outcome:
seqkit grep -nrp ";s__nosocomialis;" ribogrove_11.217_sequences.fasta -o new.fasta

Thanks, it worked. I don't know how i didn't see the output command while surfing the documentation.

For any program that prints output to the terminal (standard output), you can always redirect to a file using > or >>. Either of the two will create the output file if it doesn't exist, but > will overwrite any contents of an already existing file, while >> will append to it without overwriting. So you could just do:
seqkit grep -nrp ";s__nosocomialis;" ribogrove_11.217_sequences.fasta > out.fa


text:
Convert Cotton Probe ID to Gene Symbol
I am new to bioinformatics, my background is in Electrical Engineering. I am trying to convert Affymetrix Cotton Probe IDs to gene symbols. I have a gene expression dataset and I need the expressions for only certain genes. So in the dataset, I have a probe id such as Ghi.10301.1.S1_s_at, how can I get the gene symbol associated with it? Thanks
EDIT: This is the specific dataset I am looking at
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50770

So which file are you trying to parse? Is it GSE50770_series_matrix.txt? The link you gave has various resources associated with it .

One problem you may face is with the state of the cotton genome (Gossypium hirsutum) itself, there are no gene symbols for this genome, you're stuck with accession IDs for each probe, which the provided Affymetrix annotation file has some information on via homology to other characterised proteins and coding sequences, that's about it!  The official genome annotation is not much use either!

The first stop for decoding probe IDs for a microarray platform is usually the GEO Platform Description: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL8672
This description contains a table of annotation, and is provided (in this instance) by the array manufacturer. You can view, or download, the full table: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc=GPL8672&id=1538&db=GeoDb_blob35
GEO Platform annotations do date, and so for up-to-date annotations, download from NetAffx: http://www.affymetrix.com/support/technical/byproduct.affx?product=cotton
The most recent annotations (dated 2016-04-13) are available here: http://www.affymetrix.com/Auth/analysis/downloads/na36/ivt/Cotton.na36.annot.csv.zip.
For the example in your question, looking at this probe in the CSV:
"Ghi.10301.1.S1_s_at","Cotton Genome Array","Gossypium hirsutum","Mar 30, 2016","Consensus sequence","GenBank","Ghi.10301.1","gb:DT050231 /DB_XREF=gi:72275634 /DB_XREF=COT_DW_H07 /TID=Ghi.10301.1 /CNT=21 /FEA=EST /TIER=ConsEnd /STK=1 /UG=Ghi.10301 /UG_TITLE=Transcribed locus, moderately similar to NP_568148.1 ubiquitin conjugating enzyme ubiquitin-like activating enzyme (Arabidopsis thaliana) /REP_ORG=G. hirsutum","DT050231","Ghi.10301","Ghi.10301","---","---","---","---","---","est","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","This probe set was annotated using the Accession mapped clusters based pipeline to a UniGene identifier using 21 transcripts. // false // Accession mapped clusters // E","AI730871,DN803875,DR452973,DR457278,DT050231,DT545642,DT546719,DT546982,DT546983,DT547780,DT547781,DT555433,DT559909,DT564450,DT566924,DT569293,DT569974,DT571307,DT571318,DT571472,DT573548","AI730871 // /DB_XREF=BNLGHi8130 // gb // --- // --- /// DN803875 // /DB_XREF=G.hir-8-10 DAA bolls drought stressed751 // gb // --- // --- /// DR452973 // /DB_XREF=CM008C10 // gb // --- // --- /// DR457278 // /DB_XREF=CM053F01 // gb // --- // --- /// DT050231 // /DB_XREF=COT_DW_H07 // gb // --- // --- /// DT545642 // /DB_XREF=EST1056282 // gb // --- // --- /// DT546719 // /DB_XREF=EST1057359 // gb // --- // --- /// DT546982 // /DB_XREF=EST1057622 // gb // --- // --- /// DT546983 // /DB_XREF=EST1057623 // gb // --- // --- /// DT547780 // /DB_XREF=EST1058420 // gb // --- // --- /// DT547781 // /DB_XREF=EST1058421 // gb // --- // --- /// DT555433 // /DB_XREF=EST1066073 // gb // --- // --- /// DT559909 // /DB_XREF=EST1070549 // gb // --- // --- /// DT564450 // /DB_XREF=EST1075090 // gb // --- // --- /// DT566924 // /DB_XREF=EST1077564 // gb // --- // --- /// DT569293 // /DB_XREF=EST1079933 // gb // --- // --- /// DT569974 // /DB_XREF=EST1080614 // gb // --- // --- /// DT571307 // /DB_XREF=EST1081947 // gb // --- // --- /// DT571318 // /DB_XREF=EST1081958 // gb // --- // --- /// DT571472 // /DB_XREF=EST1082112 // gb // --- // --- /// DT573548 // /DB_XREF=EST1084188 // gb // --- // ---","---"

Your dataset was produced on a microarray. The probe IDs in the FASTA file are similar to read IDs produced from other NGS experiments like RNA-Seq. According to the workflow provided by Ensembl the first step is to align your reads to a reference. You'll only keep the reads that match the reference completely. Once the probes are aligned you'll be able to quantify the number of reads mapping to the gene that the probe was designed to attract.
I did find a file providing the database entries associated with each probe ID (one of which is the ID you provided) for the Affymetrix Cotton Genome Array.

Hi @kohlkopf thanks for your response. I don't need to quantify the reads here, since I already have expression values for different genes in the dataset. What I am trying to do is convert the Probe IDs to gene symbol so I can extract the data for a specific subset of genes I am interested in. I will add a link of the data set I am looking at on the original post.

This is demonstrably not the answer to the OPs question, the OP has not got a FASTA file of probe sequences, equally the hits < 100 step and "quantify the number of reads mapping to the gene" - would only be appropriate for a RNASeq experiment, in this case expression is via probe intensity measurements not read count, counting the number of probes would simply be a bad proxy for gene length.  This is simply an ID mapping issue.

Mapping to Ensembl won't help either as it only has the Gossypium raimondii genome were as the Affy Cotton array is based on Gossypium hirsutum.


text:
Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?
Disclaimer: This is a self answered question for documentation purpose and I adapted this from the following github gist. Especially from users terrycojones and peterjc as well as sujaikumar who raised the issue.

I have a strange situation. I have two BLAST databases (20170821 and 20170824) that both contain a certain subject sequence LC074724. When I query them both with a given query, Q, in one case (20170824) there is a very strong match and in the other (20170821) no match at all.
Why is there not even a lower quality match found although the databases are nearly the same?
BLAST version

$ blastn -version
blastn: 2.6.0+
 Package: blast 2.6.0, build May 11 2017 22:22:40

Both databases contain an identical subject sequence LC074724:
$ blastdbcmd -entry LC074724 -db 20170824 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

$ blastdbcmd -entry LC074724 -db 20170821 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

Only difference for the databases is 20170821 is just 20170824 with an extra 643 sequences. Rest is completely identical.

I'm trying to figure out why. This will be a longer read, Tl;dr at the end:
Doing a match
There is a good match of Q against LC074724 in 20170824:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170824 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
Q   LC074724.1  2581

Note that the bitscore of the LC074724 match is 2581.
But there is no match of LC074724 at all in 20170821:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
-task blastn -max_hsps 1 | grep LC074724

The option -outfmt '6 qaccver saccver bitscore' tells BLAST to print the query accession number and version, the subject accession and version, and the bitscore of the match.
What are the best hits for 20170821
So what are the bitscores of the Q matches in 20170821?
$ blastn -outfmt '6 bitscore' -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | sort -nr | head
3301
2616
2590
2590
2563
2554
2547
2544
2540
2540

The bitscore of Q against LC074724 in 20170824 (2581) would have placed it in 5th position in the top 10 matches of Q, but LC074724 is not matched at all in 20170821! The lowest bitscores matched in 20170821 is 2354.
What's the difference between the databases?
Get the ids of the subjects in each database
$ blastdbcmd -entry all -db 20170824 | egrep '^>' | cut -c2- | sort > 20170824.ids
$ blastdbcmd -entry all -db 20170821 | egrep '^>' | cut -c2- | sort > 20170821.ids

$ wc -l 2017082?.ids
    9754 20170821.ids
    9111 20170824.ids
   18865 total

20170821 has 643 sequences that are not in 20170824:
$ comm -23 20170821.ids 20170824.ids | wc -l
     643

20170824 has no sequences that are not in 20170821:
$ comm -13 20170821.ids 20170824.ids | wc -l
       0

and the two have 9111 sequences in common:
$ comm -12 20170821.ids 20170824.ids | wc -l
    9111

These numbers match with the numbers from wc -l above.
So 20170821 is just 20170824 with an extra 643 sequences.
Maybe LC074724 matches, but with a low bit score?
BLAST only shows the top 500 matched subjects by default. Let's ask it to show more:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 -max_target_seqs 10000 > 20170821-10000-matches.txt

And....
$ head 20170821-10000-matches.txt
Q   Q   3301
Q   X234A   2616
Q   CS388973.1  2590
Q   DM059402.1  2590
Q   LC074724.1  2581
Q   KJ843188.1  2576
Q   AB697496.1  2576
Q   AB697499.1  2576
Q   AB697503.1  2576
Q   AB697508.1  2576

There it is! And it's in 5th position, just as it should be (based on the bitscores discussed above).
Just to confirm, taking away the -max_target_seqs option gets no match:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
$

Three confusingly named BLAST options
I've always found the names and descriptions of the following three BLAST command-line options very confusing.
From the Formatting options section of blastn -help:

-num_descriptions =0>   Number of database sequences to show one-line descriptions for   Not applicable for outfmt > 4
  Default = `500'    * Incompatible with:  max_target_seqs
-num_alignments =0>   Number of database sequences to show alignments for   Default = `250'    * Incompatible with:
  max_target_seqs

And from the Restrict search or results section:

-max_target_seqs =1>   Maximum number of aligned sequences to keep   Not applicable for outfmt <= 4   Default = `500'    *
  Incompatible with:  num_descriptions, num_alignments

From the section names, it seems like only the latter option (-max_target_seqs) would have any effect on the search, and that the former two are just about what is displayed.
But... using either just -num_descriptions or -num_alignments with a big value also gets a match:
$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_descriptions 10000 | grep LC074724
LC074724.1  Hepatitis B virus DNA, complete genome, isolate: p621     2581    0.0
>LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621

$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_alignments 10000 | grep LC074724
LC074724.1  Hepatitis B virus DNA, complete genome, isolate: p621     2581    0.0
>LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621

So BLAST is apparently finding the LC074724 match even without the high -max_target_seqs option, but it's not displaying it unless there is a high -num_alignments or -num_descriptions value. (Keep in mind, this is the 5th best hit! Not somewhere >500)
Just to confirm, when none of the 3 options are given and no -outfmt option is either, the match is not found:
$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | grep LC074724
$

Tl;dr:
In summary, the options
 -evalue
 -max_target_seqs
 -num_alignments
 -num_descriptions

are NOT just filters applied to the results, but instead directly influence the results found.
This means you can End up with a bacterial sequence as BEST hit when using a value x, and with no bacterial sequences at all in your top20 when using value y.
To quote the NCBI user service:

We don't consider this a bug, but I agree that we should document this
  possibility better. This can happen because limits, including max
  target sequences, are applied in an early ungapped phase of the
  algorithm, as well as later.

And

The cutoffs are applied during an ungapped extension phase, and some
  alignments may improve in score during gapped extension. The safest
  approach is to relax (increase) the threshold and deal with a larger
  result set.

This was totally unexpected behaviour for me and is as far as I know is nowhere documented in the BLAST help or similar.


text:
Source for whole genome comparisons at NCBI Genomes
The NCBI Genomes database has these dendrograms for (presumably) whole genome comparisons for certain species, e.g. Pseudonomas aeruginosa  or Escherichia coli.
How were these comparisons done? Someone knows the source / paper?

They say "based on genomic BLAST". No useful information can be found, though.

I inquired into the details of the dendrograms after becoming frustrated with the lack of information. As with Ensembl, I'm sure that the folks at NCBI have a standardized pipeline that they run the sequences through to generate these dendrograms as no specific source deserving attribution was mentioned:

The tree is based on a pairwise, BLAST comparison of chromosome sequences (an assembly is used for genomes with 2 or more chromosomes). BLAST identity is normalized on sequence length. The distance map is generated by a Neighbor Joining algorithm, where the distance is the blast score.
Dr. Wayne Matten, NCBI

Thanks for inquiring. Too bad that their response doesn't contain more details or a link to the pipeline. Seems like they calculate the average PID of all alignments reported by blastn und use that as a distance measure. I wonder if they also consider the parts of each genome pair that are not covered by alignments.