text:
How do I efficiently perform a metagenome screen of “all” species?
I’ve got an RNA-seq dataset with a large proportion of environmental RNA “contamination”. BLASTing random reads reveals that much of the data comes from bacterial, plant and viral RNA. My target organism only accounts for ~5% of the RNA-seq read data.
I would like to get a more or less (more, if possible) comprehensive overview of what species are found in my sample — bacteria, plants, animals (?), but also viruses. How can I perform this?
I have been unable to find a “standard” way of performing a metagenome screen. BLAST online services all seem to be severely rate limited (certainly unable to handle uploading ~80 M reads). Installing BLAST (or a similar tool) locally is of course not an issue, but I can’t find a comprehensive database that spans all phyla — the best I’ve been able to find are databases, such as NCBI-NR, that are restricted to single phyla or classes, or to bacteria.
Ideally I’d like a workflow that I feed with my RNA-seq data, and get out a list of species with % of reads mapped.

NCBI's nr database isn't restricted to any species as far as I know.  Why do you say it's limited to  single phyla or classes or to bacteria? Wouldn't blasting against nr be good enough? And are they really that rate limited? I admit I haven't done this in a few years, but back in the day, I could use NCBI's command-line blast client which would use a remote nr DB and could happily handle several thousand sequences. Have you tried that? Perhaps by splitting your input reads?

@terdon From what I’ve read, it’s only for procaryotes. But I may be completely wrong. Regarding the rate limit, the NCBI BLAST FAQ states that they don’t support batch usage, and that large jobs shouldn’t be run through the URL API.

Nr covers a wide range of taxa as @terdon says. You can download it and build your own blast db if you have ~100gb disk space and a medium memory server. However, BLAST will be slow, I would recommend DIAMOND instead

As you point out, BLAST is not going to work well for reads from a typical RNASeq experiment.
From an RNASeq dataset, you might be able to get away with just a ribosomal screen using Silva, because the ribosome is the most abundantly transcribed gene. This will usually be the case even if polyA selection or ribosomal depletion was carried out on a sample.
However, bear in mind that there will be a lot of false results due to database incompleteness. I used Kraken when I did something like this a few years ago, and ended up with one sequence annotated as a fruit bat (where my subject was a water-living planarian flatworm). Centrifuge would be a better option now; the authors even claim that they can store the entirety of nr in a 70 GB index with Centrifuge (which would be better / more comprehensive than just a ribosomal screen).

One thing that hasn't been suggested so far would be phylogenetic placement. First option there would be placing against one larger tree roughly touching the domains you mentioned to get an overview. Second option would be having three trees (bacteria, plants, viruses) made up of relevant taxa for the type of sample you have.
This would give you a detailed and accurate view of your sample, but it does involve some work. Sadly the computational throughput for placement is still a bit behind BLAST, and not many complete pipelines exist, but in my opinion it is worth the extra mile (I might be biased. I'm working to resolve as many of these issues as possible).
If you want to look into it, I can reccomend this excellent paper by Erick Matsen. Since your samples sound vaguely soil-ish, I will also link this recent Nature paper using placement to analyze rainforest soil samples.

I think you're just looking for nr. It is absolutely not limited to prokaryotes, far from it. It is, according to the blurb on the NCBI's blast page:

The nucleotide collection consists of GenBank+EMBL+DDBJ+PDB+RefSeq sequences, but excludes EST, STS, GSS, WGS, TSA, patent sequences as well as phase 0, 1, and 2 HTGS sequences. The database is non-redundant. Identical sequences have been merged into one entry, while preserving the accession, GI, title and taxonomy information for each entry.

The link you gave in your question doesn't describe the nr database, but the collection of RefSeq non-redundant sequences. The two are not the same thing at all. The nr databases contains sequences from all available species so it should be perfect for what you're after.
Now, I admit I haven't tried this with millions of redas, but I did used to run many blast searches with several thousand sequences using the command line blast tools and the remote nr database. I just checked the blast+ docs now and see no mention of a rate limit. Which isn't to say there isn't one, for sure, but since each query is run sequentially when you give a multifasta file, there seems to reason to limit since the servers should be able to queue the jobs as needed.
So, don't use the API or roll your own, just download NCBI's blast client and try that:
blastn -query query.fa -db nr -task blastn

If that complains, you'll have to download the nr DB (well, since you want nucleotide blast, you probably want its nucleotide cousin, nt) which you can find here: ftp://ftp.ncbi.nlm.nih.gov/blast/db/. Just get all the foo.nt* files.

According to the BLAST API Usage guidelines, this isn’t really supported, even if there is no hard limit. But you seem to be right about NR, this seems to be what I’m after, after all.

@KonradRudolph I'm hoping that refers to using the API in your own tools and not to using their own standalone blast+ program. I may very well be wrong, but I'm not sure if the instructions for the API apply to them as well. I would expect them to have built their tool in such a way as to not overload the servers. I figure it's at least worth a try.

You might want to look into minihashes as a more efficient way of clustering and comparing your sequences to known species and samples. In the mash paper the authors show how it would be feasible to create minihash sketches that you could search against from large metagenomics datasets (check the "Clustering massive metagenomic datasets" section in the results):

For a large-scale test, samples from the Human Microbiome Project [36]
  (HMP) and Metagenomics of the Human Intestinal Tract [37] (MetaHIT)
  were combined to create a ~10 TB 888-sample dataset. Importantly, the
  size of a Mash sketch is independent of the input size, requiring only
  70 MB to store the combined sketches (s = 10,000, k = 21) for these
  datasets.

Moreover, the creators of sourmash (a tool using the same approach) provide precomputed sketches for the whole of NCBI's RefSeq database (includes bacterial, viral and fungal genomes). Creating a similar sketch for the whole of ENA/SRA (or a subset of it) could require some time but would allow for very fast and accurate searches.

I have had good success using the MEGAN software to detect and visualize contamination.  My input data was WGS not RNA-seq But I would check it out.  It can create a very nice weight phylogenetic tree for all the detected organisms.
https://ab.inf.uni-tuebingen.de/software/megan6

I've downloaded MEGAN but my (admittedly cursory) reading of the documentation made it seem as if the software only visualised results computed by other tools.


text:
ngs time course experiment
I have a question regarding the RNA-Seq analysis.
We have a time course data regarding a mouse model wt and mutant treated with a drug (10uM) and the taken down at different time points: 8 time points and did RNA-seq on these samples currently we have tpm values and raw gene counts from rna seq. Also we have only one sample for each time point. we would like to see for changes in genes across different time point but not sure statistical method should we use to determine this, should we consider for example a single time point wt vs mt, or see trends across wt and mt individually
I know that there are some bioc packages for NGS analysis (edgeR,
DESeq,etc). Is any package or method good for one-sample-per-time-point
data?

For "one-sample-per-time-point data" there isn't any good statistical test available. There are tricks (such as described in edgeR manual page 21), but there is no sound statistical ways to analyze these poor designs. I wonder why you choose such a poor design? And not ask a statistician or bioinformatician first before you do the (expensive) NGS?

I would recommend to do at least MDA plots to see if there is visible condition / time signal. However, I agree that 1 replicate is very bad idea. Is there chance to get more replicates?

@IanSudbery another poor soul like me I guess, definitely a duplicate need to up my search skills

There is a discussion here about how to approach this analysis using limma: https://support.bioconductor.org/p/33783/

source: Stack Exchange

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

sort: votes


text:
Biopython: resseq doesn't match pdb file
I have a PDB file, and I need to extract its residue sequence numbers (resSeq's). Based on manual inspection of the first few lines of the PDB file (pasted below), I would think that resSeq's should begin with 22, 23. However, Biopython's PDB module suggests otherwise (output attached below as well). I wonder if it's a Biopython bug or if I have problems understanding the PDB format.
ATOM      1  N   GLY A  22      78.171  89.858  59.231  1.00 21.24           N
ATOM      2  CA  GLY A  22      79.174  88.827  58.999  1.00 20.87           C
ATOM      3  C   GLY A  22      80.438  89.415  58.391  1.00 21.89           C
ATOM      4  O   GLY A  22      80.362  90.202  57.440  1.00 23.18           O
ATOM      5  N   LEU A  23      81.588  89.069  58.972  1.00 21.51           N
ATOM      6  CA  LEU A  23      82.895  89.555  58.527  1.00 20.80           C
ATOM      7  C   LEU A  23      83.288  89.020  57.162  1.00 22.41           C
ATOM      8  O   LEU A  23      82.889  87.923  56.788  1.00 22.93           O
ATOM      9  CB  LEU A  23      83.973  89.232  59.560  1.00 20.97           C
ATOM     10  CG  LEU A  23      84.225  87.818  60.062  1.00 13.32           C
ATOM     11  CD1 LEU A  23      85.448  87.888  60.939  1.00 15.24           C
ATOM     12  CD2 LEU A  23      83.035  87.258  60.829  1.00 12.21           C

...
for i in chain:
    print i.get_full_id()

OUT:('pdb', 0, 'A', (' ', 2, ' '))
    ('pdb', 0, 'A', (' ', 3, ' '))
...

Please post the complete relevant code. I.e. how is chain exactly defined?

And assuming this is a public PDB file, what is the four character identifier so people can try the example themselves to help you.

Let's check it. I copied your file fragment into a.pdb.
>>> from Bio import PDB
>>> for res in PDB.PDBParser().get_structure('pdb', 'a.pdb')[0].get_residues():
...   print(res.get_full_id())
...
('pdb', 0, 'A', (' ', 22, ' '))
('pdb', 0, 'A', (' ', 23, ' '))

So BioPython (ver. 1.66) gives 22 and 23 as expected.

Thank you! I'm soo sorry, I've found a bug in my code after 2 days of searching, and everything's working now (I forgot that I passed a pre-processed PDB to my program rather than the original one) . Thank you for your time and sorry that it's not a very useful question! Do you think I should delete it?

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


text:
Mathematical notation for formulating a rank score
This is a cross-post from Mathematics forum. As no one has answered it yet, I fIgured post here as well.
I would like to describe a transform I used to rank my data points. I have recoded my variable columns from $10$ measurements with different range to $n$ (for each measurement) where $n \in \{0,0.5,1\}$.
For columns where measurements were expressed as categorical values I simply recoded the categorical value to either $0$, $0.5$ or $1$ and for continuous variables I recoded the lower quartile as $0$, interquartile as $0.5$ and upper quartile as $1$ and finally I summed the recoded values to produce a single score $\sum n_{i1,..i10} $ for each row.
I am trying to write my method and I'd like to have it in mathematical notation. I was wondering if anyone could help me with this!

Maybe try https://stats.stackexchange.com/

What do you mean by "have it in mathematical notation"? Can you give an example of something you consider "in mathematical notation"?

Also, it would be good to include the link to the cross-post

Hi bli, Thank you for your comment. By mathematical notation I mean mathematical representation of my transform. I am including the link as you suggested!

Your question is not very clear to me. So can you please try to explain it with some example. Then we may able to help you.

Have you looked at the dataframe I uploaded in Mathematics forum?

Yes, I did. But still it not clear.

I think the problem with your question is two-fold. 1) you don't provide any reproducible data so that we can grasp the details of how you calculate this "rank" score. See this post for an example how you would make a reproducible example in R to get an idea of what we mean by "reproducible example" 2) there is not likely to be a mathematical notation for what you have done, because there is no mathematical formula to convert strings ("I", "II", "III") into numbers. One way that you could potentially solv

I think that "mathematical notation" does not necessary have to consist in "formulas". One could probably define a bijection between a set of strings and a set of numbers, in proper mathematical notation, and that would validly describe a conversion procedure.

Yeah you are probably right that there is some valid way of describing the rank score. I imagine in this case it would be something more akin to a key than notation/formula, but if you have an idea how to describe it why don't you add an answer?

As you pointed out, it is difficult to grasp the details of what the OP wants. Besides, it's been a long time since I last wrote maths formally.

Here is my solution to the problem. I am posting it here in case anyone else come up with the same idea but did not know how to formulate it in mathematical notation!

$$ \psi_i = \sum_{i=1}^9 S_i $$
where $S_i$ is the score function and defined as:
$$
S_i = \cases{
  0 & \text{if $\theta(x_i) < Q_1$ or $\theta(x_i) = $ benign/neutral} \\
  0.5 & \text{if $Q_1 \leq \theta(x_i) < Q_3$ or $\theta(x_i) = $ possibly damaging/uncertain} \\
  1 & \text{if $\theta(x_i) \geq Q_3$ or $\theta(x_i) = $ damaging}
}
$$
The $\theta(x_i)$ is the pathogenicity or conservation score for variant $x$ as defined by model $i$ and $Q$ denotes the quartile range for scores from M-CAP, CADD, GERP and Phylop models.

If Q are quartiles, then the last line for damaging shouldn't be $>Q_3$?

@Konrad and Reza : Can you look into this problem and provide some suggestions. https://math.stackexchange.com/questions/2716192/how-to-convert-a-detailed-mathematical-statistical-compuational-process-into-ma


text:
Best way to detect long insertions in bisulfite sequencing data?
I am interested in identifying indels in whole genome bisulfite sequencing data (76bp paired end). Currently, I do this by setting the -rfg and -rdg affine gap penalty scores for bowtie2 to more permissive values than the default 5+3N and mapping using the bisulfite sequencing alignment wrapper, bismark.
My question is, what values of -rfg and -rdg will allow me to identify longest possible indels without sacrificing alignment quality? Is it better to set the affine penalty to zero with a high penalty for initially opening an indel (ex. 8+0N)? Or is it better to keep the initial penalty low and having a nonzero penalty for extension (ex. 2+1N)?

You'd be best off by starting with -rfg and -rdg as is and reverting bismark's change of --score-min back to the default for bowtie2. That alone will allow for much longer indels. If that still doesn't suffice, then I'd play around more with --score-min before messing with the gap open/extend penalties. If you do need to play with those, then increase the gap open penalty and decrease the gap extension to 1. Do check that the resulting alignments aren't nonsense though!

With the bismark default --score-min, -rdg, and -rfg settings, it appears that the maximum sized indel identifiable with my read length is 3bp. I changed both -rdg and -rfg to a gap open/extension of (1,1) and the maximum length I can detect is 14bp. I will try modifying --score-min only. It's a fine line to walk, detecting the longest possible insertions while keeping the alignment clean!


text:
Multiple sequence alignment with R
I want to incorporate Multiple sequence alignment into our Shiny app, i.e. allow the possibility of selecting different nucleotide sequences, and do a multiple sequence alignment in real-time. The user would define as input different nt. sequences and then a tree will pop up in the app. Do you think this is feasible?
Any recommendations and feedback are really appreciated :)

