text:
getting KeyError on certain residues in BioPython: ATOM vs. HETATM lines
When I want to load atoms for residuum (' ', 85, ' ') in chain A for protein 3tmm with structure[0][chain_id][residue], but I'm getting KeyError exception:
Traceback (most recent call last):
  [...]
  File "/home/pok/Checkouts/bilkoviny/loader.py", line 118, in _get_residue_atoms
    return list(self.structure[0][chain_id][residue])
  File "/usr/lib64/python2.7/site-packages/Bio/PDB/Chain.py", line 74, in __getitem__
    return Entity.__getitem__(self, id)
  File "/usr/lib64/python2.7/site-packages/Bio/PDB/Entity.py", line 40, in __getitem__
    return self.child_dict[id]
KeyError: (' ', 85, ' ')

I assume this is because this residuum do not have "ATOM" lines in PDB file:
[...]
ATOM    704  CD1 ILE A  84      19.348  -4.839  18.138  1.00 47.08           C
ANISOU  704  CD1 ILE A  84     5862   6108   5917   -216  -1513    609       C
HETATM  705  N   MSE A  85      21.619  -0.130  16.236  1.00 31.89           N
ANISOU  705  N   MSE A  85     5143   3072   3900   -766  -1011    962       N
HETATM  706  CA  MSE A  85      22.277   0.675  15.226  1.00 40.76           C
ANISOU  706  CA  MSE A  85     6708   3869   4911   -972   -729   1162       C
HETATM  707  C   MSE A  85      23.146   1.786  15.818  1.00 44.98           C
ANISOU  707  C   MSE A  85     7191   4141   5757  -1268   -558   1109       C
HETATM  708  O   MSE A  85      24.257   2.045  15.338  1.00 43.69           O
ANISOU  708  O   MSE A  85     7055   3855   5692  -1611   -151   1153       O
HETATM  709  CB  MSE A  85      21.244   1.265  14.269  1.00 46.50           C
ANISOU  709  CB  MSE A  85     8043   4309   5316   -709   -948   1345       C
HETATM  710  CG  MSE A  85      21.849   1.954  13.052  1.00 54.03           C
ANISOU  710  CG  MSE A  85     9604   4929   5995   -907   -623   1638       C
HETATM  711 SE   MSE A  85      20.503   2.755  11.877  1.00 84.35          SE
ANISOU  711 SE   MSE A  85    14348   8336   9366   -486  -1038   1884      SE
HETATM  712  CE  MSE A  85      21.690   3.538  10.518  1.00 83.79           C
ANISOU  712  CE  MSE A  85    15063   7865   8907   -915   -414   2326       C
ATOM    713  N   SER A  86      22.636   2.453  16.846  1.00 43.76           N
ANISOU  713  N   SER A  86     5234   6599   4794  -1778   -946    195       N
[...]

I have been told "MSE is same as MET, it is just with selenium instead of sulphur" and that I probably need to load these atoms (I'm approximating secondary structure element position, so I want individual atom coordinates).
Is there some way how to get atoms for this residue in this protein via BioPython?

You can find the answer in The Biopython Structural Bioinformatics FAQ:

What is a residue id?
This is a bit more complicated, due to the clumsy PDB format. A
residue id is a tuple with three elements:

The hetero-flag: this is 'H_' plus the name of the hetero-residue (e.g. 'H_GLC' in the case of a glucose molecule), or 'W' in the case
of a water molecule.
The sequence identifier in the chain, e.g. 100
The insertion code, e.g. 'A'. The insertion code is sometimes
used to preserve a certain desirable residue numbering scheme. A Ser 80
insertion mutant (inserted e.g. between a Thr 80 and an Asn 81
residue) could e.g. have sequence identifiers and insertion codes as
follows: Thr 80 A, Ser 80 B, Asn 81. In this way the residue numbering
scheme stays in tune with that of the wild type structure.

The id of the above glucose residue would thus be ('H_GLC', 100, 'A').

So in your case it would be ('H_MSE', 85, ' ').

Thank you! Please do you know if there is a way how to get something like (*, 85, *) - I mean, if I do not care about hetero-flag and insertion code at all.

@jhutar ignoring biopython hetero-flag would make sense, as it is biopython-specific workaround for incorrect PDB files; I don't know if it's possible. Do not ignore insertion codes, they are part of the sequence number.


text:
Hba-a2 has two copies on UCSC genome browser
I was browsering the UCSC viewer of the genome at this region.
UCSC genes label two genes as Hba-a2 while GENCODE thinks the left one is Hba-a1
Which annotation is correct?

Gencode is correct (when in doubt, always assume Gencode/Ensembl is correct when it differs from UCSC). What presumably happened is that the UCSC annotation pipeline noticed that Hba-a1 and Hba-a2 are nearly identical (they produce the same protein and have minimal nucleotide difference) and just misannotated.

As a general rule, it is usually safe to believe Gencode over RefSeq and RefSeq over UCSC (not because UCSC is bad, but because it isn't manually curated).

source: Stack Exchange

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

sort: votes


text:
How to show only my read on IGV?
I have a large BAM file on my IGV. I'd like to visualize a read alignment. I have the read name, and I don't care anything else in the file.
IGV has an option "Select by name" (see screenshot). The problem is that IGV doesn't hide the alignments not matching the name. I'd still have to manually navigate on IGV with my mouse to find my read.
That was too cumbersome. What I'm  doing now is make a new alignment BAM file. For example,
samtools view -h L.V.51.bam | grep '@\|HWI-D00119:228:CB08RANXX:8:2116:4095:45403' | samtools view -bS - | samtools sort - > new.bam
I'd load the new alignment file on my IGV.
If I have 20 reads, I'd have to create 20 new BAM files... Is there anything I can do on IGV that I'm not aware of?

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

Filter for reads matching the regex passed to -i option:
[h] for help: grep -i '72789|57688|70785|56299|27116|12086'

(Optionally, show sam records as strings)
[h] for help: print

Example result:

For command help type grep -h and for more filtering options see also awk -h.
Hope this helps and feel free to post questions and open issues on github!

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


text:
Error in .checkedCall : subset indices out of range
After normalization the way it is described here:
https://bioconductor.org/packages/3.7/bioc/vignettes/scran/inst/doc/scran.html
Using quickCluster() function and running computeSumFactors() I am trying to run calculateQCMetrics() the way it is described here:
https://f1000research.com/articles/5-2122/v2
sce <- calculateQCMetrics(sce, feature_controls=list(Mt=is.mito))

And I am getting an error:

Error in .checkedCall(cxx_colsum_subset, exprs_mat, is_feature_control) :
    subset indices out of range

Does anybody know how to correct it? What is going wrong here?

Just posted the question and found out why the error was happening. I had two datasets: raw, and filtered one with half of the rows as raw. I did normalization on the filtered one. Then I calculated is.mito using the unfiltered, raw dataset. So, it produced indexes that were outside of the range of the dataset.


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:
Improve scRNA-seq dataset for further analysis
I got a dataset from C.Elegans scRNA-seq paper:
GSM2599701_Gene.count.matrix.celegans.cell.Rdata in GSE98561_RAW.tar
The dataset is 40 000 x 68 000, where rows represent genes and columns - cells. So, I took it and tried to process myself to build an scRNA-seq pipeline. Here is what I did:

I filtered out the genes that have zero counts in all of the cells and the dataset was reduced to 29 000 x 68 000
I removed all of the cells with counts < 100 in all of the genes - the dataset became 29 000 x 66 000
Then, because the dataset was too big to run normalization even on the cluster with 120Gb RAM (because there are multiple distinct types of cells, first clustering needs to be done), I selected just even columns and ran normalization with 29 000 x 33 000 dataset (UMI_count):
library(scran)
library(scater)
sce <- newSCESet(countData=UMI_count)
clusters <- quickCluster(sce)
sce <- computeSumFactors(sce, clusters=clusters, positive=TRUE)

After running the code above I decided to check whether the data is fine and so I ran:
> summary(sizeFactors(sce))
Min.    1st Qu.  Median  Mean   3rd Qu.    Max.
0.0000  0.0000  0.0000  0.0717  0.0000 33.3900

I also ran PCA on the normalized dataset and it looks like that:

It seems to me that the normalized dataset is terrible and I need to do some further processing before I could run further analysis. What else could I do to improve it? How to filter it? There are no spike-ins and maybe 200 mitochondrial genes. The approach described here does not work, probably because the majority of the cells have low number of genes expressed:

I tried removing low-abundance genes after normalization, but it seems that most of them are going to be removed:
>ave.counts <- rowMeans(counts(sce))
>keep <- ave.counts >= 1
>sum(keep)
109

Should I filter all of the cells - columns - for up to 500 total gene expression count instead of 100? Is it a good idea? I can't think of anything else.

Did you look up at the unnormalized data (via PCA to check if the normalization didn't handle these dataset correctly)? What kind of cells are there? Are they from several cell lines, stages, or disease?

If you look at the samples section of that GSE, you can find that there are "c.elegans cells (L2 stage)", "HEK293T, HeLa S3 and NIH/3T3 cell mixture ", and "frozen fixed HEK293T and NIH/3T3 cell mix", so you might need to subset

What happens if you drop the cells in the bottom 25 or 50% of reads sequenced (and maybe a few of the top ones too)? I suspect that anything that tightens the range of size factors is going to lead to cleaner results.

Several cell lines, several dozen maybe. That is why need to do clustering first. These are all of the cells in C.Elegans. If you look into the raw compressed file, there you can find the dataset (see the top of the question). Before normalization PCA looked the following way: https://bioinformatics.stackexchange.com/questions/3115/filter-out-pca-outliers-automatically        It definitely became better after normalization, if you look at the variance % of the PC1 and PC2, but not effective still

@DevonRyan I am not sure what you mean regarding size factors. Could you please elaborate on this?

You used the sizeFactors() function in your post, I mean its output.

I usually use a minimum of 200 genes/cell for relaxed filters (500 if possible).
For genes I use presence in at least 10% of the cells.
A few other filters you can try are:

High percentage of reads from only top50 genes. This means you haven't captured a lot of genes and the diversity of the cell won't be useful.
High percentage of mitochondrial genes. This has been reported to be linked to bad quality cells. The idea is that mRNA is depleted because of early cell death (or other outer membrane injury) leading to having only mito mRNA left.
Minimum UMI per cells. Similar to minimum genes but more dependent of the protocol used.
UMI per gene/cell. This can be interesting to detect a potential bias from amplification and is also dependent on the protocol.

The best would be to compare all your samples using violin plots with those different filters and choose your thresholds from there.
A quick example (using Seurat):
Pre filtering:

Post filtering:

I like that top 50 fraction metric. Any idea why it's so clearly bimodal? Does it correlate with percent mitochondrial?

Most of them yes. I haven't looked at them specifically but I'm pretty sure a few other cells are also filtered with this one.

I don't think you can conclude that the dataset is terrible based on that PCA. Depending on the specific protocol, each scRNA-seq dataset is going to be very different. Unlike bulk RNA-seq where all the samples are going to be of very similar quality, individual cells will be highly variable. For example, the most basic QC metric is the number of reads. If some of your bulk RNA-seq samples have 10X more reads than other samples, that's a cause for concern. In scRNA-seq experiments, that's not at all surprising. The PCA exposes the major sources of variation. In scRNA-seq experiment, those are probably going to be technical rather than biological, so the output may not make much sense. Furthermore, even if all the cells are of perfect quality, PCA is not really optimal for sparse data with many dropouts like you see with scRNA-seq.
Regardless of PCA, those library sizes seem very small. Based on your plot, there are barely any cells with over 10K reads. I understand this is C elegans, but that seems somewhat low. That may be why you are having trouble with normalization. It would help with understanding the data if you changed the X-axis on those plots so that you can actually see where the peak is. Right now, it looks like 0, but that can't be the case.

Given that this is single-cell sequencing results, you should be using rowSums instead of rowMeans for filtering genes. Most genes won't be expressed at a single cell level, but there should be a lot of genes that have some expression. The zero counts are probably better represented as "missing" rather than zero.
Unfortunately that feature makes running PCAs difficult, because PCAs don't work with missing values. The missing values may need to be imputed in order to generate a better PCA.
disclaimer: I don't [typically] work with single cell data, so don't have a good handle on the methods for looking at it.


text:
Approximate matches with bedtools intersect
I am using bedtools to test whether a set of variant calls produced by a new method is missing any variants from a set of known validated variants.
bedtools intersect -v -a newcalls.vcf -b goldstandard.vcf

It's working as intended for SNV calls, but it occasionally misses indels due to a single base pair discrepancy. It's usually an alignment artifact that doesn't affect the interpretation of the variant.
My first thought was to make a BED file with intervals that extend a few bp beyond each variant call, and comparing this BED file to the known variants. However, I was wondering whether approximate matching is possible with bedtools intersect?

You can use bedtools intersect but first use bcftools norm.  bcftools norm will left-align and normalize indels to get both vcf files consistent with each other.  This should resolve the matching if they are truly matching variants.

This could resolve some of the cases if the input files are not normalized. However, there could be indels that are normalized and 1 bp away from each other (i.e. close but not touching each other). In my understanding, the OP wants to consider these cases as genuine matches.

I will let the OP clarify but it seemed to me that these were matching variants just with alignment artifacts.  If that is true then they should normalize out.

You could use bedtools closest with -d option (-d: In addition to the closest feature in B, report its distance to A as an extra column.). Indels that have no match in file A would be those where the closest indel in file B is more than x bases away.

You have to be careful here that you are not matching a snp with an indel or vise versa.


text:
GSEA- parsing error for .cls file
I am creating a .cls file for my GSEA analysis and see the below error. Can anyone suggest a solution?
I create the file using TeXworks. The file looks like this:
37 2 1
#CD69neg CD69pos
CD69neg CD69pos
CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos

Error Details
---- Full Error Message ----
There were errors: ERROR(S) #:1 Parsing trouble edu.mit.broad.genome.parsers.Par ...

---- Stack Trace ----

of exceptions: 1
------Number of items found in cls data 10 is not equal to the number specified on the header line 37------

edu.mit.broad.genome.parsers.ParserException:
Number of items found in cls data 10 is not equal to the number specified on the header line 37 at
edu.mit.broad.genome.parsers.ClsParser.doSanityChecks(ClsParser.java:388) at
edu.mit.broad.genome.parsers.ClsParser._parse_genecluster_style_categorical(ClsParser.java:252) at
edu.mit.broad.genome.parsers.ClsParser.parse(ClsParser.java:226) at
edu.mit.broad.genome.parsers.ParserFactory._readTemplates(ParserFactory.java:342) at
edu.mit.broad.genome.parsers.ParserFactory.readTemplate(ParserFactory.java:293) at
edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:753) at
edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:726) at
edu.mit.broad.genome.parsers.ParserWorker.doInBackground(ParserWorker.java:52) at
javax.swing.SwingWorker$1.call(SwingWorker.java:295) at
java.util.concurrent.FutureTask.run(FutureTask.java:266) at
javax.swing.SwingWorker.run(SwingWorker.java:334) at
java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1149) at
java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:624) at java.lang.Thread.run(Thread.java:748)

