14

It's difficult to get this to go massively quicker I think - as with this question working with large gzipped FASTQ files is mostly IO-bound. We could instead focus on making sure we are getting the right answer. People deride them too often, but this is where a well-written parser is worth it's weight in gold. Heng Li gives us this FASTQ Parser in C. I ...


11

I think you can try dendextend, in this manual there is an example of coloring the branches. I don't think it is exactly like your coloring, but with a little tweaking you might get your colorscheme in there. The manual mentions an argument called color_lines for the function tanglegram(): # The `which` parameter allows us to pick the elements in the list ...


11

Simulators designed specifically for Oxford Nanopore: NanoSim NanoSim-H SiLiCO ReadSim DeepSimulator General long read simulators: Loresim Loresim 2 FASTQsim LongISLND For an exhaustive list of existing read simulators, see page 15 of my thesis, Novel computational techniques for mapping and classifying Next-Generation Sequencing data.


9

By chance, just today I've heard of a nanopore read simulator, NanoSim. It is released under a GPL license. I have never used it, though...


7

In addition to the already mentioned NanoSim, there is also SiLiCO and ReadSim (although it hasn't been updated in over 2 years, so I am not sure how relevant it is at this point considering how fast the technology is progressing).


6

The following is more than twice as fast; however, wc counts newline characters as well. We thus need to subtract the line count from the base count (using Bash): fix_base_count() { local counts=($(cat)) echo "${counts[0]} $((counts[1] - counts[0]))" } gunzip -c "$file" \ | awk 'NR % 4 == 2' \ | wc -cl \ | fix_base_count All the ...


6

I've never tried this myself, so I don't know how easy this is... One option would be to start with GMAP, which is meant to align whole transcripts against the genome. The really nice thing about this is that it can directly produce GFF3 files. You can then use that with your Ensembl GTF with cuffcompare or whatever the equivalent is in stringTie. You ...


6

pigz | awk | wc is the fastest method First off for benchmarks with FASTQ it's best to use a specific real-world example with a known answer. I've chosen this file: ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01815/sequence_read/ERR047740_1.filt.fastq.gz as my test file, the correct answers being: Number of reads: 67051220 Number of bases in ...


6

I would also recommend two very recent Hi-C visualization frameworks (with some public data available in both): HiGlass and JuiceBox.


6

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 '...


6

NOTICE: I have altered my answer slightly from the original as I have turned the original script into a pip installable program (with tests) and have updated the links and code snippets accordingly. The essence of the answer is still exactly the same. This is something I have been meaning to get around to for a while, so thanks for the prompt. I have ...


6

I like bedtools getfasta. My typical option set is bedtools getfasta -fi <reference> -bed <gff_file> -name -s. Be aware of the -s to make sure you are pulling the correct strand. I like bedtools because it is a versatile tool overall for handling bed, gff and vcf file manipulations. # bedtools getfasta Tool: bedtools getfasta (aka ...


6

Found a solution, using D-Genies, worked great. Some examples from their website: Thanks to @user172818.


6

R supports logistic regression, which would seem to be the most efficient method for tackling this question. Assuming the "Chemo" variable is the type of chemo the code would be something like: glm( (response_to_chemo == "yes") ~ BMI + Chemo + Predictor + DJANGO + gender, family="binomial") EDIT: Corrected a typo (added a missing double quote)


5

Have you tried Mauve alignment? Its pretty easy once you become familiar and has a GUI for further ease of use. Additionally there are a few online tutorials on how to re-order contigs/ scaffolds using this software. Heres one I use when in need. Mauve contains a function called Mauve Contig Mover (MCM) which can be used to a) compare an assembly to a ...


5

As per my answer to @_julien_roux on twitter: Trying to find novel transcripts within the context of an existing annotation is much less straightforward. You probably need to do a "genome-guided assembly" with Trinity and PASA: http://pasapipeline.github.io/#A_ComprehensiveTranscriptome We did something similar in much simpler organisms in our recent paper:...