What do you mean by “tree”? Trees don’t really make sense for pairwise alignments.

I agree with Konrad Rudolph, trees don't make sense with pairwise alignments. There is a pairwiseAlignment function in the Biostring package though, the writePairwiseAlignments funtion can show the alignment.

I mixed up terms here, I meant multiple sequences alignment not pairwise. I have edited the question, sorry about that!

@b.nota Pairwise alignments could be used to compute distances, and a tree could be made from the distance matrix.

@bli Yes, but does it make sense?

@b.nota I think that in theory, a multiple alignment is more apt to capture relevant homologies. So if the tree is meant to represent the phylogeny of the sequences, multiple alignment is theoretically better. That said, there might be practical reasons for which one would have to work with pairwise distances instead.

@bli, okay thanks for your thoughts. I do still agree, however, with Konrad that trees for pairwise alignments don't really make sense.

@SergioArredondo, as the question is closed I cannot provide an answer, but a relevant newer question provides a MWE that shows an outline for the situations where the nucleotides are stored in an msf format, but that can be changed accordingly as other formats can be specified https://bioinformatics.stackexchange.com/a/14160/9311 (and real time is tricky given some of the datasets can be large)

The ‹msa› package on Bioconductor does exactly that.
It doesn’t hand the results on a silver platter, though: you’ll need to read the vignette carefully to learn how to use it. But after that it’s pretty powerful.

"It doesn’t hand the results on a silver platter", that is very true, but some boiler plate code would go a long way since many request the same functionality from a same set of functions


text:
How to calculate overall reference coverage with MUMmer?
Is the MUMmer suite capable of calculating reference sequence coverage statistics for all query sequences collectively? It would be possible to achieve by parsing the output of nucmer / show-coords / show-tiling but it seems like there should be a better way. I currently do this using a sensitive read mapper, samtools depth and some scripting.
To clarify, I'd like to know the reference coverage achieved using all of the query sequences (i.e. the whole de novo assembly, in my case).

Would QUAST work on your data? It runs MUMmer and reports coverage.

Good suggestion @user172818. I'll investigate.

I believe all you need to do is to run dnadiff from MUMmer. That will run a comparison and output a number of useful metrics.

Many thanks Devon. The QRY field of the AlignedBases record in dnadiff's out.report seems to be what I was after. Although clearly a powerful tool, I was surprised how slow dnadiff was to run.