You seem to have forgotten the space after the #. It should be as in here (Just for future readers know where to find the specs). Also it should be a single line (Maybe this error was introduced while pasting here or it has been a duplicate because there are 39 samples in the 3 line and they should be 37).
Either adding 2 more samples to the initial count and removing the newline of the third line:
39 2 1
# CD69neg CD69pos
CD69neg CD69pos CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos

Or removing the third line:
37 2 1
# CD69neg CD69pos
CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos

Thank you so much! Yeah the error was about I have multiple columns rather than just one line.

Glad I could help. Don't forget to mark an answer as accepted or post yourself the solution you find most useful


text:
Parse out exon coordinates from bed file for each gene
This is my bed file for all the exon coordinates I want to take out all the exons of a given gene. Let say I have gene in chr 1 which starts from chr1 11868   12227 so I want to parse out all the exons that comes in between 11868    12227.
This is my small subset:
cat exon.bed | head -10
chr1    11868   12227   +   exon
chr1    11871   12227   +   exon
chr1    11873   12227   +   exon
chr1    12009   12057   +   exon
chr1    12178   12227   +   exon
chr1    12594   12721   +   exon
chr1    12612   12697   +   exon
chr1    12612   12721   +   exon
chr1    12612   12721   +   exon
chr1    12974   13052   +   exon

How do I parse out?
I use mostly R and bit of shell script but I'm not sure if I can use R, may a few lines of perl or shell script can help me solve my problem.

cross posted: https://www.biostars.org/p/292247/

yes i posted earlier in biostar so i didn't have any response so i posted here hoping to get some solution

There are multiple ways to go about it. On the command line you can make a 1 line BED file:
chr1   11868   12227

And then bedtools intersect with it.
In R, you could load your original BED file and use GenomicRanges:
library("GenomicRanges")
bed = read.delim("foo.bed", header=F)  # Rename this
# N.B., BED files use 0-based coordinates, I've switched to 1-based here
exons = GRanges(seqnames=levels(bed$V1),
                ranges=IRanges(start=bed$V2+1, end=bed$V3),
                strand=bed$V4)
# Again, I switched to 1-based coordinates...you don't have to do that
region = GRanges(seqnames="chr1",
                 ranges=IRanges(start=11869, end=12227),
                 strand="+")
theExonsYouWant = subsetByOverlaps(exons, region, type="within")

You'll need to either specify a strand in your region or use ignore.strand=T in subsetByOverlaps(). I hope you don't have any genes that overlap on the same strand, since they'll prove problematic.

thanks for the reply so what I want is I want to parse out the exon coordinates then  11868   12227   +   exon
chr1    11871   12227   +   exon
chr1    11873   12227   +   exon
chr1    12009   12057   +   exon
chr1    12178   12227   +   exon as i you can see 11868 to 12227 there are like 5 exons so i want only the first and the last exon to be taken and then i want to use it for divergent primer..so the way you suggested so i have to kind of define the exon i want to parse every time, is there a better way, because , i want to do the same for all the genes with their exons...

You can use findOverlaps() with something and itself.

okay will try your suggestion

I get this error "  invalid class “GRanges” object: 'mcols(x)' is not parallel to 'x'

That should probably be its own question, though my guess is that you have something the wrong length somewhere.

Via BEDOPS bedops -n and Unix I/O streams:
$ echo -e "chr1\t11868\t12227" | bedops -n 1 exon.bed - > answer.bed

Or, if you have your genes in a BED file called genes.bed:
$ bedops -n 1 exon.bed genes.bed > answer.bed

If you have your genes in some other format, like GFF or GTF, you can use gff2bed or gtf2bed, e.g.:
$ bedops -n 1 exon.bed <(gff2bed < genes.gff) > answer.bed

Or:
$ bedops -n 1 exon.bed <(gtf2bed < genes.gtf) > answer.bed

The file answer.bed will contain exons that do not overlap a gene annotation.