5

I get fairly quick results with my fastx-length.pl script, with the added bonus of being able to handle multi-line FASTQ files and displaying additional read-length QC statistics: time zcat albacored_all.fastq.gz | /bioinf/scripts/fastx-length.pl > /dev/null Total sequences: 301135 Total length: 283.902419 Mb Longest sequence: 5.601 kb Shortest sequence: ...


5

The following bit of python code should work: #!/usr/bin/env python import sys lastTranscript = [None, None, None, []] # ID, chrom, strand, [(start, end, score), ...] def getID(s): """Parse out the ID attribute""" s = s.split(";") for k in s: if k.startswith("ID="): return k[3:] return None def dumpLastTranscript(): ...


5

You can try doing standard differential expression, but I worry that the between-sample normalization will work poorly. Personally, I would do peak calling instead, followed by diffBind. You have a few tools to choose from when it comes to this. In the past, I've rolled my own methods for this using MACS2 and genomic alignments (I then converted those to ...


5

If I understand the question correctly, you'd like to plot the positions of the matches to you motif along with a gene model that shows the positions of introns and exons for the different transcripts. This can be accomplished fairly easily with ggbio: library(EnsDb.Mmusculus.v79) library(ggbio) library(biomaRt) library(stringr) library(dplyr) # plot gene ...


5

The manual suggests that each "chromosome" needs to have its own input file. You have to split the dataset by chromosomes. ShapeIT can just phase one chromosome at a time. If you are working on GWA data, you can proceed as follows: for chr in $(seq 1 22) ; do plink --file example --chr $chr --recode --out example_chr$chr ; done You will obtain ...


5

You can do all of that with khmer. For example, abundance-dist-single.py produces a file with columns: k-mer abundance, k-mer count, cumulative count, and fraction of total distinct k-mers. So for question 1 you would sum column 2. For question 2 you would just get thek-mer countassociated with ak-mer abundance` of 1. That package also provides a python API ...


5

I don't think there is a method that would estimate a genome size using raw long reads. The genome size estimates based on raw reads are done by fitting a model to kmer spectra (for instance Genomescope). The kmer spectra built from long reads are really messy due to the high error rate of long reads. That makes fitting of a model quite difficult. These ...


5

I can only speak of drug design (and even then I am terrible at turning down the jargon). In the case of drug design, this is pretty much plan C. Namely, none of compounds that entered clinical trial at the start of the year work (let's call this plan A) and none of the vaccines that are entering now clinical trial work (let's call this plan B although ...


4

I will take the liberty of giving one possible answers to my own question – but I’m very interested in other answers. One analysis type that such data enables is the analysis of transcript switches with predicted potential consequences. I myself have recently developed such a tool called IsoformSwitchAnalyzeR. IsoformSwitchAnalyzeR enables statistical ...


4

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


4

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(...


4

The gffread utility in Cufflinks package might be interesting for you. To generate a multi-fasta file with nucleotide sequences from your GFF file, then you can try: gffread -w output_transcripts.fasta -g reference_genome.fa input_transcripts.gff


4

You need to be careful of terminology. To me, a RIP-seq experiment involves a pull down, followed by a RNAseq library prep. Thus, the whole transcript is captured, not just the binding site of the protein (as in CLIP-Seq, HITS-CLIP, PAR-CLIP, iCLIP or eCLIP). Thus "peak-callers" whether they be designed for calling protein-DNA or protein-RNA binding sites ...


4

A reverse image search on google leads to the exact code used to produce the second image: [Xcoarse, Ycoarse] = meshgrid([0 1 2 3], [0 1 2 3]); [Xfine, Yfine] = meshgrid(linspace(0,3,3000), linspace(0,3,3000)); DataCoarse = [ 1 2 4 1; ... 6 3 5 2; ... 4 2 1 5; ... 5 4 2 3]; DataBicubicFine = interp2(Xcoarse, Ycoarse, ...


Only top voted, non community-wiki answers of a minimum length are eligible