text:
FASTQC overrepresented sequences after trimming
I have a set of RNA-seq samples from different experiments (Single and Paired End, depending on the experiment). I ran FASTQC in all the samples and found overrepresented adapter sequences:
I removed the adaptors (TruSeq adaptors) using Cutadapt (in addition, I removed low quality and N bases from the 3' end of the reads). After that, I ran again FASTQC and the output is the following (representative example) :

Does anyone know what is happening? Now I have an overrepresented sequence for which no sequence is provided. What does this mean?

Looks like you trimmed your sequences to nothing.  Do you have empty fastq entries in the file?

Yes, thanks for your comment, this was the problem.

As @AaronBerlin mentioned, you didn't remove reads that were completely trimmed. Next time use the --minimum-length option and set it to something reasonable, like 20. Alternatively, use "Trim Galore!", which is a wrapper around cutadapt that has more reasonable defaults.

Indeed this was the problem, I set the --minimum-length parameter to 20 and everything is fine. I didn't try Trim Galore! because I am working in a cluster and it is not installed at this moment but in the future I will try it as it seems more intuitive.


text:
Why sequence the human genome at 30x coverage?
A bit of a historical question on a number, 30 times coverage, that's become so familiar in the field: why do we sequence the human genome at 30x coverage?
My question has two parts:

Who came up with the 30x value and why?
Does the value need to be updated to reflect today's state-of-the-art?

In summary, if the 30x value is a number that was based on the old Solexa GAIIx 2x35bp reads and error rates, and the current standard Illumina sequencing is 2x150bp, does the 30x value need updating?

The earliest mention of the 30x paradigm I could find is in the original Illumina whole-genome sequencing paper: Bentley, 2008. Specifically, in Figure 5, they show that most SNPs have been found, and that there are few uncovered/uncalled bases by the time you reach 30x:
These days, 30x is still a common standard, but large-scale germline sequencing projects are often pushing down closer to 25x and finding it adequate. Every group doing this seriously has done power calculations based on specifics of their machines and prep (things like error rates and read lengths matter!).
Cancer genomics is going in the other direction. When you have to contend with purity, ploidy, and subclonal populations, much more coverage than 30x is needed. Our group showed in this 2015 paper that even 300x whole-genome coverage of a tumor was likely missing real rare variants in a tumor.
On the whole, the sequence coverage you need really depends on what questions you're asking, and I'd recommend that anyone designing a sequencing experiment consult with both a sequencing expert and a statistician beforehand (and it's even better if those are the same person!)

This is as far as I could trace back, but would love to know if anyone can find an earlier mention of 30x as a standard depth!

That is likely the earliest.  Previous technologies could not reach that depth at a reasonable cost.

Solexa Inc. sequenced NA12878 chrX to ~30x in early 2007, which later became part of Bentley (2008). This, I believe, was the first time that 30x showed up. I don't recall they had a particular reason for that. Figure 5 in the published paper was just aftermath. It does not really explain why not 25x or 35x, given that the curves between 25x and 35x in that figure are about linear.
In the abstract of Ajay et al (2011), the authors argued "the current recommendation of ~30x coverage is not adequate". Nonetheless, the discussion section seems to suggest 50–60x would be necessary with GAIIx, but 35x was adequate with HiSeq2000 plus better recent chemistry. Overall, this paper provides a more thorough analysis. The data quality at that time is also closer to data we produce today.
The required coverage is largely determined by two factors: read placement bias (e.g. GC bias) and base/mapping error rate. While GC bias has been reduced with the PCR-free protocol, base error rate has been going the downward since HiSeq2500. I guess 30x coverage would be necessary if you want to achieve the sensitivity with the older 30x data. Illumina, as a sequencing service provider, and our sequencing facility still insist on the 30x threshold.

Interestingly enough we also found in Sudbery et al, 2009 (https://genomebiology.biomedcentral.com/articles/10.1186/gb-2009-10-10-r112), that the number of continuous regions without a break in coverage was more or less linear in the number of sequenced reads between 25x and 60x when sequencing the mouse genome.

30 times coverage is not unique to this problem, but number 30 has its empirical role in statistics:

In statistical analysis, the rule of three states that if a certain event did not occur in a sample with n subjects, the interval from 0 to 3/n is a 95% confidence interval for the rate of occurrences in the population. When n is greater than 30, this is a good approximation to results from more sensitive tests.

source: Wikipedia: Rule of three (statistics)
Similarly, you can search for related questions like this one:

What is the rationale behind the magic number 30 in statistics?

In line with this, I have seen data processing in other disciplines which required n ≥ 30 for sufficient reliability of results.

In this type of community, I would expect some opponenture under the answer but I am getting only downvotes so far (curent score: +6/-3).

The point I always miss in the discussion about coverage is, that no one tells how it was calculated. Were duplicates removed? How do you count overlapping paired-end reads? As 2 or 1? Just to point out two things that have influence on the coverage.

Anyone who counts overlapping paired-end reads as 2 is doing it wrong. There's never an excuse for that.

@DevonRyan why? Overlapping parts of pair-end reads are two independent technical replicates of the same region. What is the difference to overlap of two reads that are not paired, should they also be counted as one?
Of course if you merge overlapping reads into longer sequences before an assembly, then you got to count it once, but otherwise I do not see a reason why they should be counted as one...

@KamilSJaron They represent one sequenced fragment, counting them otherwise is lying.

@DevonRyan Same fragment, but different sequencing and it is sequencing coverage, not a fragment coverage. I do not get your point.

@KamilSJaron You're being overly literal. The point of the metric is to assess how much data you have per position. Overlapping PE reads don't represent different data points, that's why they're treated as a single unit in variant and peak calling.

In my opinion the coverage should represent the number of different molecules. So counting overlapping paired reads and removing duplicates is mandatory. Otherwise the read counts is just a more or less useless number.

@KamilSJaron, if you want to utilize that information, there are algorithms that collapse the overlapping segments of such reads in an intelligent way (e.g. for each overlapping base, use the one with the highest quality and discard the other). Devon is right that they absolutely should not be treated as independent events for variant discovery, depth calc, etc.


text:
How to manipulate a reference FASTA or bam to include variants from a VCF?
I have some software which takes fastas as the input. I need to include SNVs and InDels from a VCF into the reference hg38 and then use this.
The problem is, I don't know of an algorithmically sound way to do this.

Are there any existing software packages which could do this efficiently? Is it easier to output a FASTA, or a bam (and then convert to a FASTA)?
What about if I wanted to do the same with a bedpe of germline structural variants?

Could you provide an example of what you're looking for?  I'm having a hard time picturing what you want.

@Greg The above might be unclear. I have a reference FASTA. I have variants. I would now like to manipulate the reference FASTA to include these variants. The VCF tells us there's a SNP at XX location; I would like to like a FASTA that includes these SNPs. Does this make sense?

Oh I understand now.  Here's a link to a similar biostars question that has some answers that look good.  https://www.biostars.org/p/6553/

You could convert VCF to BED via vcf2bed --snvs, vcf2bed --insertions, and vcf2bed --deletions, and then use samtools faidx by way of a wrapper script to convert BED to FASTA, e.g.:
$ vcf2bed --snvs < variants.vcf | bed2faidxsta.pl > snvs.fa
$ vcf2bed --insertions < variants.vcf | bed2faidxsta.pl > insertions.fa
$ vcf2bed --deletions < variants.vcf | bed2faidxsta.pl > deletions.fa

You need FASTA files for your reference genome, which have been indexed with samtools faidx, e.g., for hg38:
$ cd /foo/bar/baz
$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
$ for fn in `ls *.fa.gz`; do gunzip $fn; done
$ for fn in `ls *.fa`; do samtools faidx $fn; done

Once you have indexed FASTA files somewhere on your file system, you can pipe BED to the bed2faidxsta.pl script, to get out FASTA sequences.

GATK has a solution that might work for you:
FastaAlternateReferenceMaker, which : "Given a variant callset, this tool replaces the reference bases at variation sites with the bases supplied in the corresponding callset records."

Input
The reference, requested intervals, and any number of variant ROD files.
Output
A FASTA file representing the requested intervals.
Usage example
java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -L input.intervals \
   -V input.vcf \

There's a vcf2fq sub-program that was written as part of vcfutils to convert a VCF file into a fastq file given a reference sequence. Unfortunately this doesn't work properly with INDELs (it will just mask them, rather than actually converting them), so I wrote a modification to implement INDEL correction as well:
./vcf2fq.pl -f <input.fasta> <all-site.vcf> > <output.fastq>

It adds and removes sequence from the reference, as specified by the INDEL information in the VCF file. The code keeps a record of where in the reference sequence the INDELs should go, then inserts (or deletes) the INDELs, updating the reference position as necessary. It should work with things like deletions right next to insertions, but hasn't been extensively tested.

Could you explain how your perl script with the INDEL implementation works? i.e. what changes did you make for vcfutils.pl to work properly

"It should work with things like deletions right next to insertions, but hasn't been extensively tested." If you do benchmark this, let me know---this would be interesting to see

Another question: have you seen vcfutils.pl written in Python/R? I suspect these exist somewhere, as I don't want to re-invent the wheel. Is this worth a separate SO question? It would be good if your perl script was more advertised IMHO

I think an answer to that fits in with this question. The question is about manipulating FASTA files, not about vcfutils.
No, I haven't seen any python/R alternatives. I don't see much of a point in them (certainly not a python version), given that the command-line functionality already exists via vcf2fq.

You're welcome to have a go at advertising this further. I've brought it up a few times on SeqAnswers, github, and I think reddit; basically every time I see a question asking something like this.

Here's an update how this work in 2019. (The answer is quite similar to my update in this thread):
bgzip and index your vcf file.
$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz

Now use bcftools consensus:
$ bcftools consensus -f genome.fa input.vcf.gz > consensus.fa

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


text:
Annovar index out of date
I'm trying to use Annovar to annotate some variants with their CADD and FATHMM scores. I've downloaded the latest versions of the software and the databases but when I run it I get an error saying the index is out of date.

WARNING: Your index file annovar_latest/humandb/hg19_cadd.txt.idx is out of date and will not be used. ANNOVAR can still generate correct results without index file.

It still runs without the index file but it is extremely slow. How can I update the index file?

Your safest bet is to just redownload the annotation. I've not seen official documentation of the index used by annovar, but apparently it's just a text file with a single line header and chromosome-bin coordinates. Over on SEQanswers, there's a thread discussing this issue with a single perl script that allegedly recreates the index.

I just downloaded them again this morning, so that doesn't seem to work.  I have a job running the perl script fro that post to generate new index files.  If that works I'll report back.
I don't see how they could possibly be out of date because there is nothing in there other than row numbers.

I downloaded everything again and it worked.  Sometimes the simplest solution is the best.  Thanks

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


text:
How can I edit a specific FASTQ read in place, given the read ID?
I am introducing SNVs into specific samples in order to estimate false negative rates for a variant calling pipeline. I know reads can be simulated but I would actually prefer to use the real data so as to keep everything else equal.
Given a read ID, I want to edit a single basecall (e.g. the 12th base) for just that read from within a large FASTQ file containing millions of reads.
example, i want to change the 12th base ('C') in read 31027 to a 'T':
@70630 1:N:0:ATCACG
GAAGGTCCATGGATAATACTCAATTTTCCACAACAGCTTTTGTACTCTAGATCATTGATATTTACCAAAAGTCACTCAAACTCATCCTATGCATAATTCTAGTCCACCAATCATGATATGATGGAGAACATGGTTGTAATCAGGAAGACAG
+
DDDDDHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIGHHIIIIIIIIIIIIIIIIIIIIIIHIIHIHIIIHHIGIIIIHHHHIIH@
@31027 1:N:0:ATCACG
CAAAAGTCACTCAAACTCATCCTATGCATAATTCTAGTCCACCAATCATGATATGATGGAGAACATGGTTGTAATCAGGAAGACAGATAAAGCAGCAGACCAAAAGTAATCTGAGAAATTATATTTGAATCACTCAGATATACATCAAATA
+
DDBDDHC@GHFHHIIE@CEHHIHEFFCGFCH?HHEGFCCEHHCGH@HFHHHCHHIIIHEHIIHII@CECHIHIIIEECCDGEHFFHHHHEHHFHHGGIHHHDGHFHIIIGHHHHHHEHHIIHIIIGGHHHCHHHCHIIIHHIEH@GHIHIC
@87319 1:N:0:ATCACG
CAATTAAGCTTTGGCAACGGTGGTCAAGATGAGATGCATATGGAGATAATAACTAAAAGTCAATCGAGACTCATCGTATGCATATTTCTAGTCCATCGATCATGAAATGATAGGATAGCTAGAATGAAAAGTAAATTTCCAGAAGGTCCAT
+
DDDDDIIIIIIIIIIIIIHIHHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHHIIIIIIIIIIIIIIIIIHIIHIIIIIIIIIIIIHIIIIIIIIIHIHIH

Of course I could stream through the entire FASTQ file rewriting everything to a new file until i get to the read of interest, do the edits to that read, write it out to the new file, and then keep going until the end. BioPython would work fine for this as would other approaches.
However, is there an efficient way to do this without the full read/write stream? Could this be done with sed/awk ? What about an indexed fastq file?
cheers
David

You'll first need to determine the appropriate line number, which you can do with grep -m1 -nw "@31027" foo.fastq. After that, note that you can provide a line number to sed:
sed -i '123456s/CAAAAGTCACTCA/CAAAAGTCACTTA/' foo.fastq

That will do the replacement only on line 123456 and edit the file in place (the -i option).

I was not aware of the -n option for grep. This should do the job

@dkainer just for completion's sake, you can get the line number of the line(s) matching the string foo with sed -n /foo/= or perl -ne 'print $. if /foo/ or awk '/foo/{print NR} as well.

If your fastq file is simple, if each of your reads only has a single line of DNA, you could do:
awk '{
        if($1 == "@31027"){
            a=NR+1
        }
        if(NR==a){
            split($0,s,"");
            s[12]="T";
            for(i in s){
                seq = sprintf("%s%s",seq,s[i]);
            }
            $0=seq
        }
    }1;' file.fastq


text:
How to adapt the fgseaL function to perform rapidGSEA computation of gene ranks across 9 different phenotype labels?
I wish to adapt the r language function fgseaL, https://github.com/ctlab/fgsea , to perform rapidGSEA, https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-016-1244-x , computation of inter-class  deviation per gene and the subsequent gene rank sorting operation on 9 different phenotype labels as illustrated in the diagram immediately below:

I thought of applying the R-language rank() function on the Expression Data Matrix D. If that is not correct, what sequence of R language commands should we
apply to the Expression Data Matrix D to calculate a key value sorted deviation measure  across 8 labeled human leukemia groups and a healthy labeled normal control group prior to running fgseaL?
I show below how fgseaL finds the correlation matrix between the R language variable , mat, which corresponds to the Expression Data Matrix D and the R language variable , labels , which is a vector of gene phenotype labels
    tmatSc <- scale(t(mat))
    labelsSc <- scale(labels)[, 1]

    minSize <- max(minSize, 1)

    pathwaysFiltered <- lapply(pathways, function(p) { as.vector(na.omit(fmatch(p, rownames(mat)))) })
    pathwaysSizes <- sapply(pathwaysFiltered, length)

    toKeep <- which(minSize <= pathwaysSizes & pathwaysSizes <= maxSize)
    m <- length(toKeep)

    if (m == 0) {
        return(data.table(pathway=character(),
                          pval=numeric(),
                          padj=numeric(),
                          ES=numeric(),
                          NES=numeric(),
                          nMoreExtreme=numeric(),
                          size=integer(),
                          leadingEdge=list()))
    }

    pathwaysFiltered <- pathwaysFiltered[toKeep]
    pathwaysSizes <- pathwaysSizes[toKeep]

    corRanks <- var(tmatSc, labelsSc)[,1]
    ranksOrder <- order(corRanks, decreasing=T)
    ranksOrderInv <- invPerm(ranksOrder)
    stats <- corRanks[ranksOrder]

    pathwaysReordered <- lapply(pathwaysFiltered, function(x) ranksOrderInv[x])

    gseaStatRes <- do.call(rbind,
                           lapply(pathwaysReordered, calcGseaStat,
                                  stats=stats,
                                  returnLeadingEdge=TRUE))

I found a problem with the algorithm shown immediately below.
correcttest <- data.frame(names = row.names(normal))
correcttest <- cbind(correcttest3, normal)
correcttest <- cbind(correcttest3, ALL3m)
rownames(correcttest) <- correcttest$names
correcttest$names <- NULL
correctlabelnormal <- rep(0:0, 73)
correctlabelALL3m <- rep(1:1, 122)
correctlabel <- as.vector(c(correctlabelnormal,correctlabelALL3m))
s <- apply(correcttest, 1, function(x) coef(lm(x~correctlabel))[2])
o <- rank(s)
o <- max(o) - o + 1
res <- fgseaL(df,o,correctlabel,nperm = 2000,minSize = 1, maxSize=50000)
empty data table (0 rows) of 8 columns:   pathway,pval,padj,ES,NES,nMoreExtreme,size

I found the binary phenotype labeled group fgseaL test results below looked satisfactory.
correcttest <- data.frame(names = row.names(normal))
correcttest <- cbind(correcttest3, normal)
correcttest <- cbind(correcttest3, ALL3m)
rownames(correcttest) <- correcttest$names
correcttest$names <- NULL
correctlabelnormal <- rep(0:0, 73)
correctlabelALL3m <- rep(1:1, 122)
correctlabel <- as.vector(c(correctlabelnormal,correctlabelALL3m))
fgseaL(df,correcttest,correctlabel,nperm = 2000,minSize = 1, maxSize=50000)
       pathway        pval        padj         ES       NES nMoreExtreme  size
1: Gene.Symbol 0.003940887 0.003940887 -0.2460126 -1.180009            3 45714
                                 leadingEdge
1: AKIRIN2,LRRC20,HSPA5,HSPA5,DTWD2,ZFYVE28,

Thank you for your consideration.

It might be worth to use the package provided by the authors in github. This way you don't need to implement it again

@Llopis, Thank you for the good suggestion. I will try installing it later today..

@Llopis, Can we run rapidGSEA without the CUDA enhancemnets since I am not using a computer with a NVIDIA cpu? Thank you.

I think you could be able to do so, but I am not the maintainer neither the author of the code. You could also check the source code. But that probably is better to check and ask as another question

@Llopis, How do I install the NVIDIA c++ compiler to compile cudaGSEA? I already installed the NVIDIA CUDA Toolkit. Thank you.

@Llopis, There is a Makefile in github gravitino/cudaGSEA/src. What should I do with it if I am running Windows? Thanks

@Llopis, When I perform the base install of the  NVIDIA CUDA Toolkit on Windows 8.1, I get a diagnostic message , "NVIDIA Installer cannot continue. The NVIDIA graphics driver is not compatible with this version of Windows". What should one do here to make nvcc compile okay? Thank you.

if you have a problem don't hesitate to ask a new question on a new thread.

@Llopis, May we continue this discussion in chat?

No, if you have a new question post it as such

@Llopis, Could you please comment on my answer which I just wrote below? Thank you for your help.

It is not a question where I need clarification to answer it, so I won't comment it. Please if you have a new question post it as such. I will stop paying attention to this question an answers.

You rank the fit coefficient rather than the original score matrix. So, given a score matrix, D:
D = matrix(c(22,20,9,8,46,22,18,10,3,18,3,29,2,1,5,45,43,47,17,5,14,44,21,36), byrow=T, ncol=6)
cl = c(0,0,0,1,1,1)
s = apply(m, 1, function(x) coef(lm(x~cl))[2]) # [1]
o = rank(s)
o = max(o) - o + 1 # [2]

o is then the rank of each row.
[1] This fits each row as a linear model of cl and extracts the cl coefficient.
[2] This converts the ranking to be the same as shown in the figure. I don't know if this is important, but I would assume so given how GSEA methods tend to work.

Thank you for your excellent answer. I just accepted it. Could I pass a R language nine(9) valued phenotype label vector to where you specify cl =  c(0,0,0,1,1,1) or does cl have to be a binary phenotype label? The reason I ask this question is that we have 8 human leukemia groups and 1 healthy normal control group.

The method itself is only designed for two-group comparisons, so if you have more phenotypes you'll either have to do one comparison at a time or come up with a different method (I'm not sure what at multi-factor GSEA output would even mean in such cases).

Thank you for your reply to my comment. In the original Bioinformatics Stack Exchange question above , I have added some fgseaL source code which calculates the correlation matrix  between the Expression Data Matrix D and the R language variable , labels , which is a vector of gene phenotype labels. Could I ask you to explain the purpose of this fgseaL source code extract?

@Frank Please post that as a new question, though please read the preprint linked to from github first, since I assume that goes over the method.

I will post the new question later today. Also, I will try out your rank the linear regression fit coefficient method on our data and report the results back to you later today. Thanks.

I tested your rank the linear regression fit coefficient R language algorithm on our data with fgseaL and encountered the diagnostic message Empty data.table (0 rows) of 8 cols: pathway,pval,padj,ES,NES,nMoreExtreme. Could I ask what I did wrong? My R language code is shown at the bottom of my original question. Nonetheless , you are correct about doing two-group comparisons one at a time as I demonstrate at the bottom of my original question.

That's not my algorithm, it's just how to do that step from the paper. You'll get an empty output from the R code you posted whenever everything is filtered out. Whether that makes sense or not is a completely different question.

Thank you for your comment. What is the name of the paper and page which describes that step?

@Frank It's literally the image in your post and the paper you linked to.

I don't see any mention of linear regression coefficient  in the paper I linked to. Where did you get that idea? Also, fgseaL expects the mat input to be an matrix of doubles with the raw key value to be sorted while R language rank outputs a vector of integer rank scores. Thank you.

The "inter-class deviation" == "difference between groups" == "coefficient from a linear model". You're mixing and matching methods, you'll have to figure out how to get one to work with the other.

Is it permissible to perform the two phenotype group comparison using fgseaL without scaling on the input mat or D systematically one group at a time between each of the 8 diseased groups and the normal healthy group rather than using the linear regression coefficient R language source code as illustrated in the original question? Thank you.

@Frank Please post things like that as their own question.

Why is the fgseaL pval equal to the padj of  0.003940887 when I do the two group comparison between the normal control group and the ALL3m human leukemia disease group of subjects? Thank you.

When I do res<-fgsea(df,correcttest,nperm = 1000,minSize = 1, maxSize=50000) without the phenotype label feature why do I get the error in [.data.frame(x, order(x, na.last = na.last, decreasing = decreasing)) , undefined columns selected? Thank you.

Let us continue this discussion in chat.

Could you please comment on my answer which I just wrote below? Thank you for your help.

I finally found another answer to my question. Please read this great article in May 12 2017 BioMed Central (BMC) Bioinformatics article titled Ranking metrics in gene set enrichment analysis: do they matter?.
Also, please read this blog , Diving into Genetics and Genomics: Gene Set Enrichment Analysis (GSEA) explained.
After reading these two articles, my choice for the best rapidGSEA local ranking measure is Minimum Significant Difference (i.e., MSD), because it has the best overall false positive rate (i.e., FPR) for larger sample sizes.
Finally, it is important to realize that fgseaL's phenotype labeling can hypothetically be emulated by GSEA with either one of it's sixteen possible ranking metrics or a custom ranking metric.

This weekend, I  installed the R language packages rapidGSEA and the Broad Institute GSEA , the Nvidia CUDA Toolkit 6.5-- Custom installation option  and MINGW makefiles on my son's Windows 8.1 notebook computer.
  This step was very tricky to accomplish because I had to maneuver around the "NVIDIA Installer cannot continue. The NVIDIA graphics driver is not compatible with this version of Windows" installation error message which prevents the NVIDIA compilers from being installed correctly.
  My objective is to run rapidGSEA without the CUDA enhancements since I am not using a computer with a NVIDIA GPU.
  I did this step because rapidGSEA and Broad Institute GSEA have a significantly richer set of gene set enrichment analytics than fgseaL as shown below:
The GSEA method takes five default arguments and three optional arguments
GSEA <- function(exprsData, labelList, geneSets, numPermutations, metricString,
                 dumpFileName="", checkInput=TRUE, doublePrecision=FALSE) {...}
exprsData, labelList and geneSets refer to the data obtained in the previous section. numPermutations denotes the number of permutations in the resampling test, metricString denotes the local ranking measure (one of the following):
naive_diff_of_classes
naive_ratio_of_classes
naive_log2_ratio_of_classes
stable_diff_of_classes
stable_ratio_of_classes
stable_log2_ratio_of_classes
onepass_signal2noise
onepass_t_test
twopass_signal2noise
twopass_t_test
stable_signal2noise
stable_t_test
overkill_signal2noise
overkill_t_test

Glad you got it working. Not sure what else I'm supposed to comment on here.

@Devon Ryan, If metricString denotes the local ranking measure we are requesting that rapidGSEA use when it is invoked, which metricString value should I choose so that rapidGSEA operates in a manner most closely resembling fgseaL? Thank you.

Don't try to make one work like the other. If you start making your own method then you have to get it through peer review. It doesn't appear you have the background necessary to do that. Use rapidGSEA or fgseaL, stop trying to merge them.


text:
Get canonical transcript from UCSC
I am using the following command to get all refseq genes from UCSC:
/usr/bin/mysql --user=genomep --password=password --host=genome-mysql.cse.ucsc.edu \
    -A -D hg38 -e 'select concat(t.name, ".", i.version) name, \
    k.locusLinkId as "EntrezId", t.chrom, t.strand, t.txStart, \
    t.txEnd, t.cdsStart, t.cdsEnd, t.exonCount, t.exonStarts, \
    t.exonEnds, t.score, t.name2 from refGene t join hgFixed.gbCdnaInfo i \
    on t.name = i.acc join hgFixed.refLink k on t.name = k.mrnaAcc'

That returns data in the following format (showing the 1st 5 lines):
+-------------+-----------+-------+--------+---------+-------+----------+--------+-----------+--------------------+--------------------+-------+-----------+
| name        | EntrezId  | chrom | strand | txStart | txEnd | cdsStart | cdsEnd | exonCount | exonStarts         | exonEnds           | score | name2     |
+-------------+-----------+-------+--------+---------+-------+----------+--------+-----------+--------------------+--------------------+-------+-----------+
| NR_046018.2 | 100287102 | chr1  | +      |   11873 | 14409 |    14409 |  14409 |         3 | 11873,12612,13220, | 12227,12721,14409, |     0 | DDX11L1   |
| NR_106918.1 | 102466751 | chr1  | -      |   17368 | 17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-1 |
| NR_107062.1 | 102465909 | chr1  | -      |   17368 | 17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-2 |
| NR_107063.1 | 102465910 | chr1  | -      |   17368 | 17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-3 |
| NR_128720.1 | 103504738 | chr1  | -      |   17368 | 17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-4 |
+-------------+-----------+-------+--------+---------+-------+----------+--------+-----------+--------------------+--------------------+-------+-----------+

I also want to find the accession of the canonical transcript of each of the genes returned by the command above. Those seem to be stored in the knownCanonical table:
$ /usr/bin/mysql --user=genomep --password=password --host=genome-mysql.cse.ucsc.edu -A -D hg38 -e 'select * from knownCanonical limit 5'
+-------+------------+-----------+-----------+------------+--------------------+
| chrom | chromStart | chromEnd  | clusterId | transcript | protein            |
+-------+------------+-----------+-----------+------------+--------------------+
| chrX  |  100628669 | 100636806 |         1 | uc004ega.3 | ENSG00000000003.14 |
| chrX  |  100584801 | 100599885 |         2 | uc004efy.5 | ENSG00000000005.5  |
| chr20 |   50934866 |  50958550 |         3 | uc002xvw.2 | ENSG00000000419.12 |
| chr1  |  169853073 | 169893959 |         4 | uc001ggs.5 | ENSG00000000457.13 |
| chr1  |  169795048 | 169854080 |         5 | uc001ggp.4 | ENSG00000000460.16 |
+-------+------------+-----------+-----------+------------+--------------------+

However, there seems to be no obvious way to link the knownCanonical table to the refGene, hgFixed.gbCdnaInfo and hgFixed.refLink tables used above.
So, how can I modify my 1st query (or, if necessary write a new one) so that my results also include the accession of the gene's canonical transcript?

You can link them with the kgXref table, since kgXref.refseq == refGene.name and kgXref.kgID == knownCanonical.transcript. Since it seems that knownCanonical.transcript is what you want anyway, you don't even need to join on it:
mysql --user=genomep --password=password --host=genome-mysql.cse.ucsc.edu \
-A -D hg38 -e 'select concat(t.name, ".", i.version) name, x.kgID, \
k.locusLinkId as "EntrezId", t.chrom, t.strand, t.txStart, \
t.txEnd, t.cdsStart, t.cdsEnd, t.exonCount, t.exonStarts, \
t.exonEnds, t.score, t.name2 from refGene t join (hgFixed.gbCdnaInfo i, hgFixed.refLink k, kgXref x) \
on (t.name = i.acc and t.name = k.mrnaAcc and t.name = x.refseq) limit 10'

The output is then:
+-------------+------------+-----------+-------+--------+---------+--------+----------+--------+-----------+--------------------+--------------------+-------+------------+
| name        | kgID       | EntrezId  | chrom | strand | txStart | txEnd  | cdsStart | cdsEnd | exonCount | exonStarts         | exonEnds           | score | name2      |
+-------------+------------+-----------+-------+--------+---------+--------+----------+--------+-----------+--------------------+--------------------+-------+------------+
| NR_106918.1 | uc031tla.1 | 102466751 | chr1  | -      |   17368 |  17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-1  |
| NR_107062.1 | uc031tlm.1 | 102465909 | chr1  | -      |   17368 |  17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-2  |
| NR_107063.1 | uc032cta.1 | 102465910 | chr1  | -      |   17368 |  17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-3  |
| NR_128720.1 | uc032dmn.1 | 103504738 | chr1  | -      |   17368 |  17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-4  |
| NR_036051.1 | uc031tlb.1 | 100302278 | chr1  | +      |   30365 |  30503 |    30503 |  30503 |         1 | 30365,             | 30503,             |     0 | MIR1302-2  |
| NR_036266.1 | uc033cjs.1 | 100422831 | chr1  | +      |   30365 |  30503 |    30503 |  30503 |         1 | 30365,             | 30503,             |     0 | MIR1302-9  |
| NR_036267.1 | uc032csz.1 | 100422834 | chr1  | +      |   30365 |  30503 |    30503 |  30503 |         1 | 30365,             | 30503,             |     0 | MIR1302-10 |
| NR_036268.1 | uc032hiw.1 | 100422919 | chr1  | +      |   30365 |  30503 |    30503 |  30503 |         1 | 30365,             | 30503,             |     0 | MIR1302-11 |
| NR_026822.1 | uc001aak.4 |    654835 | chr1  | -      |   34610 |  36081 |    36081 |  36081 |         3 | 34610,35276,35720, | 35174,35481,36081, |     0 | FAM138C    |
| NR_106918.1 | uc031tla.1 | 102466751 | chr1  | -      |  187890 | 187958 |   187958 | 187958 |         1 | 187890,            | 187958,            |     0 | MIR6859-1  |
+-------------+------------+-----------+-------+--------+---------+--------+----------+--------+-----------+--------------------+--------------------+-------+------------+

If you want the protein ID then join knownCanonical on x.kgID.


text:
Secretome and membrane receptor profiles
I am looking the secretome profile and the membrane receptor profile for a given cell type.
In my specific case, this should be the secretome and outer membrane receptor profiles of dorsal root ganglion.
What I've done in the past is taken proteomics datasets or RNAseq datasets and used this as a reference cell surface proteome or secretome.
I wasn't sure what the best way to filter secreted proteins or membrane receptors were, so I used the Human Protein Atlas and this turned out to be a big headache. A possible alternative method might be to use DAVID, but there must be a more efficient method via the command line.
My questions are these:

Is there a "consensus" database containing protein expression
profiles for different tissue types in human?
If not, what is the best way to filter out the secreted and transmembrane proteins by gene symbol?

Is there a "consensus" database containing protein expression profiles for different tissue types in human?

https://www.ebi.ac.uk/gxa/home

If not, what is the best way to filter out the secreted and
  transmembrane proteins by gene symbol?

Personally I like to run a local copy of SignalP, as it will yield consistent results. However bioinformatic approaches tend not to capture highly unusual transmembrane proteins or membrane associations (e.g.: only one leaflet etc.). http://www.cbs.dtu.dk/services/SignalP/
Depending on the organisms, and scope of your analysis, further possibility could be to visit supplemental tables of recent publications containing somewhat supervised computational predictions on secretion, e.g.: http://science.sciencemag.org/content/356/6340/eaal3321

I am looking the secretome profile and the membrane receptor profile for a given cell type.

The above links don't inform about cell-type specific retention of potentially secreted proteins. - but their combination will yield cell-type specific RNA and protein levels of transmembrane and secreted proteins.

I don't think there is a database with expression data for your particular cell type. I am afraid your "headache" approach doesn't sound so bad...
So get expression data e.g., from GEO for your cell type, and define which genes/proteins you consider as "expressed" (e.g., minimal RPKM value).
Then, you can use GO annotation to find your genes involved with membrane or secretion (or what ever). You can use biomaRt in R for this.
For example, you want to know all "membrane" genes.
> library(biomaRt)
>
> ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
>
> membrane.genes <- getBM(attributes=c('hgnc_symbol',
> 'ensembl_transcript_id', 'go_id'),
> filters = 'go', values = 'GO:0016020', mart = ensembl)

And now you want to know it only from genes that are expressed in your ganglion cells (with dummy example of two genes).
> gene.list <- c("KIR2DL3", "GPR1")
>
> membrane.genes.list <- membrane.genes[membrane.genes$hgnc_symbol ==
> gene.list,]
>
> membrane.only <- membrane.genes.list[membrane.genes.list$go_id ==
> 'GO:0016020',]

Note that this will only retrieve genes directly annotated to GO:0016020 (membrane). For example, the gene PLSCR5 is not in the list despite this being annotated to plasma membrane (GO:0005886). Also, PM may be a better choice for exploration of cell-type specific surface receptors.

I try to give an example of how to use biomaRt for GO annotation. GO annotation is not a holy grail, but just a tool that can be helpful in your research sometimes.

Agreed, it's a good answer, just pointing out this for other people coming by this question and using biomaRt for other GO queries.

You can get expression data from databases and for membrane proteins, run a transmembrane protein prediction method like CCTOP (http://cctop.enzim.ttk.mta.hu/).


text:
How to convert BED to GFF3
I would like to convert a BED format to GFF3.
The only useful tool that I could find via a google search seems to be Galaxy, and I do not feel very comfortable with online tools, plus the webserver is currenlty under maintenance.
Does anyone knows about a command-line tool that can handle this conversion?
Edit: here are some lines of my BED file:
$ head -4 last_minion-r7_sort.bed
211000022278137 175 211 8e5d0959-7cdb-49cf-9298-94ed3b2aedb5_Basecall_2D_000_2d 42  +
211000022279134 0   503 e8a9c6b8-bad2-4a7e-97d8-ca4acb34ff70_Basecall_2D_000_2d 69  -
211000022279134 24  353 e258783d-95a3-41f5-9ad5-bb12311dbaf4_Basecall_2D_000_2d 45  -
211000022279134 114 429 26601afb-581a-41df-b42b-b366148ea06f_Basecall_2D_000_2d 100 -

The bed file thus has 6 columns as for now: chromosome, start coordinate, end coordinate, read name, score, strand. This file was obtained from conversion of MAF format (as output of alignment of RNA-seq reads to reference genome, using LAST) converted to SAM using maf-convert, then to BAM using samtools, finally to BED using bedtools.
The aim of my conversion is basically to convert SAM -> GTF, for post-processing. Since there is no straightforward way to do this, I am going through steps, the only way to do this in my knowledge is : SAM -> BAM -> BED -> GFF3 -> GTF but for now I am stuck in the BED -> GFF3 part.

There isn't necessarily a very meaningful way to do that. Can you post a few lines from your BED file? Also, what's the end goal?

Yes, it would be helpful to see some of the BED. Is it BED12 or BED6, what does it represent? In GFF3 who have to say what your intervals represent.  No reason it can't be done though. Oh, and are you happy with code?

@DevonRyan I have edited the question, hopefully answers your comments.

@AminaEchchiki While I could provide code (or you could just use the tool in Galaxy), I don't think the information from the SAM file is being accurately transmitted to BED in a way that's needed to produce a useful GFF3 or GTF file. I strongly recommend you write a tool to directly produce a GTF/GFF from either SAM or MAF (presuming none exist).

if this is an organism that has splicing, your convsion to bed is unlikely to be accurate. I will post BAM -> GTF code and BED to GFF3, but let me know if your reads are spliced/

@IanSudbery yes my reads are spliced (D. melanogaster)

Right, in which case you should generate your bed files with bamtobed --bed12, but really you should go direct from BAM to GFF/GTF.

@IanSudbery thanks for the tip, will do that

Galaxy has API and API-consuming libraries (such as BioBlend) that will allow you to interactively script against it without opening the graphical interface at all.
However you can also take almost any tool out of Galaxy and use it independently since everything is open source. The converter you mentioned is available as a Python script here and the tool 'wrapper' which you can use to understand how to invoke the python script is next to it.

To answer the question as asked,  for people googling.
For BED6, in python:
#contigs.tsv contians chromosome names and lengths in two columns
for line in open("contigs.tsv"):
   fields = line.strip().split("\t")
   print fields[0], ".", "contig","1",str(fields[1]), ".", "+", ".", "ID=%s" % fields[0]

for line in open("my_bed_file.bed"):
   fields = line.strip().split("\t")

   # note: BED is 0-based, half-open, GFF is 1-based, closed
   start = str(int(fields[1]) + 1)
   print fields[0], "bed", "interval", fields[1], fields[2], fields[4], fields[5], ".", "ID=%s;parent=%s" % (fields[3], fields[0])

For bed12, in python:
#contigs.tsv contians chromosome names and lengths in two columns
for line in open("contigs.tsv"):
   fields = line.strip().split("\t")
   print fields[0], ".", "contig","1",str(fields[1]), ".", "+", ".", "ID=%s" % fields[0]

for line in open("my_bed12.bed"):

   fields = line.strip().split("\t")
   contig = fields[0]
   # note: BED is 0-based, half-open, GFF is 1-based, closed
   start= int(fields[1]) + 1)
   end=fields[2]
   name=fields[3]
   score=fields[4]
   strand=fields[5]
   print contig, "bed", "interval", str(start), end, score, strand, ".", "ID=%s;parent=%s" % (name, contig)

   block_starts = map(int,fields[11].split(","))
   block_sizes = map(int, fields[10].split(","))

   for (block, (bstart, blen)) in enumerate(zip(block_starts, block_sizes)):
      bend = start + bstart + blen
      print contig, "bed", "block", str(start + bstart), str(bend), score, strand, ".", "ID=%s_%i;parent=%s" %(name, block, name)

To convert BAM to GTF, which is the best way to get a file to compare with cuffcompare:
import pysam

bamfile=pysam.AlignmentFile("my_bam_file.bam")

for alignment in bamfile.fetch():

    if alignment.is_unmapped:
        continue

    contig = bamfile.get_reference_name(alignment.reference_id)
    name = alignment.query_name

    if alignment.is_reverse:
        strand = "-"
    else:
        strand = "+"

    for start, end in alignment.getblocks():
        # note: BED is 0-based, half-open, GFF is 1-based, closed
        print contig, "BAM", "exon", str(start + 1), str(end), "0", strand, ".", 'gene_id "%s"; transcript_id "%s";' % (name, name)

NB, this will also work with a SAM file as long as it is headered.

As a note for future googlers, if your bam file comes from a pipeline producing unmapped BAMs, the code as it is will return error ValueError: file header is empty (mode='rb') - is it SAM/BAM format?. As a workaround it is necessary to add options check_header=False, check_sq=False while reading the bam.
@IanSudbery  might be useful to add it as a comment to your code

@Amina Echchiki, this error is not because your BAM contains unmapped reads, but because your BAM has no header. BAMs without headers are not valid BAM files. I am surprised fetching the contig name worked if the file has no header, even if you set check_header=False, check_sq=False because it uses the header to do so.

Bioconductor makes this so easy. It does the coordinate conversion on import.
library(rtracklayer)

## import the bed file
bed.ranges <- import.bed('regions.bed')

## export as a gff3 file
export.gff3(bed.ranges,'regions.gff3')

And people wonder why R is so popular for bioinformatics...
Also if needed, you could go straight from BAM file to gff3 in R as well.

As long as you have enough memory to hold the whole bed in memory at once (not a given if it comes from a BAM) this is indeed an easy option.

I recently developed bed2gff to quickly convert .bed files to a gff3 format, a tool written in Rust. Could be of help here!

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.

If anybody needs a simple bed to gff converter, simple in the sense that only takes a bed with a single type of feature, lets say peaks and translate it to GFF I recently wrote this, is Bash and as far as I know all the tools it uses should be present in most Linux / Mac operating systems.
https://github.com/UBMI-IFC/bed2gff/tree/main
It cannot deal, however with multi featured Bed files as it ignores everything from the 7th column and beyond.

You could also use agat_convert_bed2gff. agat is generally a very useful tool when needing to convert from / to / between gff formats. See here how to use it to convert from bed to gff. https://agat.readthedocs.io/en/latest/tools/agat_convert_bed2gff.html


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

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

The output python:
$ python

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

The command line output from pip isn’t relevant for the question, I’ve thus removed. pip --version would have been more relevant. Also, be sure to read the formatting help for questions, especially the part about preformatted text (which applies to console output).

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

But:
Your Python installation is a bit screwed up: there are Python 2 and Python 3, which are unfortunately incompatible.
On your system python and pip seem to be aliases for Python 3, which I’d strongly recommend against (lots of tools will break). Instead, they should alias Python 2. Python 3 should be invoked via python3/pip3.

That's a pretty sweeping statement. The command python defaults to python3 on my two Arch systems and various Ubuntu boxes we use at work and I've never had issues with it. What makes you say that using python2 as the default is better? And on what systems?

@terdon It's from a PEP recommendation. And failing to do so, while not catastrophic, definitely breaks many scripts that assume python = Python 2. Not a fan of Arch's/Ubuntu's decision to ignore that unnecessarily.

Well, I'd argue that using a shebang line calling python instead of python[23] is a very bad design decision since the python devs, in their infinite wisdom, made the two versions incompatible. Either way, your advice is dangerous since, presumably, the OP has a working system with python3 as the default and changing that will break lots of tools while leaving it as is will not.

@terdon You realise that there are Python scripts that predate the existence of Python 3, right? The all will, by necessity, have python in their shebang line. And fixing the system shouldn't break anything but, yes, it will require reinstalling all modules.


text:
Major variability in speed of BLAST between proteins
I'm running blast+ searches for a list of proteins using specific parameters.
This is my search
 blastp -query protein.fasta -remote -db refseq_protein -out
protein_refseq_blast.txt -outfmt '6 qseqid sseqid qstart qend length evalue' -
entrez_query 'txid2 [ORGN] or txid4751 [ORGN]' -evalue 200000 -max_target_seqs 10000

For some proteins, this runs moderately quickly, ~5-10 mins. For other proteins, this command never stops running, even after an hour or so. I'm wondering if this variability in speed is due to any part of this choice in parameters, is simply variability in speed for the BLAST servers. Is it common for BLAST searches to take this long? Is there a way to run all of them in parallel using these parameters if this is the case?
Thanks.

FWIW, I remember experiencing very long waiting times when running remote blast from Biopython.

You have your evalue set incredibly high -evalue 200000, which might affect the time it takes to return a query.

The variability could be simply because the proteins are very different sizes but even 5-10 mins is very long for a single protein, as fugu says, your evalue is crazy! Probably reducing the evalue to a 'normal value' (<1) will resolve this

Try changing you evalue threshold to something like 0.001.
If you blasted a sequence against a database containing 10,000 sequences and got a hit with an evalue of 1, this would mean that, given the size of your database and length of protein, you would expect to find a match as good as (or better than) the hit returned once purely by chance. If the evalue threshold in your query was set to 0.01 (for example) then this hit would not be returned.
If you got a hit with an evalue of 200,000 this would mean you would expect to see 200,000 matches with a similar score simply by chance. If your evalue threshold were set to 200000, then all of the hits with an evalue <=200000 would be considered significant. If you set it to 0.001 then only those hits with evalues <= 0.001 would be considered significant.
This means that (depending on the size of the data base) almost every sequence in your query set will be considered as significant given your threshold, which will slow down your search hugely.
See the the section on E values here

The evalue associated with a blast HSP (High-scoring Segment Pair) is a measure of how often you would expect to find such an HSP in a database of this size by pure chance. For instance, if you blast the sequence ACTG against the entirety of the NCBI's nr database, you would expect several hundred thousand hits by pure chance. That is, several hundred thousand sequences that can be aligned to your query but are in no way related to it (not actual homologs). Those hits don't actually mean anything, so you use the evalue to limit the results to those that are unlikely to be found by chance and, therefore, more likely to have biological meaning if hits are found.
So, setting the evalue to 200000 means you are interested even in HSPs whose characteristics (score, length, etc) mean you would expect to find <=200000 HSPs with such characteristics by pure chance. This evalue is pointless.
Because your e-value is so high, BLAST will report absolutely anything and everything it finds, and the vast majority of your results will be meaningless and irrelevant. If you want a permissive e-value, you could use 1 or even 10 but 200000 is way too much. You usually only want negative e-values, cases where you wouldn't expect to find even one such HSP by chance. For highly dissimilar sequences, I would recommend you use 0.1 or something of that scale.
Changing the evalue will speed things up considerably. That said, you should always expect different times for different queries since the time it takes to find hits will always depend on the number of hits that can be found. If your sequence can be aligned to many of the sequences in the database, it will take longer to fetch all results.


text:
Software recommendations - DNA composition
I'm looking to generate data based on the DNA composition of a region of my genomes (data is incomplete genomes from HiSeq runs in fasta format). I'm looking for software which will give me sliding windows for GC content, GC skew, codon bias etc. I would like to use something which I can use with or plug into R.

There are more R packages available that calculate GC content, for example Ape's GC.content() function.
For example:
> library(ape)
> data(woodmouse)
> GC.content(woodmouse)
[1] 0.3873347

With a sliding window is in Biostrings package.
> library(Biostrings)
> DNA <- DNAString("ACTGAAACCGTGGCAGTTTGAC")
> letterFrequencyInSlidingView(DNA, view.width=4,letters="CG")

For codon usage bias, I think cubfits or sscu packages might be useful.

Could you please [edit] your question and give an example of how the OP can use this package, what kind of results it would return, the type of input data etc? As it stands, this is more of a comment than an answer.

According to ape gc.content will only give a percentage. I'm looking for the overall distribution.

The definition of GC content is a percentage.


text:
IGV faster loading remote bigWig files
I noticed remote bigWig files (http URL) take sometimes about 1 minute to load on IGV. Once they loaded, jumping from one region to another can also take several seconds.
My understanding is that bigWig files have an internal index for different chunks of different sizes. Can the index be modified so that the speed of loading on IGV is improved? Does this need to happen at the point of creation of the bigWig file, or can the index of an existing bigWig be updated to speed up loading?

You would need to do this when creating the bigWig file. Note that I'm not aware of anything that allows manually setting the various zoom levels (libBigWig/pyBigWig will allow you to set a max number, but that's it).
This then raises the question of whether IGV is using an appropriate zoom level to begin with. I've never looked at how it chooses which of the zoom levels to use, but I know it uses them since if you don't have any it will fail to load bigWig files.
Honestly, though, if it's taking a minute to jump to a new region then I suspect that the network is just slow or the remote server is overtaxed. Changing the internal indexing isn't going to fix that.


text:
Merging bed records based on name
I generated a file starting with the following bed lines:
$ head -6 /tmp/bed_with_gene_ids.bed
I   3746    3909    "WBGene00023193"    .   -
I   3746    3909    "WBGene00023193"    .   -
I   4118    4220    "WBGene00022277"    .   -
I   4118    4358    "WBGene00022277"    .   -
I   4118    10230   "WBGene00022277"    .   -
I   4220    4223    "WBGene00022277"    .   -

I would like to merge them based on the name field (the 4-th column), taking the min for the start and the max for the end. Other fields are expected to be the same for all records having the same name.
Expected result:
I   3746    3909    "WBGene00023193"    .   -
I   4118    10230   "WBGene00022277"    .   -

I found a potential solution based on bedtools groupby here: https://www.biostars.org/p/145751/#145775

Sample data:

cat genes.bed
chr14    49894259    49895806    ENSMUST00000053290    0.000000    ...
chr14    49894873    49894876    ENSMUST00000053290    0.000000    ...
chr14    49894876    49895800    ENSMUST00000053291    0.000000    ...
chr14    49895797    49895800    ENSMUST00000053291    0.000000    ...
chr14    49901908    49901941    ENSMUST00000053291    0.000000    ...

Example output:

sort -k4,4 genes.bed \
| groupBy -g 1,4 -c 4,2,3 -o count,min,max \
| awk -v OFS='\t' '{print $1, $4, $5, $2, $3}'

chr14    49894259    49895806    ENSMUST00000053290    2
chr14    49894876    49901941    ENSMUST00000053291    3

However:

I don't understand the groupBy behaviour (Why -g 1,4 and not just -g 4?, Why -c 4,2,3 in this order and then rearrange things using awk?)

This code doesn't work for me.

Here is what happens when I try the solution given above:
$ head -3 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 1,4 -c 4,2,3 -o count,min,max | awk -v OFS='\t' '{print $1, $4, $5, $2, $3}'
3           3746    4220

Here are attempt based on what I thought could work according to the documentation:
$ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 4 -c 1,2,3,4,5,6 -o first,min,max,distinct,first,first
I   3746    10230   "WBGene00022277","WBGene00023193"   .   -

$ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 4 -c 1,2,3,4,5,6 -o first,min,max,last,first,first
I   3746    10230   "WBGene00022277"    .   -

$ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 4 -c 1,2,3,5,6 -o first,min,max,first,first
I   3746    10230   .   -

I don't get why when I group based on the 4-th column, for which I have two distinct values, I cannot obtain two lines in the resulting output.
I understand based on the comments on the documentation page that the documentation is not up-to-date. In particular, there is a -full option that is needed if one wants all fields to be outputted. Re-reading the solution mentioned above, I think I now understand the reason for the multiple columns for the -g option and for the awk rearrangement. Hence the following attempt.
$ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 1,4,5,6 -c 2,3 -o min,max -full
    I   3746    3909    "WBGene00023193"    .   -   3746    10230

But this still doesn't give me two lines.
Are there other tools that could do what I want efficiently?

Edit: Solution
According to this answer, the problem with bedtools is that there is a bug in the latest release (2.26.0 as of august 2017). In order to have a functional bedtools groupby, one needs to get the development version from github.
With the github version of bedtools, I can now get the expected result as follows:
$ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 1,4,5,6 -c 2,3 -o min,max | awk -v OFS="\t" '{print $1,$5,$6,$2,$3,$4}'
I   3746    3909    "WBGene00023193"    .   -
I   4118    10230   "WBGene00022277"    .   -

I include fields 1, 5 and 6 in -g (besides field 4) in order to have them printed out. In my bed file, they should be the same for a given value of field 4. The awk part is needed because one has apparently not total control on the output order: the -g fields come before the -c fields.
Note (added 23/12/2020): As specified in bedtools groupby documentation, the input bed file should be sorted according to (chromosome (1) and) name (4) field(s) for this to work properly. Otherwise, duplicate (i.e. non-merged) names may appear in the output.

What do you want to do with the score and strand fields if they're different between lines, or does that never happen?

Actually, I don't care about the score field, and I would ideally set it to "." if it is not already the case. I cannot guarantee that the strand field will always be the same, but since these bed lines come from transcript annotations whose gene_id I've put in the name field, I suppose it will generally be true that for a same name, the strand will be the same. I should check this, though.

Although you don't mention it, I'm guessing you're using bedtools v2.26.0. Version 2.26.0 of groupBy has a bug in it, which you've encountered (it was fixed shortly after release, so you'll either have to use a version before the bug was introduced, or compile the current source code yourself from https://github.com/arq5x/bedtools2)
v2.26.0:
local10:~/Documents/tmp$ cat asdf.bed
I 3746 3909 WBGene00023193 . -
I 3746 3909 WBGene00023193 . -
I 4118 4220 WBGene00022277 . -
I 4118 4358 WBGene00022277 . -
I 4118 10230 WBGene00022277 . -
I 4220 4223 WBGene00022277 . -
local10:~/Documents/tmp$ groupBy -i asdf.bed -g 4 -c 2,3 -o min,max
3746    10230

v2.26.0-125-g52db654 (I.E. compiling the source code from github):
local10:~/Documents/tmp$ bedtools2/bin/groupBy -i asdf.bed -g 4 -c 2,3 -o min,max
WBGene00023193  3746    3909
WBGene00022277  4118    10230

To answer your questions:
1) You might notice that my output above gives the grouped columns first; you'll have to reorder the output via awk in order to get it back in order. As for why they chose to group on both columns 1 and 4: if you have the same name on multiple chromosomes, you may want to treat them as separate features.
2) Version differences, as stated in the first part of my answer.