text:
how to plot LD between SNPs and CNVs
I have a bunch of CNVs which I have estimated LD between each allele of CNV with all SNPs across the entire human genome. I used my R-script to do a pairwise LD analysis.
The following is the first few lines of my output for the first alternative allele of CNV and some SNPs on chromosome 2 (just as an example)
CNV_chrom  CNV_allele  SNP_chrom  SNP_position  R2
"1" "  CNV2"   "2"        " 1000"    "-0.0463767871013283"
"1" "  CNV2"   "2"        " 1029"    "-0.201478578143067"
"1" "  CNV2"   "2"        " 1044"    "-0.595441932439136"
"1" "  CNV2"   "2"        " 1045"    "-0.60363963711489"
"1" "  CNV2"   "2"        " 1130"    "0.377970958849221"

I want to plot the LD results (LD between each allele of CNV and all SNPs) ,is there any way that I could do it in R?
I have not done LD plotting before so just slightly confused.
I saw for instance LDheatmap package in R, but do not know how can I summarize the results for thousands of SNPs. Thanks for any suggestion

Summary of LDHeatmap

Produces a graphical display, as a heat map, of measures of pairwise
linkage disequilibria between single nucleotide polymorphisms (SNPs).
Users may optionally include the physical locations or genetic map
distances of each SNP on the plot.

Citation
Shin, J.-H., Blay, S., McNeney, B., & Graham, J. (2006). LDheatmap: An R Function for Graphical Display of Pairwise Linkage Disequilibria Between Single Nucleotide Polymorphisms. Journal of Statistical Software, Code Snippets, 16(3), 1–9. https://doi.org/10.18637/jss.v016.c03

There are a lot of ways to plot LD results.  Here are two links which might help.  http://csg.sph.umich.edu/abecasis/class/666.03.pdf  https://www.ndsu.edu/pubweb/~mcclean/plsc731/Linkage%20Disequilibrium%20-%20Association%20Mapping%20in%20Plants-lecture-overheads.pdf They have a lot of different examples.

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


text:
What is the ICGC normalized_read_count?
I downloaded gene expression data (exp_seq) from the ICGC file browser.
For each sample and gene, the file contains a normalized_read_count.
What is that value? I couldn't find any information on the ICGC website. The values are definitly too low for TPM.

By reading this thread on seqanswers and by comparing the data to TCGA, I figured out

raw_read_count is the read count which you use as input for e.g. DESeq2. It has been estimated using RSEM
normalized_read_count is equivalent to the scaled_estimate from TCGA. This is the estimated fraction of transcripts made up by a given gene, as estimated by RSEM. Multiplying this value with 1e6 yields the TPM.

Glad you could find it yourself! Many thanks for posting the answer

Thank you for your answer. This is an essential information for a user of this data. I wonder why this information on the ICGC website itself though. I found a page describing the columns on github https://github.com/icgc-dcc/dcc-docs/blob/master/docs/dictionary/release-20/sequencing-based-gene-expression-expseq-primary-file-p.md. It does not explicitly say what normalized_read_count are. Also the page corresponds to one of the older release (20), not the current one (28).


text:
Can second search with limited search space improve result in MS proteomics data?
First, please know that I am a novice in this field and that it is thus very likely that you will regard this question as basic or flat stupid. Please be gentle and don't be afraid to criticize. I am asking this so I can learn. Also, don't feel afraid to discuss in comments if you can't give a straight answer; all input is useful.
I am currently working in a project where we have fractionated samples in order to dig deeper into the proteome. The fractionated samples have been analysed with shotgun LC-MS and then searched with a contemporary software using a whole proteome database and to the reversed data base in order to estimate FDR. The search result is cut-off so that the highest scoring sequences are included with an estimated FDR of 1%. In this project we have 150 samples in which we can expect to find certain proteins in almost all samples, while others are unlikely to be found in any of them (but we don't know which yet). If we were to decrease the search space to include only relevant proteins (which as previously stated are not known to us yet), we would boost our rate of true positives included after FDR cut-off, as we decrease the search space which noise randomly might match to with a high score (right?). As we can't know which ones are relevant beforehand, would it be possible to determine relevant proteins with a preliminary search? That is, to do as follows:

Search against whole proteome with a certain FDR cut-off (possibly larger than 1%)
Make a new data base including only proteins that has been found in any of the samples.
Search data against new data base.

To me it seems like there is a possibility that this could increase the number of true positives that are found, but I also feel like I might be missing something crucial. Please, enlighten me: Is this methododology unbiased and valid? What risks may there be? Anything more one has to think about before attempting something like this?
The LC-MS was a top 10 shotgun MSMS using a quadrupole-Orbitrap instrument. The search was done using MaxQuant where fragmentation spectra are use to identify MS1 peaks, which intensity is registered as a quantitative datum. All files were searched together but are configured to be separated in the results file in regards to both sample and fraction. We will likely not perform additional analysis of these samples.

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


text:
How to download the whole BLAST nt database into a specific folder?
I have successfully downloaded the whole nt BLAST database into the current folder using:
wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz"

However, I would like to download it in a specific folder, e.g. output/
I have tried:
wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz"  -O output/

but it does not work. Nothing happens. My 'wget-log' is empty.
I could move to that folder and download it there but I would like to remain in the same working directory.
C.

-O is to specify the name of the output file. You want the -P option.
So, from your example
wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz" -P output/

I recommend using NCBI's tool for this called update_blastdb.pl, shipped with the blast suite (I use the conda environment from bioconda).
cd to the directory you want nt DB in and run:
update_blastdb.pl --decompress nt

Pros: if a download fails for whatever reason it is retried, downloads are md5sum checked and re-running in a folder with an existing DB will (as the name suggests) only update the relevant DB portions. More information here!


text:
Is it possible to find out the PacBio chemistry from the PacBio Sequel BAM's header?
Is it possible to find out the PacBio chemistry from the PacBio Sequel BAM's header? E.g. whether it was produced with P5-C3 or P6-C4? I checked mine but I haven't found it but maybe it's hidden in one of the tags whose meaning I don't know yet, like SEQUENCINGKIT=100-861-800?

According to PacBio Sequel BAM format manual, you can find information on the chemistry version in the header entries @RG (read group), tag DS (description), entry BASECALLERVERSION .
In case of future very PacBio-specific questions, you can contact Liz Tseng on Twitter (@magdoll), one of the most responsive PacBio developers. They also have a dedicated SeqAnswers forum, which is directly monitored by them and you are very likely to get a fast and accurate answer.

I have 5.0.0.6236 for BASECALLERVERSION. Does it mean chemistry P5-C3 was used?

This surely means chemistry is higher than P4-C2 / 2.0, which is not supported in pbbam. The logical answer would be then yes, P5-C3 chemistry (or higher) was used.


text:
Intron retention events for a condition sample with 3 replicates
I have 3 biological replicates of RNA seq data for a particular condition. I want to find out intron retention events from those biological replicates for a given condition. There is no comparison I am doing here. Just from the three replicates, I am trying to find out what introns are retained in my sample. There is only one condition here. Tools like rMATS require 2 conditions. Is there any tool I can use to find out intron retention events given that I only have one condition (with replicates).
Ideally I would like to use a method similar to that used in this paper, though I can't fully grasp their methodology. This would allow me to get p-values for enrichment of specific introns.

All introns will have some level of retention and I guess no intron will be 100% retained. In the absence of a comparator, you will need to work out how much intron retention an intron needs to be called retained. This is a human judgement, there is no automated way of deciding this.
What you need to do this is some sort of PSI calculation (percent-spliced-in). MISO will provide this, but it doesn't do replicates, so I guess you'll need to run it on each sample and then average the PSIs and determine if it surpasses some sort of threshold you have set.

I second this, but would suggest using SUPPA which handles replicates out of the box and is actively developed. As for deciding which introns are retained, I expect that most biological conditions most of them will have a low PSI (which a nice, easily interpretable metric), so I have in the past plotted the distribution and used some upper bound of the distribution to be on the safe side. That said, some condition to use as a reference would be better.

The method you linked to is interesting, though I suspect that Ian's answer is both easier and generally sufficient to get what you want.
Having said that, the general idea behind the method you linked to is to use DESeq2 to test for differentially expressed introns. To generate the control samples, they're effectively taking the number of reads across all introns in a gene and distributing them over the introns. They weight the distribution such that longer introns receive more reads (they use an effective length for this, but I'm skeptical that this makes that much of a difference). DESeq2 is then giving the p-values you're after. Unless you can get the authors to send you there code, you'll have to write a bit of code to do all of this...but it shouldn't be too bad.
BTW, you can use the GenomicFeatures package from bioconductor, specifically the exonsBy(..., by="gene) follows by gaps().

Hi Devon,
Thanks for your help.The paper compares insilco replicate counts with actual counts and that is what is input for DESeq2. From your explanation, for a gene, they took all the reads which map introns and then lets say we get 1000 reads for all the possible introns for this gene and then for each intron, they give a weightage like 1000/(len of intron). I am looking at GenomicFeatures and I have already created the intron file coordinates which you mentioned at seqanswers. If possible and if you have time can you guide me in writing the above code. Thanks for your help always
Regards

It's not 1000/(len of intron), but rather 1000/(len of intron / sum of all intron lengths in gene).

I only have time to provide vague guidance on the code. You can use findOverlaps() (from GenomicRanges in bioconductor) for the counts. Aside from that it's mostly split() and mendoapply() that you'll need.

so that would give us insilco replicate1 count for a particular intron and similarly insilco replicate2 count for 2nd replicate and then we compare these 2 with actual counts we are getting in rep1 and rep2 via DESeq2

Yup, exactly. Once you have the counts its just down to using DESeq2 as normal.

D'oh! I've actaully done this before (that'll teach me not to look at the links in the question). I used feature counts to count the reads in introns, then a python script to calculated the effective length from mappability tracks downloaded from UCSC, and then used an R script to find the "detained introns".  See these gists: https://gist.github.com/IanSudbery/df92ebf4f84d76027394750a6516db1b , https://gist.github.com/IanSudbery/eb36ba5d59b4455f9cd7b2085447eb3d. I offer no guarantees with these, they were pulled from my archives.

It doesn't work too bad. Better on nuclear RNAseq than cytoplasmic, but it is very sensitive. Finds introns with only a small amount of retention.

IRFinder can be used to detect the IR events and can easily reach your goal.
Here is an example. IRFinder can calculate the IR ratio to detect IR event. mimirna.centenary.org.au/irfinder/example1.html. (You just need run IRFinder for one condition and each replicate will output a folder). Two conditions can be used to detect the differential IR events. If you have only one condition, you can select the IR events from the common events of the three replicates or select IR events exist at least two of three replicates based on the main output files( named IRFinder-IR-dir.txt). The IR events can be selected by filtering the IR ratio, Intron coverage etc

Hi Renhua, this answer doesn't seem to really address the issue. Could you please elaborate how IRFinder helps when there is one condition with replicates? How would that be done with IRFinder? Many thanks for your time


text:
Assign cell types to groups of cells based on their gene expression profiles
I have large filtered, normalized dataset of scRNA-seq data of C.Elegans species. Rows are genes (10 000), columns are cells (66 000). Let's say that I got 40 different groups of the cells based on their expression profiles, how could I now assign cell types to these groups? I would expect to get some reference dataset for C.Elegans and compare it with what I got, but I honestly have no idea of where to get such a dataset, and whether there are any R packages that do that. Any advice would be greatly appreciated.

I think the analysis will involve two steps:
1) Compile a reference list of cell types, with a list of genes expressed in each of them, based on the literature. C. elegans has only ~1000 cells, so I'm sure it's been well researched. Mark which genes are just expressed in that particular cell type, and which are uniquely expressed, this latter will be more useful markers.
Resources to do that:

http://openworm.org/science.html#other
http://www.wormatlas.org/
http://www.wormbase.org
Pubmed.
Your lab's knowledge.

This will involve mostly manual work, but maybe you can find some list.
For example, platelets express VWF. And activated platelets express SELP.
2) Use one of the single-cell RNA-seq packages (e.g. Seurat) to cluster cells and to find the "marker" (DE) genes for each cluster (vs. all other clusters).
For each cluster, look up the marker genes in the reference list and infer the cell type.
You may need to infer the cell type from the functions and descriptions of the genes.
Also see Retrieve detailed gene descriptions
Keep in mind that you will have many confusing doublets in the dataset expressing characteristics (genes) of two cells.
You may also need to re-do t-SNE and clustering on a subset of cells, e.g. blood/immune cells, in an iterative, hierarchical way.