To actually merge the file:
Make sure to run this with a version other than v2.26.0 (As Devon Ryan writes in the comments, you may want to add column 6 to -g to make it strand-specific):
./bedtools2/bin/groupBy -i asdf.bed -g 1,4 -c 2,3,5,6 -o min,max,first,first \
     | awk -v OFS='\t' '{print $1, $3, $4, $2, $5, $6}'
I   3746    3909    WBGene00023193  .   -
I   4118    10230   WBGene00022277  .   -

If you include 6 in -g 1,4 then you benefit by not merging genes on different strands. UCSC sometimes has those and they really aren't the same gene and shouldn't be merged together. You don't need 1 in -c, or 6 if you add it to -g.

You could do this with the CGAT toolkit:
cgat bed2bed --method=merge --merge-by-name -I bed_with_gene_ids.bed
Installing such a massive package might be overkill for this task though.

It happens that cgat is already installed on my computer (though I forgot for what purpose). I tried the command you suggest, and I end up with a duplicate of I 3746 3909 "WBGene00023193" . - . Granted, there were duplicate lines in the original bed. But is this behaviour expected?

Also, If I run this on the whole file and not just the first 6 lines, after a while, the program fails on TypeError: '<' not supported between instances of 'Bed' and 'Bed'. I'm upgrading cgat to see if the error persists.

I reported the issues here: https://github.com/CGATOxford/cgat/issues/347

Your first problem is not as far as I am aware the intended behaviour. And I'm pretty sure the second is not intended. I suggest you file a bug report.

We crossed posts!

You can do this easily with Hail.  Hail primarily uses BED files to annotate genetic datasets (see the last annotate_variants_table example), but you can manipulate BED files using Hail's general facilities for manipulating delimited text files.  For example:
$ cat genes.bed
I   3746    3909    "WBGene00023193"    .   -
I   3746    3909    "WBGene00023193"    .   -
I   4118    4220    "WBGene00022277"    .   -
I   4118    4358    "WBGene00022277"    .   -
I   4118    10230   "WBGene00022277"    .   -
I   4220    4223    "WBGene00022277"    .   -

The Hail script (python code):
from hail import *
hc = HailContext()
(hc
 .import_table('genes.bed', impute=True, no_header=True)
 .aggregate_by_key('f0 = f0, f3 = f3',
    'f1 = f1.min(), f2 = f2.max(), f4 = ".", f5 = "-"')
 .select(['f0', 'f1', 'f2', 'f3', 'f4', 'f5'])
 .export('genes_merged.bed', header=False))

The result:
$ cat genes_merged.bed
I   3746    3909    WBGene00023193  .   -
I   4118    10230   WBGene00022277  .   -