Do you have information on which cells are from the same worm? Once you have a reference list, you can try and identify cells of just one worm.

This is an old post, but we have recently published a web-based tool called CIPR (cluster identity predictor) that enables users to upload their data along with a custom reference file for annotating unknown clusters. You can probably obtain relevant reference datasets from GEO.
You can read the manuscript in BMC Bioinformatics for more information about how the algorithm works.


text:
Is this a low complexity region in our human genome?
I have a screenshot where many of my reads were aligned to a region that I suspect a low complexity region. Although you can't see, all those reads are clipped in the cigar strings. Sample cigar strings: 64S32M29S, 74S32M18S etc... Consequently, the actual sequence mapped to the genome is less than the read length.
I have a feeling that my alignments are bogus because of the complexity and clipping, but I'm not sure if the complexity is indeed low. All I see is a few "T" between a bunch of "A". How to define a low complexity region?
Q: Is this a low complexity region in the human hg38 genome? Would that be fair if I report "the alignments are likely be errors due to the low complexity and alignment clipping"?

Yes, that's a short low complexity region wedged between a SINE and an snRNA. More importantly, your alignments have MAPQ of 0 (that's why they're filled with white in IGV), which will happen if they map equally well to multiple locations. Without looking at the sequence one can use that alone to determine that these are not trustworthy mappings.

Oh, I didn’t know about white meaning MAPQ=0. Very useful!

Yes, this is a low complexity region.
Regions are considered low complexity (or having a simple sequence) when they contain an abundance of a single base, or an abundance of short tandem repeats. The simplest is a tandem repeat of a single base (e.g. AAAAAAAAAAAAAAAA), also called a homopolymer. In your case there is a 3-unit tandem repeat of TAAAAA, and surrounding As (with a few intermingled Ts). Depending on the genetic history of this region, it might be that those other As and Ts are derived from the same TAAAAA unit, with some slipping that has occurred during replication to create a cryptic (i.e. non-perfect) repeat structure.


text:
Find number of possibilities
We have a column which has following format : (1, 1, 1, 0, 1).
Here is the definition of what means two columns are compatible.

Using the notation Oi to denote the collection of rows possessing a 1
  in their i-th element, we conclude that one of three possibilities
  must be true: Oi ⊆ Oj , Oj ⊆ Oi , or Oi and Oj are disjoint (the first
  two cases include the possibility that Oi = Oj). If columns i and j
  satisfy this condition, then we call them compatible

The question is : How many possible columns of length 5 are compatible with the column (1, 1, 1, 0, 1)?
You can read more here where problem is described in page 14-16.
Thanks in advance !

I'm voting to close this question as off-topic because it is not a question about bioinformatics.

@AlexReynolds, I attached a link and you will see that is about bioinformatics.

In the document there are also algorithms for finding genes tree.

That the topic it can be applied is in bioinformatics doesn't make your question about bioinformatics. You can edit explaining what is 1 and 0 and which is your bioinformatic question. But it seems a nice problem

Let $n$ be the number of rows (in your case $n=5$) and $|O_j|$ the cardinality of $O_i$.
We can calculate the number of compatible columns as sum of the number of $O_i \subseteq O_j$, $O_j \subseteq O_i$ and $O_i, O_j$ disjoint.
number of $O_i \subseteq O_j$:
Here we need to find the number of supersets of our $O_i$. That means we need to find the number of sets which contain the indices in $O_i$ and for each index not in $O_i$, we have two choices: It is contained, or it is not contained in $O_j$. Therefore the number of supersets is $2^{n-|O_i|}$, which is 2 to the power of elements not in our set.
number of $O_i, O_j$ disjoint:
Here all elements of $O_i$ are fixed not to be in $O_j$ and for each index not in $O_i$, it can again be in $O_j$ or not. That means that the number of possible sets is again $2^{n-|O_i|}$.
number of $O_j \subseteq O_i$:
Here we need to find the number of subsets of $O_i$ and this is the number of elements in its power set, which is $2^{|O_i|}$.
So finally, this means that the number of compatible columns can be calculated as: $2^{n-|O_i|} + 2^{n-|O_i|} + 2^{|O_i|} - 2$.
Why -2? We have counted 2 sets twice:

The empty set, which is a subset of $O_i$ and also counted in the number of disjoint sets
and $O_i$, which we counted in $O_i \subseteq O_j$ and in $O_j \subseteq O_i$

This means for your examle, we have $16 + 2 + 2 - 2 = 18$ compatible sets.
Note: The all zero vector (which corresponds to $O_i=\{\}$) is a special case for which the formula above does not work, as all disjunct sets and all supersets are counted twice. For this special case, the number of compatible columns is simply $2^n$.

Thanks for your answer. How do you calculate |Oi| ?

$|O_i|$ is the number of elements in $O_i$, so for your column $|O_i|$ is 4 because $O_i = {1,2,3,5}$

In my example why is not 5 ? Cardinal({1,1,1,0,1})=5

The cardinality of your column is 5. But in your definition: "Using the notation Oi to denote the collection of rows possessing a 1 in their i-th element,", $O_i$ is the set of rows with a 1. In your column there are 4 rows with a one, so the cardinality of $O_i$ = 4

Thanks a lot for your explanation.Now I understand. It means that for (0,0,0,0,0) the answer is 2^5+2^5+2^0-2 =63 ?

The number of compatible columns can not be larger than the total number of sets with 5 elements ($2^5$) :). (0,0,0,0,0) Is indeed a special case for which my rule above breaks, as all disjoint sets and all supersets are counted twice. So for This special case you have to hard code the result to $2^n$. For all other cases, the formula should work.


text:
ethnicity check either from bam or vcf files
What tool could I use to check the ethnicity of a human bam or vcf file? I would like to use the results as a QC check to know whether a given sample or set of samples match the ethnicity information annotated in the meta data or not.
Seen so far:

Similarity:

The examples work, but starting from a GRCh38 .vcf file (from GATK UG) gives me an error:
Read    mytest.vcf   :1 Individual(s)Exception in thread "main" java.lang.NumberFormatException: For input string: ""
        at java.lang.NumberFormatException.forInputString(NumberFormatException.java:65)
        at java.lang.Integer.parseInt(Integer.java:592)
        at java.lang.Integer.valueOf(Integer.java:766)
        at Individual.get_Individual(Individual.java:73)
        at similarity.main(similarity.java:80)

NGSrelate:

It won't work with GATK UG .vcf files straightaway, it needs to reprocess .bam files using the ANGSD pipeline:

VerifyBAMID:

This tool does not give ethnicity checks but only genotype checks. Running all-vs-all vcfs vs bams can create a semi-matrix of similarities.

The main difficulty here is the use of GRCh38.  Unfortunately, despite the fact that it's more than four years old, the major ethnicity-labeled public datasets (1000 Genomes, gnomAD when allele frequencies are enough) still aren't available for that reference.  It is necessary to perform a liftover operation, or just use overlapping rsIDs and hope for the best.
Let's suppose you go with overlapping rsIDs, and a list of those rsIDs, one per line is in 'rsids.txt'.  (Due to strand flips between reference builds, you may want to restrict this list so that all allele codes match, and it's also reasonable to throw out A/T and C/G SNPs.)  Then, the following process would work:

Download ADMIXTURE (https://www.genetics.ucla.edu/software/admixture/download.html ), plink (https://www.cog-genomics.org/plink/1.9/ ), and/or chr1-chr22 of 1000 Genomes phase 3 (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ ) if you don't already have them.
Use plink to extract the overlapping higher-MAF rsIDs from each 1000 Genomes phase 3 VCF ("plink --vcf ... --extract rsids.txt --maf 0.05 --make-bed --out ..."), and then merge the resulting per-chromosome filesets ("plink --merge-list ... --out merged_phase3_subset").
ADMIXTURE prefers a dataset with around 100k variants.  Unless you have a small list of overlapping rsIDs (in which case you had better use another method), you should have more than that.  plink's LD-pruning function is a good way to select a subset for ADMIXTURE's use:
"plink --bfile merged_phase3_subset --indep-pairwise 500kb 1 0.2; plink --bfile merged_phase3_subset --extract plink.prune.in --make-bed --out admixture_data"

Adjust the 0.2 threshold as necessary to keep the right number of variants.

Run ADMIXTURE in unsupervised mode ("admixture admixture_data.bed 5 -j8"; adjust the -j parameter depending on the number of processor cores).  This generates an admixture_data.5.P file with population allele frequencies, and an admixture_data.5.Q file with sample population assignments.  Verify that the .Q file corresponds to the 1000 Genomes phase 3 'superpopulations'.
Convert your data to plink-format if necessary, keeping just the overlapping rsIDs and 'downgrading' to GRCh37 coordinates.  Make sure your variants are sorted in GRCh37-coordinate order, and the allele order also matches that in your reference dataset ("--a2-allele admixture_data.bim 6 2" during your final --make-bed operation will do the trick); otherwise ADMIXTURE won't do the right thing.  Then run ADMIXTURE in projection mode ("cp admixture_data.5.P my_data.5.P.in; admixture -P my_data.bed 5").  my_data.5.Q will then have the ethnicity estimates you're looking for.


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:
Pathifier bioconductor package recommends to use technical noise for an argument. How to determine it?
I'm running the pathifier approach against C2 pathway curated
database for a specific microarray dataset.
As I was reading the documentation of the pathifier, in order to configure it properly for my dataset, I saw that for the min_std argument they suggest to use the technical noise.

min_std:   The minimal allowed standard deviation of each gene. Genes
with lower standard deviation are divided by min_std instead of their
actual standard deviation. (Recommended: set min_std to be the
technical noise).

But how am I going to calculate it ? At first run I just used the value 0.2254005 from the example they provide with that package, but this might not be the right for my dataset.
Update:
So from here it says: The minimal standard error (min_std) was set as the first quartile of the standard deviation in the data. How am I supposed to calculate this first quartile of the std? Any idea? Should I calculate the sd() for each gene and then find the first quartile?
What do you think about this function?

Answer from @llrs, converted from comment:
Probably is a question for the maintainer but I would guess that is a variance of the gene in the control dataset. So, if in normal samples the gene has a standard deviation of 0.5 this is the expected too in the tumoral cells. calculate the sd only on the control samples. So if the healthy/control samples have such variation the altered ones should have these or more according to the article (if I understand it correctly).
std <- apply(data[, controls], 1, sd, na.rm = TRUE); quantile(std, 0.25) should be much faster and efficient
But  I would check it in support.bioconductor with the authors/maintainers.


text:
What is the meaning of the Hmmer profile?
I want to understand what the hmmer Profile used in hmmer.
consider the profile
HMMER2.0  [converted from 3.1b2]
NAME  aRNH_profile
LENG  121
ALPH  Amino
RF    no
CS    no
MAP   yes
NSEQ  34
DATE  Tue Apr 26 21:24:57 2016
XT       -8455     -4  -1000  -1000  -8455     -4  -8455     -4
NULT        -4  -8455
NULE       656  -1722     98    419   -333    475  -1125    239    250    947  -1073   -271    -50   -338    114    451    113    430  -2131   -717
HMM        A      C      D      E      F      G      H      I      K      L      M      N      P      Q      R      S      T      V      W      Y
         m->m   m->i   m->d   i->m   i->i   d->m   d->d   b->m   m->e
           -6      *  -7932
     1    -214   -916   -565    694    862  -1448    512   -838    115  -1042    170   1348  -1481    369   -158   -290   -159   -625    -52   1675     1
     -    -210   -336    220    -37   -342    377     73   -669    210   -511   -732    404    423     70     27    439    205   -414   -161   -176
     -     -18  -6890  -7932   -894  -1115   -701  -1378    -18      *
     2     285    990  -2213  -1678   2305  -1819   -698    219  -1289   -140    545  -1167  -1835   -973  -1322   -424   1025    762   -559    377     2
     -    -210   -336    220    -37   -342    377     73   -669    210   -511   -732    404    423     70     27    439    205   -414   -161   -176
     -     -18  -6890  -7932   -894  -1115   -701  -1378      *      *
[....]

what exactly does each component mean? By definition a HMM $\lambda$ consist of an alphabet $V$ in my case twenty amino acids. Different states $S$ in our case main, insertion and deletion (beginning and end). And two probability matrices of state transitions and emission in certain state.
now if I look at the profile I can see that it has the allowed transitions
m->m   m->i   m->d   i->m   i->i   d->m   d->d   b->m   m->e
and all the amino acids, but what does each row in the table mean. I checked the manual from the documentation and it says rather ironically

If you were to look at this file (and you don’t have to – it’s
  intended for HMMER’s consumption, not yours), you’d see [part of my output]

are these number translations of the probabilities in $[0,1]$  into another interval $[-x,+y]$? what is the meaning of the three rows under A 1
-124
-210
 -18

Is this the transition matrix? It would be great if someone could explain this  or point to the correct reference.

The tutorial in the documentation does indeed state what you quoted but a full description of the file format is given in Section 8 of the documentation, starting on page 106.
To summarise, as you noted, the number columns A-K is the alphabet size (the number of amino acids). The number of rows underneath is the number of nodes, which corresponds to the number of match states. The three numbers (-214, -210 and -18 in your example) are: the match emission, insert emission, and state transition values.
Caveat: The HMMER profile HMM format changes frequently between HMMER versions (even between minor versions) so there is no guarantee this information is correct for all the older file formats.

Quoting the relevant part of the documentation on page 109:

main model section
All the remaining fields are mandatory. The first
two lines in the main model section are atypical.‡ They contain
information for the core model’s BEGIN node. This is stored as model
node 0, and match state 0 is treated as the BEGIN state. The begin
state is mute, so there are no match emission probabilities. The first
line is the insert 0 emissions. The second line contains the
transitions from the begin state and insert state 0. These seven
numbers are: B → M1, B → I0, B → D1; I0 → M1, I0 → I0; then a 0.0 and
a ’’, because by convention, nonexistent transitions from the
nonexistent delete state 0 are set to log 1 = 0 and log 0 = −∞ = ‘’.
The remainder of the model has three lines per node, for M nodes
(where M is the number of match states, as given by the LENG line).
These three lines are (K is the alphabet size in residues):
Match emission line
The first field is the node number (1 . . . M). The
parser verifies this number as a consistency check (it expects the
nodes to come in order). The next K numbers for match emissions, one
per symbol, in alphabetic order. The next field is the MAP annotation
for this node. If MAP was yes in the header, then this is an integer,
representing the alignment column index for this match state
(1..alen); otherwise, this field is ‘-’. The next field is the CONS
consensus residue for this node. If CONS was yes in the header, then
this is a single character, representing the consensus residue
annotation for this match state; otherwise, this field is ‘-’. The
next field is the RF annotation for this node. If RF was yes in the
header, then this is a single character, representing the reference
annotation for this match state; otherwise, this field is ‘-’. The
next field is the MM mask value for this node. If MM was yes in the
header, then this is a single ’m’ character, indicating that the
position was identified as a masked position during model
construction; otherwise, this field is ‘-’. The next field is the CS
annotation for this node. If CS was yes, then this is a single
character, representing the consensus structure at this match state;
otherwise this field is ‘-’.
Insert emission line
The K fields on this
line are the insert emission scores, one per symbol, in alphabetic
order.
State transition line
The seven fields on this line are the
transitions for node k, in the order shown by the transition header
line: Mk → Mk+1, Ik, Dk+1; Ik → Mk+1, Ik; Dk → Mk+1, Dk+1. For
transitions from the final node M, match state M + 1 is interpreted as
the END state E, and there is no delete state M + 1; therefore the
final Mk → Dk+1 and Dk → Dk+1 transitions are always * (zero
probability), and the final Dk → Mk+1 transition is always 0.0
(probability 1.0).

Thanks, Chris. I missed the later chapters. Sorry that was too quick of me.


text:
Where can I find gene expression data on one of the cell lines in NCI-60?
I'm trying to retrieve the most expressed genes from a cell line from the NCI-60, lets say cell line NCI-H23, but I'm not sure where to find the gene expression data. My main goal is to create a network using that data and iRefIndex to be loaded into Cytoscape. Does anyone know where I can find the correct data or guide me in the right direction?

CellMiner portal has plenty of data about NCI-60 cell lines, inclusive microarray gene-expression.
Also, NCI-60 is a subset of the ~1000 cell lines covered by the GDSC project for which you have multiple types of measurements (e.g. copy-number, mutation, gene-expression).

Cosmic Cell Lines Project provides gene expression data for 1000+ cell lines (including NCI-60 panel together with NCI-H23 cell line). For each gene in a cell line, you can get Z-score and over/under-expression status calculated across all tumor samples diploid for that gene.


text:
Where did the BIND website go?
I was trying to reach the BIND (Biomolecular Interaction Network Database) but the website is down.
Where can I find something similar? Has this platform moved somewhere or what?

Answer from @heathobrien converted from comment:
looks like that link has been dead for years:
https://www.biostars.org/p/16763/
Answer from @jason-aller converted from comment:
The biostars article mentions this:
http://baderlab.org/BINDTranslation


text:
How to demultiplex pair-end fastq reads with barcode 2 in the identifier line?
I have multiplexed pair-end fastq reads with dual barcodes. The issue is that one barcode is present in the header and one is present at the beginning of the read. I need a method to demultiplex this data, but in order to assign a read to an individual, both barcodes are required, as there is overlap between the barcodes. It seems there are packages available to demultiplex using header ID or in-line barcodes to demultiplex, but not both.
example reads:
@700819F:525:HT235BCXX:2:1101:1139:2144 1:N:0:ATCACGAT
CGAATTGCAGATTTTTTCTGAATAAAGCAGTGCAATAAAATTCCCCGCAAAAACACTTNANNNGNNNNNNNNNNNNNNNNNNNNANNNNNNGNTAATAAA
+
GGGGGIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#<###.####################.######<#.<<GGGI
@700819F:525:HT235BCXX:2:1101:1212:2172 1:N:0:ATCACGAT
AAGGATGCAGGGCATCTCCCTCAGGCTGCGCTCTATCGAAGTCATCCCAGAATTAGATTCCGACCACAGACCAGTCTTAGTCAAACTAGGACCCGAGTGT
+
GGGGGIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIGIIIIIIIIIIIIIIIIIIIIGIIIIIIGGIIIGIIIIIIIIIIIIIGGIGIGIIIGIGI
@700819F:525:HT235BCXX:2:1101:1110:2173 1:N:0:ATCACGAT
GGTTGTGCAGAAAGAGTTGCTGATAAACTTAGCCATGCAGAACAGAATTATGAGTTAGAAGTATGTATATATATACCAATCACTATATCAACCCATTACC
+
<G.G<GGIIIG.GA.GAGAAGG<AA<A<.<<<GA<.<<.G<.G<<A<GGAA....<G.G..<<.GA..<A.<GG<<<.<..<<.GG.A..G..<<<.<<G

One of the barcodes is the sequence found at the end of the header line (ATCACGAT). The second barcode is found at the beginning of the sequence (CGAAT in read 1). Both are required for assignment to an individual.

Can you define where in the read the barcode should be found?

One of the barcodes is the sequence found at the end of the header line (ATCACGAT). The second barcode is found at the beginning of the sequence (CGAAT in read 1). Both are required for assignment to an individual. This is ddRADSeq data, if that is important.

Ok that make sense.  I didn't know the expected length of the barcode in the read.

You could try something like that:
#!/bin/sh

fq_in="${1}"

bioawk -c fastx '{print substr($comment, 7, 15)}' ${fq_in} | uniq \
    > "barcode1.txt"
bioawk -c fastx '{print substr($seq, 1, 5)}' ${fq_in} | uniq \
    > "barcode2.txt"

for bc1 in $(cat "barcode1.txt")
do
    for bc2 in $(cat "barcode2.txt")
    do
        fq_out="${bc1}_${bc2}.fq"
        bioawk -c fastx -v bc1=${bc1} -v bc2=${bc2} \
            'substr($comment, 7, 15) == bc1 && substr($seq, 1, 5) == bc2 {print "@"$name" "$comment"\n"$seq"\n+\n"$qual}' \
            ${fq_in} \
            > ${fq_out}
    done
done

This might be slow. You can a little bit of time if you directly provide lists of barcodes for the for loops.
It could be more efficient to write a custom program that can write to multiple files. It happens that I wrote something that does this in python just yesterday, but only reading the inline barcode. Modifying it to also use the in-header barcode is likely simple.

Since you've found separate programs to demultiplex by Illumina barcode and in-line barcode (I presume this is actually just the sequence from the cut site), just run one and then the other on its output. Typically the Illumina barcode would have been used to demultiplex the files before you even received them (if not, your sequencing facility is incompetent), in which case you just need to demultiplex by the in-line barcode

The program extract from UMI-Tools will extract the barcode from the read and add it to the ID. It will also handle pair-end data so that even if the barcode is only on one of the reads, it will be added to the ID of both reads.


text:
What is the detection limit of protocol in scRNA-seq
In Shapiro et al., when discussing about loss of molecules as source of error in single-cell sequencing, it is written that:

Another source of error is losses, which can be severe. The detection limit of published protocols is $5$–$10$ molecules of mRNA.

I used to think that the detection limit of a protocol was the average maximum number of molecules detected across cells in a scRNA-seq experiment. However, by working on the dataset provided by The Tabula Muris Consortium et al. (pp) I noticed that RNA raw counts often spikes to more than 10 molecules.
What is the actual meaning of detection limit of a protocol in scRNA-seq?

The detection limit refers to the ability to detect a given gene/transcript. In this case, they're saying that a given transcript requires 5-10 copies in order to be likely detected (in other words, the loss rate is 80-90%, which is directly related to a previous question of yours).

Exactly, in that question I misunderstood the concept of detection limit. Now it's clear, thanks.


text:
Plot percentage of genome covered
Given a aligned bam file (wgs/bwa-mem) what steps do i need to perform to generate a plot as seen below (from this paper):

We are trying to see how much of genome was covered using pacbio and illumina for cow.

You mean remove this left plot?

The general steps are:

Install deepTools
Run plotCoverage on your BAM files
Remove the left-most plot in photoshop/gimp/imagemagick

updated my question

Yes, I mean the left plot in the example, since it looks like you want % >= the coverage value on the X axis.

Devon gives a practical answer using the great deepTools package. These should be the theoretical, lower level steps behind the plots:

Compute read depth at each genomic position, using for example samtools mpileup with appropriate flags to accept/reject reads. You get a table like:
Pos       Depth
chr1:1    0
chr1:2    10
chr1:3    11
...
chrX:100  10

For each unique value x in the Depth column count how many genomic positions have depth >= x. You can use Unix tools like sort ... | uniq -c | awk ... for this. You get a table like:
Depth   N bases with depth >= x   % bases with depth >= x
0       (all genomic positions)   100 %
1       1000                      90 %
2       995                       85 %
.
.
.
N       2                         0.1 %

Now just plot column Depth (x) vs % bases with depth >= x (y) with R or you favourite plotting tool.

Did I get it right?

yes, not sure bust does picar collectwgsmetrics do something similar?


text:
Gene expression data in GC-RMA format?
I'm having a hard time retrieving the most expressed genes in a cell line. One of the databases I've found has values for a lot of genes in 60 different cell lines and it seems the data is obtained using GC-RMA. The data is available here in the file NCI60_U133A_20070815.raw.csv.
My problem is what does the values actually represent? Would it be correct to assume the genes with high values are the ones most expressed in the corresponding cell line?
I've found a short explanation on why there is no unit definition on the page here.
Here is a snapshot of the file:
-            huh-7    huh-7      HEK293  ...
1007_s_at    31.727    55.551    360.731
1053_at     402.633   441.622    557.782
117_at        9.968    10.178      9.565
...

These values come from a microarray so they represent the number of sequences that where attached to the microarray's probes (1007_s_at, ...).
In general yes, a higher value means more expression. The normalizations (GC-RMA) applied are to correct for several effects. They help for testing the probability of the given probe to be a false positive.


text:
Counting reads within the intron
I have paired-end RNA-seq samples and a list of intron coordinates in bed format. I want to count the intronic reads such that:
1. The reads should overlap the intron by atleast 25 bp
2. If both the reads of the pair are within the intron, count it as 1.

What tool can I use to compute it.
Thanks

How are your reads stored? Are they BAM or are they BED. If they are BED how is the pairing encoded?

Hi,
My reads are in bam file, I can convert them into bed file. I am trying to follow what you did with your written scripts. Although I am having issues following how you calculated effective length under null model. I am trying to first make a non overlapping set of introns just like DEXseq manual says for exons

You could use BEDOPS bedmap to map reads to introns, using 1) the --count operator to do counting of reads overlapping by your criteria; and, 2) the --indicator operator to do a true/false operation, where reads are contained entirely within the intron.
For instance, to count reads that overlap introns by at least 25 bases, use --count to count and --bp-over 25 to specify the minimum overlap threshold:
$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed > answer.counts.bed