I aggregate over chrom and name so this solution won't merge entries on different chromosomes. The select is necessary to reorder the fields because aggregate_by_key places the keys being aggregated over first.
Disclosure: I work on Hail.

$ cut -f4-6 in.bed | sed 's/\t/_/g' | sort | uniq | awk -F'_' '{ system("grep "$1" in.bed | bedops --merge - "); print $0; }' | paste -d "\t" - - | sed 's/_/\t/g' | sort-bed - > answer.bed

Given your sample input:
$ more in.bed
I   3746    3909    "WBGene00023193"    .   -
I   3746    3909    "WBGene00023193"    .   -
I   4118    4220    "WBGene00022277"    .   -
I   4118    4358    "WBGene00022277"    .   -
I   4118    10230   "WBGene00022277"    .   -
I   4220    4223    "WBGene00022277"    .   -

The answer.bed file:
$ more answer.bed
I   3746    3909    "WBGene00023193"    .   -
I   4118    10230   "WBGene00022277"    .   -

Sorting with sort-bed is useful at the end, so that you can pipe it or work with it with other BEDOPS tools, or other tools that now accept sorted BED input.
Streaming is a pretty efficient way to do things, generally.

How this works
Here's the pipeline again:
$ cut -f4-6 in.bed | sed 's/\t/_/g' | sort | uniq | awk -F'_' '{ system("grep "$1" in.bed | bedops --merge - "); print $0; }' | paste -d "\t" - - | sed 's/_/\t/g' | sort-bed - > answer.bed