To filter for introns where reads are contained entirely within the intron, use --indicator to get introns that meet the criteria of full containment of reads, via --fraction-map 1:
$ bedmap --indicator --echo --fraction-map 1 --delim '\t' introns.bed reads.bed | awk '($1==1)' | cut -f2- > answer.indicator.bed

Once you have answer.counts.bed and answer.indicator.bed, you can compare those two results with bedmap --indicator --exact to see where there are commonalities (true/1) and differences (false/0), and increment the final count where there are identical introns, e.g.:
$ bedmap --indicator --echo --exact --delim '\t' answer.counts.bed answer.indicator.bed | awk '{ $NF += $1; print $0; }' | cut -f2- > answer.final.bed

If your reads are in a different format, i.e. BAM, you can use bam2bed in a process substitution, replacing reads.bed with <(bam2bed < reads.bam), e.g.:
$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed <(bam2bed < reads.bam) > answer.counts.bed

Etc.
You can also just convert separately, so that you only need to do it once:
$ bam2bed < reads.bam > reads.bed
$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed > answer.counts.bed

Etc.

Thanks Alex,
Can you explain what is --indicator and --fraction-map doing in generating answer.indicator.bed. I don't get this part?

I suggest reading the documentation, as it goes into much more detail: http://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html

You can almost do what you want using featureCounts from the subread package. The only problem is that featureCounts requires GTF rather than bed.
If you had a GTF of your introns, you could do:
featureCounts -a introns.gtf my_bam.bam -f -p --minReadOverlap=25 -o counts.tsv