We start by cutting columns 4 through 6 (id, score and strand), replacing tabs with underscores, sorting and removing duplicates:
cut -f4-6 in.bed | sed 's/\t/_/g' | sort | uniq

What we get out of this is a sorted list of "needles" — one for each ID-score-strand combination: an ID-needle — that we can use to grep or filter the original BED file.
This list is piped to awk which, for each ID-needle, runs grep against the original BED file and pipes the subset to bedops --merge -, which merges overlapping intervals.
Note that merging only works for overlapping intervals. Merging is not necessarily the same as returning a min-max pair, and this pipeline will break if there are intervals that do not overlap. But you could modify the awk statement to process the input intervals and return the minimum and maximum interval coordinates, if that is really what you want, by tracking the min and max values over all intervals that come into awk, and printing a final interval with an END block.
The system command prints the merged interval on one line. The following print $0 statement prints the needle on the next line:
awk -F'_' '{ system("grep "$1" in.bed | bedops --merge - "); print $0; }'

We take each pair of alternating lines and re-linearize them with paste. This result now contains four columns: the three columns of each merged interval, and the ID-needle.
We then use sed to replace underscores with tabs, so that we turn the ID-needle back into three, tab-separated ID-score-strand columns:
paste -d "\t" - - | sed 's/_/\t/g'