The problem with this then, is getting the GTF when you currently have a BED.
One way to get the GTF would be to use CGAT.
cgat gtf2gtf -I ensembl.gtf.gz --method=merge-exons -S genes.gtf.gz
cgat gtf2gtf -I genes.gtf.gz --method=exons2introns -S introns.gtf

or you could put this together into one pipeline:
cgat gtf2gtf -I ensembl.gtf.gz --method=merge-exons  | cgat gtf2gtf --method=exons2introns -S introns.gtf

This will first merge all the exons of a gene together to make a single, collapsed transcript. It will that invert this to find that parts that are in the introns of this collapsed transcript. What you end up with is parts of the gene that are never exonic in any of the transcripts.
alternatively, you could convert your bed regions into GTF:
cgat bed2gff -I introns.bed -S introns.gtf --as-gtf


text:
Read counts from BAM file
I have few BAM files which are generated from Ion Torrent Server (ampliseq) aligned to hg19 genome. I want to extract read counts from the bam files and I know that "featureCounts" can be used for this. But before getting the read counts, to view the bam file I gave the following command:
samtools view -H IonXpress_input.bam

Output:
@HD VN:1.4  SO:coordinate   GO:none
@SQ SN:NR_039978    LN:984
@SQ SN:NM_130786    LN:1766
@SQ SN:NM_138932    LN:9293
@SQ SN:NM_033110    LN:2678
@SQ SN:NM_000014    LN:4678
@SQ SN:NM_144670    LN:5192
@SQ SN:NR_040112    LN:1201
@SQ SN:NM_017436    LN:2106
@SQ SN:NM_016161    LN:1771
@SQ SN:NM_015665    LN:1854

I guess these are from header. After these I see something like following:
@RG ID:JODXM.IonXpress_022.1    CN:S5TorrentServerVM/MT205524   DT:2017-12-07T12:04:52+0100 FO:TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGA KS:TCAGTTCGAGACGCGAT    PG:tmap PL:IONTORRENT   PU:s5/540/IonXpress_022 SM:SNU449_TEAD4

Haven't seen a bam file like this before. Any idea what are these and how I can get read counts?

Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR_* or NM_* accession codes with your samtools view -H code. With samtools idxstats file.bam you'll get the reads per chromosome (or per nucleotide sequence in your case).
For use in featureCounts my advice is to get the fastq files, and align it yourself to hg19 (or hg38 to be more up to date).
If your data is properly aligned to hg19 your output should be something like this:
samtools view -H aligned2hg19.bam

@SQ SN:chrM LN:16571
@SQ SN:chr1 LN:249250621
@SQ SN:chr2 LN:243199373
@SQ SN:chr3 LN:198022430
@SQ SN:chr4 LN:191154276
@SQ SN:chr5 LN:180915260
@SQ SN:chr6 LN:171115067
@SQ SN:chr7 LN:159138663
@SQ SN:chr8 LN:146364022
@SQ SN:chr9 LN:141213431
@SQ SN:chr10    LN:135534747
@SQ SN:chr11    LN:135006516
@SQ SN:chr12    LN:133851895
@SQ SN:chr13    LN:115169878
@SQ SN:chr14    LN:107349540
@SQ SN:chr15    LN:102531392
@SQ SN:chr16    LN:90354753
@SQ SN:chr17    LN:81195210
@SQ SN:chr18    LN:78077248
@SQ SN:chr19    LN:59128983
@SQ SN:chr20    LN:63025520
@SQ SN:chr21    LN:48129895
@SQ SN:chr22    LN:51304566
@SQ SN:chrX LN:155270560
@SQ SN:chrY LN:59373566

So, I confirmed with the people who gave me these bam files. These are aligned to hg19 genome only.

Please show your output to these people.

But you mean the bam file I have doesn't have any chr, read count and other details? I was thinking the above posted one are from header and there is another command to get the readcounts and other details.

It probably has read counts, but not per chromosome but per @SQ sequence (which are RNA sequences). With samtools idxstats you'll see the read counts per @SQ.

Ok. I want to use featurecounts in R to get the read counts. I have the bam files and hg19.fasta and hg19.fasta.fai files. May I know how I should give the command to get the counts for all the samples. Thanq

Maybe you misunderstood, but you need to correctly align to hg19 first in order to do this in featureCounts. From your output it is clear that it is not aligned to hg19, but to a bunch of RNA sequences. Please discuss with your people that your data is not properly aligned to hg19 (show your output, this is proof), or try to align it yourself from fastq. When you have done that, you'll need a proper gtf annotation file from UCSC or ENSEMBL which is necessary for featureCounts. Good luck!

Your options are to either align to the genome (as @b.nota said) or filter using samtools however you like (e.g., for MAPQ) and then get counts with samtools idxstats.

When you are told they are aligned to hg19, I guess they mean tat are aligned to the hg19 version of the transcriptome, not hg19 itself. In this case you done need featureCounts to get read counts (and indeed, it would be hard/impossible to do so). samtools idxstats should give you want you are looking for.

Yes, thanks a lot for all. And yes it is aligned to hg19 transcriptome. And will give a try with idxstats.

One question. Do you I need to sort the bam files? samtools idxstats what are the options I should give in the command to get read counts?

Yes you'll need to sort bam files first. samtools idxstats aln.sorted.bam will give you the results, see here.

Yes, I saw the manual and tried like this and endup with an error.
samtools sort -T /tmp/input.sorted -o input.sorted.bam input.bam. This gave me [W::bam_hdr_read] EOF marker is absent. The input is probably truncated; [E::bgzf_read] bgzf_read_block error -1 after 0 of 4 bytes;
[bam_sort_core] truncated file. Aborting.

The tmp dir has 6 input.sorted.nnnn.bam files.

Please start a new question for that error.


text:
Duplicate gene symbol handling in GEO gene expression data
I have downloaded a gene expression data from GEO database (GSE3268) in which in some of its rows there are duplicate gene symbols.
For example TP53 exists in two rows with different expression values and different GenBank Accession IDs.
They have the following Nucleotide Title:

Human p53 cellular tumor antigen mRNA, complete cds
Homo sapiens tumor protein p53 (TP53), transcript variant 1, mRNA

I want to map these data with edge list of gene symbols to make a network. Which expression value should I consider for such duplicate gene symbols for mapping?

Since this is a microarray (U133A), what you're seeing is that each probe is associated with one or more transcripts. The general strategy when dealing with microarrays is to use RMA (or fRMA) to summarize to probesets and then collapseRows to summarize to gene-level information. Since this is all typically done in R, you can find some further discussion on the bioconductor site. I would suggest going through this process and then only using the result of that in building your network (whether you're happy with sticking to transcript-level quantifications or gene-level quantifications for building your network is up to you).

I have a problem with using fRMA since the there is no supplementary .cel file in GEO database for this dataset. How can I use R on .SOFT files to summarize gene information?  Another point is that in one of the subsets of my data set, (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM73388), it seems that RMA have been used and data is collected in .SOFT file, but again there is duplicates in it.

If RMA has already been used then you needn't run it a second time.

Very useful answer! In EuropePMC I found 63 publications mentioning the use of collapseRows(). Many  don't report the setting they apply. Some high impact publications reported as follows: "MaxMean" (default) = Gandal (Science2018), Raznahan (PNAS2018), Langfelder (NatNeurosci2016); "maxRowVariance" = Langfelder (NatNeurosci2016)


text:
sorting BAM file error using samtools
I have few bam files and would like to get read counts using
samtools idxstats

[Data is aligned to hg19 transcriptome]. To use that command I need a sorted bam file. So to sort them I gave the following command.
samtools sort -T /tmp/input.sorted -o input.sorted.bam input.bam

This ended up showing:
[W::bam_hdr_read] EOF marker is absent. The input is probably truncated.
[E::bgzf_read] bgzf_read_block error -1 after 0 of 4 bytes
[bam_sort_core] truncated file. Aborting.

And the tmp dir has 6 input.sorted.nnnn.bam files.

The cause of the error is that your input file is truncated or /tmp is running out of space. If you can do samtools view -H input.bam without error (reading the header also checks for the magic number at the end of the file) then it's like the that /tmp is filling up.

I did "samtools view -H input.bam" there is no error. It ended up showing all the header things which I asked here [https://bioinformatics.stackexchange.com/questions/3080/read-counts-from-bam-file/3081?noredirect=1#comment4889_3081]

Then either /tmp is full or otherwise one of the files there is somehow corrupt (you can test them with samtools view -H too).

Ok. looks like that file is corrupt. And I tried with other bam file and showed - [bam_sort_core] merging from 24 files... and I see a sorted.bam file. I used "samtools idxstats sorted.bam" - This is what I see - [bam_idxstats] fail to load the index.

Don't forget to samtools index sorted.bam

Ok. So, after the indexing I used idxstats and it showed 4 columns without any headers.                                                                                NR_039978 984 38 0
NM_130786 1766 1 0
NM_138932 9293 37 0
NM_033110 2678 2 0
NM_000014 4678 4 0
NM_144670 5192 18 0
NR_040112 1201 2 0
NM_017436 2106 76 0
NM_016161 1771 0 0

They are (in order): contig name, contig length, number of mapped reads, number of unmapped reads (this will typically be 0).

Ok. Thankyou. So, for DE analysis what I need is mapped reads. how do I get the gene names? Can I use edgeR for DE analysis using the mapped reads?

You mapped to the transcriptome, so they're the first column.

Yes, I see that. Is there any way to get the gene id/ transcript id from this results. And may I know the answer for DE analysis edgeR?

If you want to get annotation (like gene symbols) coupled to your accession numbers, please read the biomaRt user guide. If you have questions after reading it, please start a new question in stead of discussing it in the comments.


text:
Import a tab-separated file with differing numbers of elements in each row; prokka output
I am using prokka to annotate a bacterial genome:
prokka ecoli.fa
Prokka is outputting a tab-separated file (called PROKKA_12142017.tsv) with differing numbers of elements in each row:
locus_tag    ftype    gene    EC_number    product
BOFHCHGE_00001    CDS    hypothetical protein
BOFHCHGE_00002    CDS    rdgB_1    3.6.1.9    dITP/XTP pyrophosphatase
BOFHCHGE_00003    CDS    hypothetical protein
BOFHCHGE_00004    CDS    hypothetical protein
BOFHCHGE_00005    CDS    hypothetical protein

Does anybody know how I can get round this to import into a dataframe in R?
Is there a way to ask prokka to add something in the gene column?

Please make sure to send a bug report to the prokka authors, since this should be fixed on their end.

Author here - this was a bug which was fixed in August:
See the issue in Github.
You will need to use the latest Github HEAD version until the next release.

You don't really need to get around anything, your example can be loaded correctly with read.delim(). The rows with missing values are filled in with blanks. You may, however, prefer the readr package, which handles this a bit more elegantly (it'll tell you which values were missing and not do the annoying string->factor conversion).

Thanks, useful but it puts 'hypothetical protein' under the gene column when it actually belongs in the product column. The gene column should be empty for hypothetical proteins.

Oh, I thought that was intentional. In that case prokka has a bug and you'll need to test for "" as the product and just assign accordingly.

And I get 'Sec-independent protein translocase protein TatA' for EC number, so that's definitely wrong. I never see an empty element before a non-empty element.

You can use the stringsAsFactors = FALSE to prevent converting strings to factors in read.delim/read.csv...

@Llopis Quite true, but the readr package has that as the default :)

As Devon said, read.delim can deal with this perfectly well as long as the file is properly formatted (so missing fields still end with the field separator, like \t\t). For example, with this input:
Field1  Field2  Field3
Row1field1  Row1field2  Row1field3
Row2field1      Row2field3

Or, to show the fields more clearly:
Field1\tField2\tField3
Row1field1\tRow1field2\tRow1field3
Row2field1\t\tRow2field3

R can import:
> df<-read.delim("file",header=TRUE, sep="\t")
> df
      Field1     Field2     Field3
1 Row1field1 Row1field2 Row1field3
2 Row2field1            Row2field3

Alternatively, if you really want to add something to the missing columns, you can use:
$ awk -F"\t" '($3==""){$3="none" }1;' file
locus_tag   ftype   gene    EC_number   product
BOFHCHGE_00001  CDS hypothetical    protein
BOFHCHGE_00002  CDS rdgB_1  3.6.1.9 dITP/XTP    pyrophosphatase
BOFHCHGE_00003  CDS hypothetical    protein
BOFHCHGE_00004  CDS hypothetical    protein
BOFHCHGE_00005  CDS hypothetical    protein

Or, to add none to any empty field, not just the 3rd one:
awk -F"\t" '{for(i=1;i<=NF;i++){if($i==""){$i="none"; }}}1;' file > newfile

I fear that the .tsv file is mal-formatted then. df<-read.delim("file",header=TRUE, sep="\t") will misimport and EC column will end up having product names, gene having hypothetical protein, etc.

@charlesdarwin if it is indeed malformatted, yes. However, I would be very surprised if that were the case. Using empty fields as described above is standard practice. You can check easily enough though. FInd one of the lines with a missing field (using 12 here as an example), and then run awk 'NR==12' | od -c and make sure you see the \t\t where the empty field is. If you don't see that, this becomes a text parsing problem. I suggest you post a new question with some example lines and we'll see what we can do.

my 6th line has a missing EC number. My file is file.txt and my data frame is called df, what is the command I should run in the terminal with awk?

@charlesdarwin (sorry, I just saw this for some reason) it depends. If there is a \t\t for the empty field, you can do awk -F"\t" '($4==""){$4="none" }1;' file.txt > newfile.txt. If that doesn't work, you might need to ask a new question, or come into the site chatroom and ping me and I'll see if I can help.

newfile.txt is a mess. The element 1,1 is now BOFHCHGE_00001 CDS hypothetical protein none. I've written to the developer. I don't think it should be that difficult. I am using the GFF file instead. It is well formatted and contains all the information, too. Thanks for your help.

@charlesdarwin that's probably for the best. Sorry, but it's really hard to parse something without having the original file to be able to see exactly what's needed.

very easy to reproduce, download the E. coli genome from here: ftp://ftp.ensemblgenomes.org/pub/release-32/bacteria//fasta/bacteria_91_collection/escherichia_coli/dna/Escherichia_coli.HUSEC2011CHR1.dna.toplevel.fa.gz, run prokka on it (conda install prokka) and you'll get a malformed .tsv file.

@charlesdarwin OK. Sorry, but I don't use prokka and I don't want to install it for this. If you want to paste a few lines of the output file in chat or in your question, I can try to help.

Here is a link to the offending file: http://cristian-riccio.ch/wp-content/uploads/2017/12/PROKKA_12142017.tsv

@charlesdarwin oh wow. Yeah, that's definitely a bug, the fields are wrong and since different lines are missing different fields, I dont' see how you can change it unless you do so manually. That's really impressively bad design.

Since the missing values don't have field separators, the easiest option is to rearrange the columns so that the ones with missing values are at the end. They will then be the ones that are filled with missing values when you import them into R. You can use use fill=TRUE with read.table() or use readr, which will complain about missing columns, but fill them in with missing values automatically:
awk -F'\t' -v OFS='\t' '{print($1,$2,$5,$3,$4)}' file.txt |perl -pe 's/\t+/\t/g' >newfile.txt

With readr:
read_tsv("newfile.txt")

# A tibble: 5 x 5
       locus_tag ftype                  product   gene EC_number
           <chr> <chr>                    <chr>  <chr>     <chr>
1 BOFHCHGE_00001   CDS     hypothetical protein   <NA>      <NA>
2 BOFHCHGE_00002   CDS dITP/XTP pyrophosphatase rdgB_1   3.6.1.9
3 BOFHCHGE_00003   CDS     hypothetical protein   <NA>      <NA>
4 BOFHCHGE_00004   CDS     hypothetical protein   <NA>      <NA>
5 BOFHCHGE_00005   CDS     hypothetical protein   <NA>      <NA>

The perl regex is necessary because awk adds tabs after the missing values, which defeats the point of reordering the columns

I did as you said but the data frame is still wrong. On line 6, gene and product are interchanged.

Since you only provided 5 lines in the example, I can't comment on the results for the 6th line

all. I know this question was asked 2 years ago, but I wrote a script in R, with a little bash, to parse the tsv file in the old version of prokka. In the new version of prokka this problem is fixed.
First use bash to get the rows where ftype is CDS.
grep "CDS" all_prokka_annotations.tsv > all_prokka_annotations_CDS.tsv

Next we can parse the TSV with the CDS lines in R.
prokka_annotations = readLines("all_prokka_annotations_CDS.tsv")
prokka_annotations_split = strsplit(prokka_annotations,split="\t")
# first get the protein product name.
prokka_annotations_product = sapply(prokka_annotations_split, function(x) x[length(x)])
# now the I have all protein product names remove them
prokka_annotations_split = sapply(prokka_annotations_split, function(x) x[-length(x)])
#next get locus tag and ftype
locus_tag_cds = sapply(prokka_annotations_split, function(x) x[c(1,2)])
locus_tag_cds = t(locus_tag_cds)
# now remove locus tag and ftype
prokka_annotations_split = sapply(prokka_annotations_split, function(x) x[-c(1,2)])
# next get EC number
ec_number = lapply(prokka_annotations_split, function(x) x[grep(".",x,fixed=T)])
ec_number = as.character(ec_number)
ec_number[ec_number=="character(0)"] <- NA
## now get the gene names if available. to do this remove the EC numbers. whats left are gene names
prokka_gene_names = lapply(prokka_annotations_split, function(x) x[!grepl(".",x,fixed=T)])
prokka_gene_names = as.character(prokka_gene_names)
prokka_gene_names[prokka_gene_names=="character(0)"] <- NA
prokka_parsed_tsv_correct = cbind(locus_tag_cds,product=prokka_annotations_product,gene=prokka_gene_names,EC=ec_number)

```