The output is now a six-column BED file, but it is ordered by the sort order we applied on ID-needles further up the pipeline, which we don't want. What we really want is BED that is sorted per BEDOPS sort-bed, so that we can do more set operations and get a correct result. So we pipe this to sort-bed - to write a sorted file to answer.bed:
sort-bed - > answer.bed

Thanks for the answer, it works, and I think I understood how. Maybe some explanations about the different steps could be useful.

If you're 100% sure that everything but the start and end positions will be the same for all lines sharing a name, you could just do it yourself. For example, in Perl:
$ perl -lane '$start{$F[3]}||=$F[1];
              if( $F[1] < $start{$F[3]} ){
                 $start{$F[3]} = $F[1]
              }
              if( $F[2] > $end{$F[3]} ){
                 $end{$F[3]} = $F[2]
              }
              $chr{$F[3]} = $F[0];
              $rest{$F[3]} = join "\t", @F[4,$#F];
              END{
                 foreach $n (keys %chr){
                  print "$chr{$n}\t$start{$n}\t$end{$n}\t$n\t$rest{$n}"
                 }
              }' file.bed
I   3746    3909    "WBGene00023193"    .   -
I   4118    10230   "WBGene00022277"    .   -

I was hoping that an efficient tool already existed and would avoid me re-inventing the wheel in a slow scripting language.

@bli absolutely, that makes much more sense. I just figured this is simple enough, so I may as well give a scripting solution. But yes, this will be slow and is also very naive so it will break if your files are slightly different.


text:
Melt p-values for CpG sites mapping to the same gene
I have some data I am working with, and I am curious if I am able to combine p-values from a paired t-test for CpG sites in the genome using Fisher's Method to get one p-value for each unique gene. Linked here is the Wikipedia page for Fisher's Method. I understand that an assumption of the method used is that each individual p-value being combined must be independent. I'm relatively new to biostatistics, so I'm curious if using CpGs from the same gene would violate this assumption.

Why do you want to melt CpG values? Do you want to know if an specific gene/region is more phosphorylated than others? What are your hypothesis (you tell about a t-test for the CpG sites? Are you testing if they are phosphorylated vs not or more than the rest of the genome or...?

In a nutshell, your intuition is correct: they are not independent.

Methylation levels have high local correlation, so Fisher's method would be problematic. Having said that, you have no reason to use Fisher's method after a paired t-test. A paired t-test will give you a single p-value per gene, which is what you want. Do be sure to only include CpG with some minimal coverage in both group.

What is typically done in methylation analysis is to assess the islands of methylations.
Check this workflow, in the section linked it uses some predefined islands for instance. I am no expert on this area but you could asses if the islands or certain regions are more methylated than expected.


text:
Given a Genomic Ranges of SNPs, how to inject these SNPs in genome via BSGenome?
Let's say I have the genome hg19 loaded into R via BSGenome
library("BSgenome")
hg19genome = getBSgenome('BSgenome.Hsapiens.UCSC.hg19', masked=FALSE)

I then have a list of SNPs loaded as a GRanges object, gr
library(GenomicRanges)

 > gr
 GRanges object with 212 ranges and 3 metadata columns:
        seqnames               ranges strand |     width     REF    ALT
           <Rle>          v <IRanges>  <Rle> | <numeric>
    [1]        1 [86099032, 86099032]      * |         1     C      T
    [2]        1 [86099033, 86099033]      * |         1     C      A
    [3]        1 [86099199, 86099199]      * |         1     G      A
  ....

Is there a straightforward way to inject these SNPs into hg19genome?

By inject you mean that the hg19genome should have the alternative base instead of the reference base for those positions?

@Llopis Yes, that's what I mean

This seems relatively complicated given the structure of a BSGenome object.
The creator of the package answered this question previously on the Bioconductor support forums:
https://support.bioconductor.org/p/86665/#86757

We don't provide an easy way to inject arbitrary SNPs in an arbitrary
  BSgenome at the moment. However, it should not be too hard to forge
  the BSgenome ... package. First you would need to
  compute the sequences of the mutated chromosomes (you can
  use replaceLetterAt for this), then write them to a 2bit file (put
  them in a DNAStringSet object and call rtracklayer::export on it),
  then use that 2bit file to forge the BSgenome ....
  package (see the BSgenomeForge vignette in the BSgenome package for
  how to do this).
  - Herve Pages

Here is another similar question with more info about replaceLetterAt: https://support.bioconductor.org/p/26199/

I noticed those. Question: can you think of another algorithmically efficient way to accomplish this? One way would be to convert the GRanges into a BED, then used this to create a new FASTA with the SNPS, and loaded this into BSGenome and then used getSeq(). That's very complicated though

I should mention that this is a section in the official documentation---the creators appears to have "not gotten around" to implementing this. See " Injecting known SNPs in the chromosome sequences" here: https://bioconductor.org/packages/devel/bioc/vignettes/BSgenome/inst/doc/GenomeSearching.pdf


text:
Tools to create annotated table of variants from VCF
The problem: I have a VCF file, a reference genome, and a bunch of annotations for the reference (genes, repeat regions, etc.) as GFF or BED files.
What I would like is a tool that takes all of this as input and outputs a tab- or comma-delimited table containing as much information as possible. Potential columns in the output include:

Type of variant (SNV, structural, etc)
Details of variant (e.g. reference base, variant base, coverage, position, etc)
Annotations overlapping variant
Annotations near variant (e.g. is it just upstream of a gene)
If it appears in a coding region, does it change the amino acid.

There are many tools that do something like this. But to the newbie (like me) it is not clear which tools are worth starting with. Since the majority of tools take a little effort to get working in the first place, my question is:
What tool comes closest to doing what I have described, is known to reliable, and is likely to be maintained for the next few years.

I have used FEATnotator and I think it can provide all of the columns you would like to see. There are many output files generated, but the consolidated output has the following columns:

Chromosome
Position
Column_3
Consensus_Allele
Annotation_Signature
Reference_Base
Alternate_Base
Transition/Transversion SNP_Type
Premature_STOP_Gained   STOP_Lost
START_CODON
STOP_CODON
SPLICE_SITE
InterGenic
Gene_Body
Intron
Exon
Coding
UTR
Transcription_Start_Site
Nearest_gene
Distance

here is some example output records from using a VCF, reference genome and GFF annotation file (sorry about the crappy formatting...there are a lot of fields!):
Chr01   28177   T       G       NA      T     G TRANSVERSION    NA      NO      NO      NO      NO      NO      YES     NO      NO      NO      NO      NO      NO      Eucgr.A00211.v2.0       16054
Chr01   44876   C       G       gene:Eucgr.A00211.v2.0   mRNA:Eucgr.A00211.1.v2.0[intron]       C       G       TRANSVERSION    NA      NO      NO      NO      NO      NO      NO      YES     YES     NO      NO      NO      NO      NA  NA
Chr01   46819   A       G       gene:Eucgr.A00211.v2.0   mRNA:Eucgr.A00211.1.v2.0   CDS:Eucgr.A00211.1.v2.0.CDS.4[AAA - K => GAA - E (MISSENSE)]   exon:Eucgr.A00211.1.v2.0.exon.4      A       G       TRANSITION      MISSENSE        NO  NO       NO      NO      NO      NO      YES   NO        YES     YES     NO      NO      NA      NA

Hail might be an option for you.
It is actively developed by a growing team at the Broad.  It is rigorously tested (continuous integration, continuous deployment, bug reports get regression tests, blah blah blah).
It was designed to solve this problem (among others).  It can import a variety of formats, including VCF, TSV, UCSC BED, JSON and interval files.  (We don't have explicit support for GFF, but we can probably handle them with general facilities.  If not, get in touch and we'll add support.)  It can call out to VEP (and soon Nirvana, Illumina's VEP rewrite).  It has general facilities to transform, filter, clean and query data.
What's more, we've curated a large collection of annotation resources (currently ~22 databases of annotations + VEP), hosted in a public bucket in Google cloud, and built an interactive query builder to select which resources you want to use.  Get in touch if you'd like us to add additional resources.  You can read more about it here.
Disclaimer: I work on Hail.

I've left the accepted answer as FEATnotator, just because this appears to do exactly what I asked in the question. But I wanted to say thanks for letting me (us?) know about Hail, which certainly sounds extremely useful. The active development and openness to feature requests is a huge plus.

snpEff is a great tool for annotating VCF files and you can add custom reference sequences.
http://snpeff.sourceforge.net/
Guide to add custom annotation files in snpEff
https://gatkforums.broadinstitute.org/gatk/discussion/50/adding-genomic-annotations-using-snpeff-and-variantannotator
There are a bunch of pre curated annotation datasets available in their database.