text:
How to safely and efficiently convert subset of bam to fastq?
Question
How can I extract reads from a bam file (produced by bwa-mem) to fastq given a list of reference sequences to filter out?
Potential difficulties

maintaining FR orientation of pair end reads (in bam all the sequences are reference sequences)
keeping R1 and R2 reads
keeping quality scores in the same encoding as original fastq (default illumina phred scores in my case)
bam can be (ana usually is) sorted by coordinates

Almost there solutions

Sujai's perl solution in blobology does exact opposite - getting reads from list of references (so I could just reverse the list). The disadvantage is that the script outputs an interleaved fq file; requires unique names of mates, otherwise R1/R2 information is lost.
samtools + grep them all from fastq files

create a list of read names that do not map to filtered scaffolds. (cut will extract just read names, uniq will collapse pair end read names if they are the same). Then grep read names from fastq files and remove -- separator between hits
samtools view foo.bam | grep -vf list_of_scaffols_filter \
  | cut -f 1 | uniq > list_of_reads_to_keep

grep -A 3 -f list_of_reads_to_keep foo_R1.fq | grep -v "^--$" > foo_R1_filtered_bash.fq
grep -A 3 -f list_of_reads_to_keep foo_R2.fq | grep -v "^--$" > foo_R1_filtered_bash.fq

filter bam & picard tools

Or I could do just the filtering part and use Picard-tools (Picard.SamToFastq), but as usual I am avoiding java as much as I can. I guess
samtools view foo.bam | grep -vf list_of_scaffols_filter \
  | java -jar picard.jar SamToFastq INPUT=/dev/stdin \
  FASTQ=foo_R1_filtered_bash.fq SECOND_END_FASTQ=foo_R2_filtered_bash.fq

The first solution does not really work for me since I do not want to rename all the reads in bam file and I want to keep the R1/R2 information (since R1 and R2 have different error profiles). Both solutions 2 and 3 I find bit clumsy and I am not sure if they are general, I might get some unexpected behaviour if one of reads is mapped second is not... They both relay on the same filtering step.
I was wondering about some pysam solution. I guess it will be much slower, but at least it will be much clearer and perhaps more general. Something like in Convert Bam File To Fasta File - there is pysam solution for fasta (not fastq), almost there...
Story
I have very fragmented reference genome. Some of scaffolds them are too small to works with, and some of them are contaminants (identified using blobtools).
I want to separate reads that are mapping to different groups to separate contaminants, short scaffolds and scaffold that will be used for downstream analysis. The reason is that if we remap all the reads to filtered reference (0.7 - 0.8 of original genome), the most of them (0.95 - 0.99) will still find a place where they map, therefore there is 0.2 - 0.3 of misplaced reads that will obviously have to affect downstream analysis, like variant calling.
This filtering idea is based logic that if the filtered duplicated genomic region will contain some small differences, they will attract their reads (and if I filter them I will improve variant calling) and if they will be exactly same, they will get reads assigned at random, so there is no harm in doing that.

What do you want to do if only one read in a pair maps to a scaffold/contig that you want excluded? Honestly, I'd just name sort the BAM file and write a bit of python to do the conversion and filtering, but that's just me being lazy.

I am not sure yet, but I guess better will be to keep all pair where at least one of the pair maps to any non-filtered scaffold.
In fact if R1/2 have same names, this is exactly how solution 2 is going to behave. If the names will be different I will filter different number of R1 and R2 reads.

I'm not aware of any pre-made program to do this, so I wrote one for you. This will take a BAM file with any ordering and produce properly ordered gzipped fastq files with the filtering as you requested. Internally, this iterates over all of the entries in the BAM file (ignoring secondary/supplemental entries and those where both mates map to your filter list), store the properly oriented sequence/quality/read name in a buffer, and then dumps that buffer entry to disk once the mate is found. This should be reasonably performant (hey, it's python, to don't expect too much), though if you happen to have indexed BAM files then one could think of ways to make this run faster.
Do check the output, since I've only run one test.

You were 3 minutes faster and your code looks way nicer. I will try it out... Thanks a lot.

The only real benefit to my code is that it avoids jumping around in the BAM file looking for mates and that the BAM file doesn't have to be sorted. If you DO have sorted files, then I think it'd be faster to make a list of contigs that aren't excluded, add * to that, and then iterate over it as in my program. That'd save more time.

Ok, I wrote a bit brute pysam / BioPython parser that uses index of bam to get proper order of read pair for R1 / R2 files and bitwise flag. It should not be too difficult to add more sophisticated filtering rules now.
#!/usr/bin/env python3
# 1. arg - indexed bam file
# 2. arg - list of headers to filter
# 3. arg - name pattern for the output reads

import os
import sys
import pysam
from Bio import SeqIO, Seq, SeqRecord

samfile = pysam.AlignmentFile(sys.argv[1], "rb")
header_set = set(line.strip() for line in open(sys.argv[2]))
base = sys.argv[3]

out_R1 = base + 'R1_filtered.fq'
out_R2 = base + 'R2_filtered.fq'

with open(out_R1, mode='w') as R1, open(out_R2, mode='w') as R2:
    for entry in samfile:
        if entry.is_read1 and not entry.reference_name in header_set:
            # pait pair
            entry_R2 = samfile.mate(entry)
            # do some sequence gymnastics with R1
            seq_R1 = Seq.Seq(entry.seq)
            if entry.is_reverse :
                seq_R1 = seq_R1.reverse_complement()
            # do some sequence gymnastics with R2
            seq_R2 = Seq.Seq(entry_R2.seq)
            if entry_R2.is_reverse :
                seq_R2 = seq_R2.reverse_complement()
            R1.write('@' + entry.qname + '\n' + str(seq_R1) + '\n+\n' + entry.qqual + '\n')
            R2.write('@' + entry_R2.qname + '\n' + str(seq_R2) + '\n+\n' + entry_R2.qqual + '\n')

samfile.close()

There are some ugly bits, please feel free to make it nicer.


text:
Formula for k-mer coverage
Let $C$ be base coverage, $R$ is the length of reads and $K$ is the length of $k$-mer. Then $k$-mer coverage $C_k$ can be computed as $C_k = C\cdot(R - K + 1)/R$.
Could someone please explain why is this equation valid (I'm mostly confused as why it is divided by $R$)?
Source: Velvet manual Section 5.1

$C_k$ is defined as the number of reads containing a k-mer. The fraction of a read available to contain a k-mer is $(R-K+1)/R$, which is the number of k-mers in the read divided by its length. That times the nucleotide coverage ($C$) is then the expected coverage of the k-mer.

I was still puzzled from the answers, so I tried to calculate with all the steps. I take this definition "$C_k$ is the number of reads containing a k-mer." and corresponding definition for coverage ($C$): "$C$ is the number of reads covering a base".
Coverage is $C = \frac{T \cdot R}{L}$, where $T$ is total number of reads, $R$ is read length and $L$ is length of genome. Given the $C_k$ definition, $C_k = \frac{T (R - K + 1)}{L-K+1}$, where $R - K + 1$ is just number of kmers in a read, and $L-K+1$ is number of kmers in a genome. Then,
$$C_k = \frac{T (R - K + 1)}{L-K+1} = \frac{T (R - K + 1)}{L-K+1} \cdot \frac{R}{R} = \frac{R - K + 1}{R} \cdot \frac{T \cdot R}{L - K + 1}$$
since $L >> K$, we can approximate $L - K + 1 \approx L$, then we reduce the expression to
$$\frac{R - K + 1}{R} \cdot \frac{T \cdot R}{L} = \frac{R - K + 1}{R} \cdot C$$
which is the formula for $C_k$.

source: Stack Exchange

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

sort: votes


text:
Do any publicly available databases detail protein structure and functional domains?
I am interested in finding a database that takes a gene or protein name as input (possibly with the option to specify transcript) and gives information about the protein's functional domains in terms of either specific residue ranges, genomic coordinates, etc. Ideally, the database would include 3D conformational information as well as the functional domains. Does such a database exist? It seems that this information can be inferred for a given protein using various publications and references but I have yet to find a single collection containing these data all in one place.

You know about PDB I guess?

@Chris_Rands yeah PDB is great, but to me it seems a bit more like you need to know what you're looking for (ex. EGFR T790M with some kind of binding agent present). I'm looking more for almost a literal diagram or table of the gene detailing relevant functional and confirmation information in a more simplified way. Definitely a source I will be using though.

Pfam might help too.

Some of this information (at least some domains, active sites, etc) is available from UniProt.
If you want to download their whole database, you can search without specifying any terms and then click the Download button.

This seems like exactly what I was looking for, for some reason I haven't come across this in my searching. Seems very simple and to the point, also a plus that you can reference the db locally. Thanks!

EnsEMBL also has this. Search for your gene of interest, choose your transcript, go to the page of its protein product(s), and select "Domains & Features" from the right-hand menu (using human p53 as an example):
Domain source       Start End  Description                                            Accession InterPro
PANTHER             3     331  -                                                      PTHR11447:SF6-
Pfam                5     29   p53 transactivation domain                             PF08563   IPR013872
PANTHER             3     331  p53 tumour suppressor family                           PTHR11447 IPR002117
PRINTS              116   142  p53 tumour suppressor family                           PR00386   IPR002117
PRINTS              158   179  p53 tumour suppressor family                           PR00386   IPR002117
PRINTS              213   234  p53 tumour suppressor family                           PR00386   IPR002117
PRINTS              236   258  p53 tumour suppressor family                           PR00386   IPR002117
Prosite_patterns    237   249  p53 tumour suppressor family                           PS00348   IPR002117
PRINTS              264   286  p53 tumour suppressor family                           PR00386   IPR002117
Pfam                95    288  p53 DNA-binding domain                                 PF00870   IPR011615
Superfamily         97    287  p53-like transcription factor DNA-binding              SSF49417  IPR008967
Gene3D              94    297  p53/RUNT-type transcription factor DNA-binding domain  2.60.40.720IPR012346

In addition to UniProt and EnsEMBL as suggested above you might like:

InterPro, functional analysis of proteins by classifying them into families and predicting domains and important sites.
ECOD, a classification of protein domains by evolutionary relationship.
Pfam, a collection of protein families.
CATH, protein structural domains.
Genome3D, which combines various resources to structurally annotate proteins.

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

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


text:
Fasta Sequence Identifier format?
I have been analyzing some virus DNA from the NCBI databases. The fasta sequences that I recieve from them have header lines that look like this:
>gi|61393989|gb|AY848686.1|

I think that the second number is the GenBank ID, but I have several questions:
1) What does the "gi" in the first ID stand for?
2) Is this a standard format for storing a series of IDs associated with sequences? That is, should I expect to see this format elsewhere and be prepared to parse fasta headers like this in code that I write?
3) It seems that AY848686 is also a valid GenBank ID, and I see the .1 missing/present in GenBank IDs in various places. What does the .1 mean at the end of GenBank IDs and why is it there/does it need to be there?
Thanks in advance for you answers.

The gi is an abbreviation for "Genbank identifier".
This is a pretty standard convention used by data stored in NCBI databases. There used to be a pretty comprehensive description of the conventions used at NCBI (I wouldn't say it was a standard or specification, just convention) here, but this page is no longer available it seems.
The AY848686 label is an accession value, and the appended number is a version number.

See this page for more info.
The short answer is: these are all internal conventions. None of it needs to be there, except perhaps to ensure NCBI's internal tools work as expected.

UPDATE: I found a non-official copy of the old docs here. Look for the section starting with "FASTA Defline Format". I've copied it here to prevent further problems with bitrot in the future.
FASTA Defline Format
--------------------
The syntax of FASTA Deflines used by the NCBI BLAST server depends on
the database from which each sequence was obtained.  The table below lists
the identifiers for the databases from which the sequences were derived.

  Database Name                         Identifier Syntax

  GenBank                               gb|accession|locus
  EMBL Data Library                     emb|accession|locus
  DDBJ, DNA Database of Japan           dbj|accession|locus
  NBRF PIR                              pir||entry
  Protein Research Foundation           prf||name
  SWISS-PROT                            sp|accession|entry name
  Brookhaven Protein Data Bank          pdb|entry|chain
  Patents                               pat|country|number
  GenInfo Backbone Id                   bbs|number
  General database identifier           gnl|database|identifier
  NCBI Reference Sequence               ref|accession|locus
  Local Sequence identifier             lcl|identifier

For example, an identifier might be "gb|M73307|AGMA13GT", where the "gb" tag
indicates that the identifier refers to a GenBank sequence, "M73307" is its
GenBank ACCESSION, and "AGMA13GT" is the GenBank LOCUS.

"gi" identifiers are being assigned by NCBI for all sequences contained
within NCBI's sequence databases.  The 'gi' identifier provides a uniform
and stable naming convention whereby a specific sequence is assigned
its unique gi identifier.  If a nucleotide or protein sequence changes,
however, a new gi identifier is assigned, even if the accession number
of the record remains unchanged. Thus gi identifiers provide a mechanism
for identifying the exact sequence that was used or retrieved in a
given search.

For searches of the nr protein database where the sequences are derived
from conceptual translations of sequences from the nucleotide databases
the following syntax is used:

                     gi|gi_identifier

An example would be:

        gi|451623           (U04987) env [Simian immunodeficiency..."

where '451623' is the gi identifier and the 'U04987' is the accession
number of the nucleotide sequence from which it was derived.

Users are encouraged to use the '-I' option for Blast output which will
produce a header line with the gi identifier concatenated with the database
identifier of the database from which it was derived, for example, from a
nucleotide database:

        gi|176485|gb|M73307|AGMA13GT

And similarly for protein databases:

        gi|129295|sp|P01013|OVAX_CHICK

The gnl ('general') identifier allows databases not on the above list to
be identified with the same syntax.  An example here is the PID identifier:

        gnl|PID|e1632

PID stands for Protein-ID; the "e" (in e1632) indicates that this ID
was issued by EMBL.  As mentioned above, use of the "-I" option
produces the NCBI gi (in addition to the PID) which users can also
retrieve on.

The bar ("|") separates different fields as listed in the above table.
In some cases a field is left empty, even though the original
specification called for including this field.  To make
these identifiers backwards-compatiable for older parsers
the empty field is denoted by an additional bar ("||").

Thanks for the quick and accurate answer. As a follow up question- is there some reason why NCBI has accession numbers as well as an ID numbers?

@JonDeaton history, mostly. The GIs are far older. Also see the "Historical Note" section in the link Daniel gave you.

@JonDeaton For parsing be aware that NCBI have phased out GI numbers from fasta headers since Sept 2016. I couldnt find a new source for this info but please see this note from their blog https://ncbiinsights.ncbi.nlm.nih.gov/2016/07/15/ncbi-is-phasing-out-sequence-gis-heres-what-you-need-to-know/

That is a GenBank gene ID and it (apparently) stands for "GenInfo Identifier". The full GenBank ID is gi|61393989.

Yes, this is the standard format for sequences provided by NCBI. However, you can have all sorts of extra details appended to that. They (NCBI) generally use | as a field separator, so you could try to write code that parse that, but I really wouldn't recommend it. The FASTA headers have traditionally been free. People can put anything they like there apart from a \n, so I really wouldn't put any effort into writing code to parse these headers.
If you need a short name for the sequence, just take the first N characters or, if it matches ^..|.+?|, take the characters until the second |. Seriously though, don't bother parsing this it is very rarely helpful and is almost certain to break as soon as you use a different database.

The AY848686.1 is a GenBank accession. The gb indicates that it's a GenBank accession and the .1 means it is the 1st version of that accession.

As for why there are two (source):

Unlike the gi number system, in which sequence identification numbers were not necessarily consistent across the databases (e.g., GenBank and EMBL could each assign their own gi number to a sequence), the new system is designed to ensure consistency. It is also designed to show a relationship between a sequence identification number and the accession number of the record in which it is found. In contrast, GI numbers are assigned consecutively and bear no resemblance to the accession number. Finally, the new system allows the assignment of alphanumeric protein IDs to proteins translations within nucleotide sequence records. The protein IDs contain three letters followed by five digits, a period, and a version number.


text:
Reject reads with low quality bases from a Bam file through pysam
I have a code below:
def FilterReads(in_file, out_file):

def read_ok(read):
    """
    read_ok - reject reads with a low quality (<5) base call
    read - a PySam AlignedRead object
    returns: True if the read is ok
    """
    if any([ord(c)-33 < _BASE_QUAL_CUTOFF for c in list(read.qual)]):
        return False
    else:
        return True

_BASE_QUAL_CUTOFF = 30

bam_in = pysam.Samfile(in_file, 'rb')
bam_out = pysam.Samfile(out_file, 'wb', template=bam_in)

for read in bam_in.fetch():
    if read_ok(read):
        bam_out.write(read)

This code works fine by

Rejecting reads having a base with phred quality score below five
But it first takes a BAM file
And then creates a filtered BAM file for further analysis but this takes a lot of time.

So is there a way to reject these reads using pileup in pysam so that I may not have to create a file then sort it and again read it?
Or can I modify this code to perform the same function

Why are you sorting the file prior to iterating over it?  Your function doesn't require it to be sorted and since you are not using the random indexing pysam doesn't require it either.

if …: return False else: return True is an anti-pattern. Just write return not …. Or, better yet, logically invert the condition; in your case, you could write return all([ord(c) - 33 >= _BASE_QUAL_CUTOFF for c in list(read.qual)]) (this follows from straightforward application of de Morgan’s laws).

This is probably not the cause of the slowdown, but you can avoid two list creations per read by refactoring line 9 as any(ord(c)-33 < _BASE_QUAL_CUTOFF for c in read.qual):

I don't think pileup is going to save you any time here, as it has to load the same data that you are loading with the fetch method.

For can do this by accessing the basecall qualities from PileupRead.alignment. For example:
_BASE_QUAL_CUTOFF = 30
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            qual = read.alignment.query_qualities
            if all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual]):
                # do some stuff

This involves checking every read at every position. You could memoise this task:
def check_quality(read):
    try:
        return passed_quality[read.query_name]
    except KeyError:
        qual = read.alignment.query_qualities
        passed = all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])
        passed_quality[read.query_name] = passed
        return passed_quality[read.query_name]

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff

Or, using functools.lru_cache (thanks @winni2k):
from functools import lru_cache

@lru_cache(maxsize=1000)
def check_quality(read):
    qual = read.alignment.query_qualities
    return all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff

FYI: There's a decorator in the library that does this and more! https://docs.python.org/3/library/functools.html#functools.lru_cache


text:
creating graph of distribution of blast hits on the query sequence
I'm attempting to visualize the results of my BLAST search in a way similar to the graphical display of the distribution of blast hits from the web BLAST.
For example, from my BLAST search:

except I would like the graph to contain information including information for the distribution of all 6139 hits across my query sequence.
Since this is a large number of sequences to visualize a distribution, I am considering using a score of the number of hits on a certain region, and having an output more similar to:

If I were to run my blast using the BLAST+ command line tool, and read the results into R, which parts of the output should I plot in order to recreate the graphical display of the distribution? Ideally I would like to make a reusable object that can make an equivalent plot for any BLAST input.
Thanks.

How do you get the score of number of hits on regions? I think you should get the score first and plot the score and region in R. If you want aliginment visulization, maybe this link help

What format would you use for the blast output? One of the standard ones, or are you open to defining your own? And would this be plotting the number of hits at each position of your query or their score? Or maybe their overall identity? What exactly is "Conservation" on your example axis?

@terdon I'm not wedded to any specific type of blast output, and am open to defining my own for this problem. I think that I would choose to plot the number of hits at each position instead of score. this link provides a description of the method used to make the conservation plot in my post in detail.

But the number of hits at a given position is useless information if you don't take into account the quality of the relevant HSPs. You seem to want to plot conservation, if so, that's the wrong way to do it. For instance, you could have hundreds of hits at position 12, but if those hundreds are all tiny and very badly conserved, they are irrelevant.

I don't want to plot conservation. I want to know which regions of a given human protein are the most similar to bacterial sequences. If that is my question, does size/conservation necessarily matter?

You probably want to include query start (qstart) and query end (qend) in your blast output.
Something like this:
blastn -query your.fasta -out blast.out.txt -db your.db -outfmt '6 qseqid sseqid qstart qend length evalue'

In R you can take the "qstart:qend" from each line for density plot.
There are many ways in R to plot the densities of these start and end amino acids.
Let me show an example with a small data frame:
qstart <- c(200, 300, 250, 400, 500)
qend <- c(300, 450, 400, 600, 650)

df <- as.data.frame(cbind(qstart,qend))

aa <- vector()
i=1
for(i in 1:5){
aa <- append(aa, c(df[i,1]:df[i,2]))
i+1
}

hist(aa)
dens <- density(aa)
plot(dens)

Thanks @b.nota I am giving this a shot right now and will see how it turns out. This may be exactly what I am looking for.

I'm currently working on making the density plot and having some issues. I have a new column with "qstart:end" for each row in my data frame, but when I try to make a density plot with this column R throws an error. How can I get around this?

Great man! There are many ways in R how to get the data and make a plot from it. I'll give you a small example in the edit of my answer. Good luck!


text:
When performing differential expression analysis, should genes with low read counts be removed before or after normalization?
I have RNA seq data which I've quantified using Kallisto.  I'd like to use tximport to transform the read count data into input for EdgeR, following the R code supplied in the tximport documentation:
cts <- txi$counts
normMat <- txi$length
normMat <- normMat/exp(rowMeans(log(normMat)))
library(edgeR)
o <- log(calcNormFactors(cts/normMat)) + log(colSums(cts/normMat))
y <- DGEList(cts)
y$offset <- t(t(log(normMat)) + o)
# y is now ready for estimate dispersion functions see edgeR User's Guide

I would ideally follow this up with EdgeR's exactTest.  However, I'd also like to remove genes with low read counts using code such as this:
filterCounts <- function(counts) {
  cpms <- cpm(counts)
  keep <- rowSums(cpms > 1) >= 1
  counts[keep,]
}

The only issue is that I can't tell whether I should remove the genes with low read counts before or after I normality the counts while importing with tximport in the first code.  I'm thinking that it would be best to do this before, editing out the abundance, counts, and length rows of the txi data frame which correspond to lowly expressed genes, but I would like this to be confirmed before I proceed.

The more genes you have the more robust the scaling factor is (among the reason why one doesn't normalize to ERCC spike-ins without a compelling reason), so I suppose in theory it's better to filter after determining the scale factor. Having said that, I'd be surprised if the results changed much either way. Unless you end up filtering out a LOT of genes the scaling factors should be fairly robust to the change.

Since edgeR uses a robust estimator for the scaling factor (TMM?), you’re right, it shouldn’t make much difference.

tximport's documentation implies that I need to supply my own offset when moving data into edgeR.  If I create an offset matrix using my above code and then remove a large portion of lowly expressed genes, will the previously created offset matrix be incorrect?  Do I even need to supply an offset in this case as I think edgeR can calculate it automatically.


text:
What database do gene IDs starting with sp_v2 refer to?
I am a student trying to analyze GEO2r datas for one of my courses.
The IDs given in the output are different for different series. I need to convert all of them to a similar format.
In this process I encountered the following type of ID which I don't know its origin:
”ID"    "logFC"
"SP_v2 4634"    "-0.9897758"
"SP_v2 3382"    "-0.8391782"
"SP_v2 4210"    "-1.1693583"
"SP_v2 2117"    "-1.0504727"
"SP_v2 3488"    "-0.9756444"
"SP_v2 1128"    "-0.8289103"
"SP_v2 2735"    "-0.8629999"
...

Each one of the rows represent a single gene.
My question is that what is this ID? The GEO accession is GSE97750.

Looks like some home made accession code, from which GEO set is it?

The GEO accession is GSE97750

Yeah it is a custom array. Annotation can be found in the platform (GPL22166) data table. See answer by sjockell.

This is possibly a wild goose chase, but a lot of searching led me to this sample:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM2536852
Which appears to exhibit accessions of the right format. This in turn leads to the platform:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL22166
Which is a custom spotted cDNA array, which the metadata seems to suggest is a human platform.
The full annotation table gives gene symbols and Entrez IDs for most of the probes on the array:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc=GPL22166&id=53508&db=GeoDb_blob144
To take the example of the first line of your results:
ID          ORF     Entrez gene SEQUENCE
SP_v2 4634  ETNK1   55500       AAAGCAGCTTCATCTTTCAAAATTGATTTGCTCTGGTTTT

The Entrez Gene record for that ID matches up:
https://www.ncbi.nlm.nih.gov/gene/?term=55500


text:
Documentation for Transfrag class codes (cuffcompare)
Is there any extensive documentation or description for each class of Transfrag class codes as reported by the cuffcompare tool in the Cufflinks package?
Official doc might not be the best.
E.g. what does contained mean (class code: c)? Or, what is a generic exonic overlap (class code: o)?

You can get a bit more information by looking at the source code comments. Contained is for single exon transcripts that are primarily contained (>80% overlap) within an exon of an already annotated transcript. Generic exonic overlap is for cases where a transfrag exon overlaps that of an annotated transcript. The strand doesn't seem to matter in such cases.


text:
Extract nanopore read ID & start times from fastq file
I have a fastq file from minION (albacore) that contains information on the read ID and the start time of the read. I want to extract these two bits of information into a single csv file.
I've been trying to figure out a grep/awk/sed solution, but without success.
E.g.
@93a12f52-95e5-40c7-8c3e-70bf94ed0720 runid=17838b1d08f30a031bf60afabb146a8b0fba7486 read=12217 ch=492 start_time=2017-07-04T06:42:43Z
CTATTGTCCCCTGCCGCTGCCCCTCCTGCTACGCCCCACTGCCTCACCAGCCGTTACGGTCGCCCCCCATCGCATGCCTTTACACACACACTTCTTTACACATGCTATCTTCCC
+
""*)&$-.,-(#"&'$%''+16#"$##&)%%/"+*(*(&#"&'%1"+)#)""%$$#&&'1%"'8>MJ<#'&%'.2'.$(&#'()'&&%'$('%"%%%..$#"&"#+&,*$%"#"
@ff37e422-a25f-404c-8314-ef1733f9c30c runid=17838b1d08f30a031bf60afabb146a8b0fba7486 read=8432 ch=200 start_time=2017-07-04T06:56:41Z
CGATGGCCGTATGCTTTTGTTATGAAGCGAAAAGCTGCTCGCTTCTCTAGATAATAATGATGTGGCGAAAACGCTATGCGATTCGTTGACATACATGATGGCGGATTTATCTACCACTTTGTGGCATGCTTTTCTCGCCAGATAATGGAATGTTTCTCTGGCGGTAATGGATAGTATCAAATCTCACTAGCCCATTCTATAAGCGCATCCGCATGCACTAGTTCTTGATTCGATCGTCCTCTAGCATGTTCGAAGAAGATAGCATTCACTATCATCATCGCCTCAGGTAAGTTTATTCGGTTGGGGCGTGTGAAGGCAAACACCTTGTTGTCCAGTAAGTTTTCAGTTACTATAACTTAAAGTCGCACATGAATCTAGTCTCCTATTCCCCACCCATGATCCCACTCACACATTTCTACAGAGATGTGGTTAGAAATTTTCATATTAGGTCAGCTTTGACTCAATAAGACATAATTCTTCACTGAATGACTTTTTAAGAACCACCAGGACCAGAGAGAACCAAGAGAGTGGTACCTCTTAAAACACAATAAAGTGATTCAGCCTTAGCCATTGGATTCTGGAGGACCTTGAACCATGTGGGAAGCAGCTCAGGGTGGCCATGTACTATACTGGCGGGTAAGCTTCTGGAGTGCTAGGTTCTTTTTGTCTTTTCTTAAGCATTGCCGCCAGTTGATTGGGTTTTGAACATAAAATAATGCGCCACCAGCAATTCCAGATTTGTTCCTACGGGATAGATTTGTTCAGTTCTAGCATTATGCTTCACTAACCAGATGCGGGCCCTAAGTCCTTCACTTGGAATATTGGATTGGATCATGAGAATATTCTGTCTGAAGCTCGTCATTAATTTTGTTACAAAATAGAGCTTTTTGACTGGAAGTACCACCATACGTGTTCTCAAACTTCAGCATTTTTAGAACTTCCCACGGCATCTTGACCCTTTTCACAGCATGGATAGTCAGGCAGCAGTGAACTTTGTGACTCTTTAATGCCTTCACTTTTCTCTCAGTTTCCCCGCCTTGCGTTATCTTTACTCGTCTTGGGACTTTTATCCCAATGCCAGCCTTCTACCCTGAGACCTCAGTGGGTCATCATCCCAGCCCGGGACATCTCATCCCATCATTTATGGGCTGTTGTGTTTTTTTCAAAACCTAGCCCTCTCAGGAGGAGGAGGAGTGGGAGTCAGTTCAGTGAGGAGGATTAGGATGATCTGAAATGTAAGCACATATAAGCGAAGCACTTATTTTGGGTTGGGTCCTCACGGTGGACATAAGATCGCCTTATGTGTTTAGTAAGCCATTCCTAGCTCTCAATGGCGTGATTACATAGAAGCGTGAGGGATCAGTCCTATGGAAGACTAGGAAGTAAATGAACAAAATATATTAACCATAGAAGTCTCATGGGTCGCTGTAGCCAAAAGATTAACACTTTTGACTACATTGTGGTTTTAGGCATTGAAACAAAAACTTTGAGTCTCCTAAACAAATGAATGGAAAATAGTAGCGAACTTCGATTCCTAACATTAAATCTAGAAATAGCAAGTTAGTTTAAAGACTTTATTTAGCTTTGCTTGCTATAATGAAAACCTTGCCTCCCGGTCGGGGCCATTGTGCCTGAAGCTAGCTTATTGTCTCCTCGAGCTCCCAGCTTCAGCAACTCCTTTTGAAGCGTTTGTCTCAGCTTGGATCTTCAGCAGCTCTTGGTGGCTCTTTTGAGCTAGCTCCTCTGAGATCTTGTATTTGGTAGGTCGCTTAGTCATAGTACTTTTCTTTTAACACCCTTCAGCTCTACGATTACATTTGGTTTTGTGGATATCATAATGATTGATGTGAAGATACATTGTACATGTG
+
"#""#"##"'"&%&?JP7+80)'&*+&&(,.>3&*(#%(')*&'*#(-$)"&,63;?844&&#$'3+.>@9;-/259...&:"$"/*"#*#%(.&%)+,/76/+'4B93:;70.1)'1:)*#&()(+'03589/++&)'$-:;@,B5(9+.JAU;7-+&,1212.+329NPSDHRHH6),&)%)+).-*$,<1(.+&$$2&,-0+.&'//"4(%"+#,/)<56%*%*$<4*/;DI+6>((-%*).+--;2-'04-062<@-<>6CP)#.0/+,#(*2+#"$$-2'5(%%&'5+)$$%8+0;>.@9*(,((,*-2/49G9/1/3./.0CQ93E+/-D-'67J*)($-2+'*/&*)+%&'092<C9/5()$&%(.,+-'**)'2/.(//6<59-5-)"**40:A4+834)-(*+;K0-,&)+%)(#-&+&&))+$>9.6<?A;3'/85$/--..)/3DP99HJ9."2;I?@C/5%.4<0+0+&$-05C46208AFQE80"%356+02*9I73//(-*/9706192*/)%)(&'($1749C*KKA*/#"/*)5($'#+7@=$%)/&+%##*&&;.),#(4.,2/'-G03%""#,(2)./,%%1;7#+%$'(03'5)+&.%1)&(%'#*,$3#''-3,),$"),'+,,*$%'%$33ILADPPY@JW\O2-8'9[Y@--*.,""$%$"%#"65',0/4(6.*43B1.#&)+),(6BI729=9&#$.)#,(,(/-286D&),'C;P=)'-'#$)*91*$%(.3LH31,($/(+)'$$*#+'((('*)00:7C1E24(%+$&-,$'++)<)'**"&''"(*(.14++$5&*,%?.F1<10"&%$,)&&*(&$%$/.""+3"#"#"*%5*->87FQ>AKY<9855,&(*@>HAZRQ)'%$"--&/#%.+,/A9F8-1'+"6--*34P-9;=997;7<SUU8((-+.:784/%).&)#*'(*/%7H(1JI-+&+)&9.,/N)"%"&))22/&,)%)'+*/17:8$''(*.*..(/.(%+*&22(+)7GM99A@9)0*1%/470((&((%(937'15084#,($"&'%"'4.,+5<7($*0+*291$*&$%%('+/)'.1-(-.7=K/'%"'*)03<:/<=1(>.((-9-'+>&&&%';,)69V<'1:>FNXA+01:7*#"$'$.%&4BQT7/&03++)'&'*3E.+0(*6..*<.)9*)))'%#(*.&&'"*-IU=&,-))-.,->'G--6..90067%(,0/*%9)&#'""#$,5H/$/#-/&2(.0.+;>9--7.3@+2(/4.*0&,-&'/*%0-**-#,/2&-$*)'+#"*)-*##&$3):C:D6H+',./2?<0>)%#&&3::9;8$#%#('$#$',,,*-346@=.%&&(&%('*+&*&1)'&-..#'+2*)/6499=-2*,0,82,$)($:2%*##%$,.5=-5/((#)2./('($@=4,)%1+176:AHQT,(,)0(()@2'+(18G?($)(*'.)&&+*12SD7,18I2)()*+C:01ETZIPDUOB+.55K7B:A5)*/('')/-1=7E<IGR=43@A68;K/:8.'--)BGEE5-0/5+,2''DXBLW7)#DHTC5AL82+--,&#);+).,'3)/.336/&'&(,+8E346-/%1-)1*$+),@-&-%'$+)5.&&%&#(,2&03&/((,(-5%%)?7I4*').+/7+&***845(*/.)(.0))<?YK3=GHH2&*)%$%#%KLY<9&&+,8;/2*./;+))-)/03$&*$+0'36<J+5/+&+(HGFFM57,59)-,0)/:**"%$&#$'&$%#)*%33VX=@FOL00(''4IGIYOG0=1,,#(*$)-0CFM?8+31E8)(+)(%'"'/""%+3/,1;+7-*''&-%98.$$/'&)&+&(&&"

Should produce

93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

awk '{if(NR%4==1) print $1, $5}' file.fastq | sed -e "s/ start_time=/, /" -e "s/^@//"

The awk command gets the first of every 4 lines, printing the first and fifth "word". sed is then used to strip the initial @ and replace start_time= with ,. The output on your example file is:
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

You could simplify your awk to awk 'NR%4==1{print $1, $5}', if you like.

Since the string start_time will only appear on the header line, or else you don't have a valid fastq file, you can simply do:
$ perl -ne '/^@(\S+).*start_time=(.*)/ && print "$1, $2\n"' file.fastq
93a12f52-95e5-40c7-8c3e-70bf94ed0720,2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c,2017-07-04T06:56:41Z

Alternatively, since you mentioned awk and sed:
$ awk -v OFS=", " '/start_time/{print $1,$NF}' file.fastq | sed 's/start_time=//'
@93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
@ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

Or, doing the whole thing in awk:
$ awk 'sub(/start_time=/,""){print $1", "$NF}' file.fastq
@93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
@ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

And if the @ annoy you:
$ awk 'sub(/^@/,"") && sub(/ .*start_time=/,", ")' file.fastq
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

And in sed:
$ sed -n 's/^@\([^ ]*\).*start_time=\(.*\)/\1, \2/p' file.fastq
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

Or, if your sed supports it:
$ sed -En 's/^@(\S+).*start_time=(.*)/\1, \2/p' file.fastq
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

Finally, since grep can't do replacements, to use it you would have to do something like:
$ grep -oP '^@\K\S+|start_time=\K.*' file.fastq | paste - -
93a12f52-95e5-40c7-8c3e-70bf94ed0720    2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c    2017-07-04T06:56:41Z

And to get the commas:
$ grep -oP '^@\K\S+|start_time=\K.*' file.fastq | paste - - | sed 's/\t/, /'
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

Albacore produces a sequencing_summary.txt file (actually TSV, not CSV) in the same directory as the workspace folder that might have the data that you want in it. Here are the fields present in that file:

filename
read_id
run_id
channel
start_time
duration
num_events
template_start
num_events_template
template_duration
num_called_template
sequence_length_template
mean_qscore_template
strand_score_template
complement_start
num_events_complement
complement_duration
num_called_complement
sequence_length_complement
mean_qscore_complement
strand_score_complement
sequence_length_2d
mean_qscore_2d

We'd thought the same. The issue is that the start_time that albacore calculates is actually the time since the start of a run. This is fine if you keep your minION switched on for a whole run. But if you reset it occasionally (which we do, because it really seems to increase yields) then the start_time column is a lot less useful.


text:
How GFF3 attributes (9th column) varies from one gene prediction algorithm to another
GFF3 files are in tabular format with 9 fields per line, separated by tabs. The first 8 fields share almost same data structure, but the 9th field varies a lot depending on feature type and gene prediction algorithm.
Presently I am trying to build my own GFF module that will parse a GFF file and return one structured data object (a python dictionary structure). For that, I need a good understanding of the 9th field of the GFF format.
So can any one please help me by providing various GFF files?

You might want to look into this python GFF3 parser, which follows biopython's style: https://github.com/chapmanb/bcbb/tree/master/gff

Have you looked into gffutils? It will probably do most of the heavy lifting for you, i.e. parsing.

The first place to start is the GFF3 specification. This is the official word on what is and is not allowed in a GFF3 file. For example, users can define arbitrary attribute keys, so long as they do not begin with an uppercase letter (these are reserved for "official" use).
But your question doesn't seem to be about what is allowed, but what is commonly used. I have a question for you: which gene predictors are YOU using? Or are you using gene annotations produced by others?
Here are a few examples.
NCBI RefSeq
If you're using GFF3 from NCBI, chances are it looks something like this.
NC_007070.3     Gnomon  gene    329235  331223  .       +       .       ID=gene14;Dbxref=BEEBASE:GB42168,GeneID:551678;Name=LOC551678;gbkey=Gene;gene=LOC551678;gene_biotype=protein_coding
NC_007070.3     Gnomon  mRNA    329235  331223  .       +       .       ID=rna27;Parent=gene14;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;Name=XM_624067.4;gbkey=mRNA;gene=LOC551678;model_evidence=Supporting evidence includes similarity to: 66 ESTs%2C 24 Proteins%2C and 99%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 113 samples with support for all annotated introns;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    329235  329459  .       +       .       ID=id117;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    329850  330082  .       +       .       ID=id118;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    330166  330301  .       +       .       ID=id119;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    330376  331223  .       +       .       ID=id120;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  CDS     329333  329459  .       +       0       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     329850  330082  .       +       2       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     330166  330301  .       +       0       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     330376  330416  .       +       2       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1

MAKER
The MAKER annotation workflow (paper, software) is a pretty commonly used gene annotation tool, and produces GFF3 output like this.
scaffold_12     maker   gene    652527  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959;Name=maker-scaffold_12-augustus-gene-0.959
scaffold_12     maker   mRNA    652527  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1;Parent=maker-scaffold_12-augustus-gene-0.959;Name=maker-scaffold_12-augustus-gene-0.959-mRNA-1;_AED=0.24;_eAED=0.18;_QI=0|0|0|0.66|0.5|0.33|3|0|218
scaffold_12     maker   exon    652527  652817  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1203;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   exon    654877  655170  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1204;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   exon    655272  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1205;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     652527  652817  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     654877  655170  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     655272  655343  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   gene    941547  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851;Name=snap-scaffold_12-processed-gene-0.851
scaffold_12     maker   mRNA    941547  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1;Parent=snap-scaffold_12-processed-gene-0.851;Name=snap-scaffold_12-processed-gene-0.851-mRNA-1;_AED=0.22;_eAED=0.22;_QI=0|0|0.25|0.25|1|1|4|661|95
scaffold_12     maker   exon    941547  941631  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1206;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    942343  942367  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1207;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    942780  942920  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1208;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    943200  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1209;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     941547  941631  .       +       0       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     942343  942367  .       +       2       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     942780  942920  .       +       1       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     943200  943236  .       +       1       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   three_prime_UTR 943237  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:three_prime_utr;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1

GeneMark, SNAP, Augustus
If we look at the output of several common ab initio gene prediction tools, none of them actually produces GFF3 by default. Here is some output from GeneMark (paper, software).
543   1   +  Initial         799051     799097         47          1 2 - -
543   2   +  Terminal        799236     799266         31          3 3 - -

544   3   -  Terminal        802357     802514        158          3 2 - -
544   2   -  Internal        802607     802685         79          1 1 - -
544   1   -  Initial         802829     802843         15          3 1 - -

Here is some output from SNAP (paper, software).
Einit   8230    8239    +       9.329   0       1       0       scaffold_12-snap.4
Exon    8848    8869    +       2.484   2       2       2       scaffold_12-snap.4
Exon    10121   10208   +       15.302  1       0       2       scaffold_12-snap.4
Exon    11361   11420   +       5.969   0       0       2       scaffold_12-snap.4
Exon    11471   11535   +       0.921   0       2       1       scaffold_12-snap.4
Eterm   12169   12187   +       18.163  1       0       1       scaffold_12-snap.4
Einit   14569   14668   +       -6.918  0       1       0       scaffold_12-snap.5
Exon    15029   15203   +       -1.023  2       2       0       scaffold_12-snap.5
Exon    16171   16307   +       -9.230  1       1       1       scaffold_12-snap.5
Eterm   16667   16698   +       9.829   2       0       0       scaffold_12-snap.5
Einit   17809   17898   -       -5.390  0       0       0       scaffold_12-snap.6
Exon    17299   17350   -       6.978   0       1       2       scaffold_12-snap.6
Eterm   17029   17129   -       -2.857  2       0       0       scaffold_12-snap.6

Both of these formats appear to be tab-delimited, but that's where the similarity with GFF3 ends.
The output of Augustus (paper, software) looks like GTF (a variant of the GFF format) once you lose all the lines beginning with a # symbol, but Augustus also has a GFF3 output mode.
# Predicted genes for sequence number 1 on both strands
# start gene g1
scaffold_12     AUGUSTUS        gene    2841    3400    1       +       .       g1
scaffold_12     AUGUSTUS        transcript      2841    3400    1       +       .       g1.t1
scaffold_12     AUGUSTUS        start_codon     2841    2843    .       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        intron  3027    3097    1       +       .       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        CDS     2841    3026    1       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        CDS     3098    3400    1       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        stop_codon      3398    3400    .       +       0       transcript_id "g1.t1"; gene_id "g1";
# protein sequence = [MAIKNAEHDLRVIVDAIEGLGLKVAPHKTEAMAFPASALCGRRGAAPPKIRLGGSSILVGSRSRWYISGHSENSSKSP
# RTEGKETTPLQQRDPLDAPLWVSGVVAHCCGGPEGQEGCPGLAAQGSDQGVLRIRDGLLCGYDGCGDHRPRPSDSSAGGGLCRP]
# end gene g1
###
# start gene g2
scaffold_12     AUGUSTUS        gene    4712    15229   0.21    +       .       g2
scaffold_12     AUGUSTUS        transcript      4712    15229   0.21    +       .       g2.t1
scaffold_12     AUGUSTUS        start_codon     4712    4714    .       +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        intron  4858    5591    0.74    +       .       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        intron  5686    15028   0.23    +       .       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     4712    4857    0.72    +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     5592    5685    0.42    +       1       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     15029   15229   0.43    +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        stop_codon      15227   15229   .       +       0       transcript_id "g2.t1"; gene_id "g2";
# protein sequence = [MGRNSHRSCCVVNCKITSAKSDCKFYKFPTAKWKINQRKMWVAAVKRQKYIKDEISHAETQTEITEVTGATKVNYANK
# KYICLLFVRTYVLRMLVDVALSNLRFSLFGIRKSLEIFGQSEKADQTRWRLPSCEMEWIESRKGKMRE]
# end gene g2
###

Parsing attributes
So...how should YOU handle attributes in GFF3's 9th column? That depends a lot on what you want to do with the data.
The most important attributes are ID and Parent, which are used to define relationships between features and subfeatures. (These relationships implicitly define a directed acyclic graph of features, although most GFF3 parsers don't directly support traversal of this graph.)
But not only do you have to handle the pre-defined attributes discussed in the GFF3 specification, you also have to be able to handle any number of arbitrary attributes whose keys you may not know beforehand.
The simplest way to handle this would be to parse the attribute column into a dictionary of key/value pairs. Once it's in this form, it's trivial to see what attributes are there and how to access them.
For example, if we have a feature with the following attributes...
ID=mRNA42;Parent=gene19;integrity=0.95;foo=bar

...we would want it in a dictionary like so.
attributes = {
    'ID': 'mrna42',
    'Parent': 'gene19',
    'integrity': '0.95',
    'foo': 'bar',
}

Python code to parse that might look something like this.
attributes = dict()
for keyvaluepair in attributestring.split(';'):
    for key, value in keyvaluepair.split('='):
        attributes[key] = value

One thing to consider is that it's possible for an attribute to have multiple values (separated by commas). It's not commonly used, but it is valid and something that should be handled or at least checked for.

Thank's @Daniel for your excellent answer.My ultimate goal is to convert any gff3 format to a python dictionary structure. can you please help me convert one gene record(One gene record may have mRNA, exon, CDS etc) into a dictionary structure?

@Arijit Parsing a single line into a dictionary is trivial, but resolving parent/child relationships involves ID cross-referencing and is more complicated. Most GFF3 parsers leave this as an exercise to the coder. I'm working on a Python library called tag right now. Feel free to check it out and see if it meets your needs. http://tag.readthedocs.io/

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


text:
GRCh38 centromeres mappability
Is the mappability of the centromeres in the GRCh38 genome reference similar to each other?
As far as I can remember when GRCh38 came out, the sequence of the centromeres was determined by a combination of sequencing data and software prediction.
Given the way the sequence of the centromeres has been determined, should we expect Illumina 2x150bp (or shorter, 2x75bp) reads to map relatively equally to all centromere sequences?

should we expect Illumina 2x150bp (or shorter, 2x75bp) reads to map relatively equally to all centromere sequences?

No. It has long been established that different chromosomes are associated with different centromeric sequences. It is sometimes possible to tell which chr a read is originated from based on its sequence.
The GRCh38 centromeres are trickier. As I remember, the centromeric sequences were computationally generated with a Markov chain (or something alike) modeled after Venter's genome. GRC can distinguish most chromosomes, but not all. Some alpha arrays are placed onto 2 or 4 chromosomes. The original GRCh38 keeps all 4 copies. When you download GRCh38 for mapping purposes, only one copy is retained; additional copies are hard masked along with PARs on chrY.
If you want to know more, see this paper and this.

I doubt it, unless you are asking if mapping would be similarly bad for all centromeres. Here are some repetitive structures (probably not centromeric) that I've found in the nanopore reads for "human" sample NA12878, produced by the nanopore-WGS consortium:

These structures are consistent in that they repeat lots of times, but the internal patterns can be quite different. Here are a few more:

Given that centromeres need to be uniquely linkable to a single chromosome, it would make sense to me if the centromere internal structure is unique to each chromosome.
It's possible to have highly-repetitive structures which don't map to each other. While I haven't dug too deeply into the human reads, I've looked at the 5 most-compressible regions in an assembled rodent parasite (Nippostrongylus brasiliensis) genome, and found no internal similarities between them:

One of the issues with assembly from Illumina reads is that these highly-repetitive regions get collapsed into a single repeat (or at best a region up to twice the fragment length). With internal repeat units that have identity of over 98%, assembling the true sequence is very difficult, even with an exact knowledge of paired read separation. Even if this were possible, it can be impossible to correctly place a read because multiple internal units could be identical (or similarly different) to the sequenced read.


text:
Difference between CPM and TPM and which one for downstream analysis?
What the difference between TPM and CPM when dealing with RNA seq data?
What metrics would you use if you have to perform some down stream analysis other than Differential expression for eg.
Clustering analysis using Hclust function and then plotting heat map to find differences in terms of expression levels, correlation and pca
Is it wrong to use TPM for such analysis, if yes then when does one use TPM versus CPM.

You can find the various equations in this oft-cited blog post from Harold Pimentel. CPM is basically depth-normalized counts, whereas TPM is length-normalized (and then normalized by the length-normalized values of the other genes).
If one has to choose between those two choices one typically chooses TPM for most things, since generally the length normalization is handy. Realistically, you probably want log(TPM) since otherwise noise in your most highly expressed genes dominates over small expression signals.

in a given case if one would trim adapters from paired end sequenced rna seq data, this would result in different read length, as you say above TPM is length normalized does it mean that this difference in read length is taken into consideration?

@novicebioinforesearcher No, transcript/gene length, which correlates (somewhat) with counts and therefore will tend to drive clustering unless handled in a reasonable way.

Rob Patro also wrote a pretty good article about the topic as well: http://robpatro.com/blog/?p=235

[link] (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8220791/)

Neither CPM nor TPM are well suited here, because neither performs robust cross-sample normalisation (see the blog post Devon linked to).
DESeq2 provides two robust log-space normalisation methods for downstream analysis, the regularised log (rlog), and the variance stabilising transformation (vst). The DESeq2 vignette explains how to use these for things like hclust.

On a more general note, CPM does not account for transcript length differences, while TPM does. If the choice is between TPM and CPM I would therefore use TPM. However, if you are only comparing the same transcripts across experiments, the transcript length is actually invariant so it doesn’t matter (but CPM is still not a good cross-experiment normalisation).

I am confused now , so what is the use of TPM, why does one produce it and when or where do you use it? , in other words which tools, analysis in RNA seq would you use TPM if everything revolves around using counts and pushing it through DESeq2

Tools produce TPMs because they don’t have the information (= the other samples) necessary to perform cross-sample normalisation. Lacking that, TPM is the best they can do. TPM is also useful for within-sample comparisons: It can give you an accurate estimate of how much genes are expressed in a given sample relative to each other.

Plz Correct me if am wrong here, given an experimental design different cell types from a normal mouse say 4 cell types (3 replicates each), sequenced using same library prep but may be at different times. Aim would be to check for set of cell type specific transcripts you would use TPM where as if you want to add a significance parameter (need a pvalue) one would use raw counts based analysis?  I guess the confusion for me is when we use the word  "expression" . People use TPM and call it expression, also use raw counts and call it expression

Both are estimates for expression, given the data. Your use-case sounds reasonable although I would generally prefer to determine “cell type specific transcripts” by comparing different cell types, rather than based solely on a single sample. Which would imply performing differential expression analysis.

Which would imply performing differential expression analysis on?

@novicebioinforesearcher On whatever data sets you wish to compare. It rarely (if ever!) makes sense to describe a gene as being cell type specific without saying “as opposed to these other cell types”. For example, a cell type specific gene may nevertheless be lowly expressed: as long as it’s entirely absent in other cells, it’s cell type specific. This is in fact often the case. You therefore cannot characterise many cell type specific genes without comparing different cell types.

Is VST/RLOG normalised for transcript length? If not is it possible to get length normalised VST/RLOG?

@rmf No, they do not normalise for transcript length; both functions merely change the distribution shape of the counts into something closer to linear. For those purposes where you’d use rlog/vst, accounting for transcript length, is normally not important. You could however apply further transformation (let’s call it rlog-TMP), if you have an application where you need cross-sample as well as within-sample normalised values.

Would you say it's ok to do something like (vst/length)*(10^6) and use that for heatmaps where I actually want to compare expression of one gene to another gene.

@rmf Yes but in a heatmap you usually scale by (gene) row anyway (either explicitly or the plotting function performs the division internally) so the normalisation per transcript length will be strictly a no-op.


text:
Why are some of the gene peptides returned by biomaRt missing an asterisk in the end?
Note: this question has also been asked on BioStars
I am getting peptides using biomaRt library in R for the gene 'BRCA1'. 27 different aminoacid sequences are returned, and 12 of those have an asterisk (*) in the end, while 15 do not. What does that mean? I know that asterisk is used to indicate stop codon. So, does that mean the peptides without an asterisk do not have a stop codon? Is that possible at all? Are all of those valid peptides or not?
Below is the R code I used:
library(biomaRt);
mart = useMart('ensembl', dataset='hsapiens_gene_ensembl');
seq = getSequence(id='BRCA1', type='hgnc_symbol', seqType='peptide', mart = mart);

MDLSALRVEEVQNVINAMQKILECPICLELIKEPVSTKCDHIFCKVLLCCPSWSTVVRS* is an example sequence with an asterisk in the end. MDLSALRVEEVQNVINAMQKILECPICLELIKEPVSTKCDHIFCKSLQESTRFSQLVEELLKIICAFQLDTGLEYANSYNFAKKENNSPEHLKDEVSI is an example without an asterisk.

If it's 3' incomplete that means the evidence used to create it was a fragment. Here's the evidence used to construct BRCA1-214 ENST00000477152.5, a 3' incomplete.
You can see that there's a full length cDNA from EMBL, AK307553.1, which was used to create this model.
The sequence was mapped against the genome sequence to create the transcript model. When searching for the translation, it was found that the open reading frame was started, but there was no stop codon. This suggests that the cDNA AK307553.1 is actually a fragment, that the mRNA was broken or cleaved before it was reverse-transcribed sequenced. We display it in Ensembl in the hope that this will lead someone to identify the full length transcript that it represents.
Many people choose not to work with these incomplete transcripts. You can do this by filtering by Gencode Basic in biomaRt.

If you look at the original data at Ensembl you'll notice that most of these are labeled "CDS 3' incomplete" and have a TSL (transcript support level) of 1, which is as low as it goes. It seems likely that that is simply an incomplete annotation. I'm not surprised that there are a bunch of these for BRCA1, since there was a long time when there were a LOT of legal issues surrounding working on it.

@DevonRyan TSL1 is the most reliable, TSL5 is the least. Also, OP, there are some times of day when you might get an answer instantly, but mostly I would assume a full working day is a reasonable amount of time to wait before impatiently crossposting everywhere. You don't know where in the world the one or two experts who can answer your question are, and what time they'll get onto the web.

@EmilyPerry Gah! Thanks for letting me know about my TSL mistake. How the heck did an annotated "incomplete" transcript get a TSL of 1?

Still won't let me comment. Devon, TSL and incomplete are independent. It's TSL1 because there is at least one non-suspect piece of evidence that covers the whole thing. If all we've got is evidence of a fragment, then that's all we've got, but the evidence of that fragment is strong.

@EmilyPerry Sorry, that's an annoying stackexchange thing, you need 50 rep to comment. Anyway, thanks for clarifying that TSL and "3' incomplete" are independent annotations!

@EmilyPerry I am not a bioinformatician or a biologist, so I don't really understand what a "non-suspect piece of evidence that covers the whole thing" means. Does it mean that this sequence with no asterisk in the end was found as a protein synthesized by the gene? Probably not since there is no stop codon, right? Or does it mean that an mRNA that matches this peptide sequence was measured? The mRNA without a stop codon would be degraded, though.. Sorry if this is a dumb question, but I'd really appreciate some clarification on that.


text:
How to filter intervals (reads or genomic coordinates) that have the exact same 5' or 3' ends?
Say I have reads that overlap some genes that produce small RNAs, but I want only those reads that start at exactly the TSS of the loci. In other words, reads whose 5' end match the 5' end of a genomic feature.

   5'....3'
...+++++++... Gene
...0000000... R1
...00000000.. R2
..000000000.. R3
..00000000... R4

The output should be R1 and R2.
I believe this is a fairly common bioinformatics operation but somehow it doesn't seem to be an option in the general use tools I looked at (py-bedtools, samtools, bedops).
A solution I thought of would be:

reduce the gene coordinates to its 5' end, and
using bedtools bedClosest to annotate the distance between reads and genes, and finally
select reads that overlap (d=0) and write them to a (bam) file.

This involves quite a bit of wrangling data with pybedtools or awk/bash, and I wonder is there is a more elegant solution to this. Solutions in R/Bioconductor are also welcome.

This method won't work as overlapping segments will also have a distance of 0.  You will have to clip the read coordinates as well as the gene coordinates.
From the bedtools documentation:
"-d In addition to the closest feature in B, report its distance to A as an extra column.

The reported distance for overlapping features will be 0."

Presuming you have deepTools installed:
#/!usr/bin/env python
from deeptoolsintervals.parse import GTF
import pysam

anno = GTF(["some/gtf/or/bed/file", "or/more/than/one/if you like/"])
bam = pysam.AlignmentFile("something.bam")
for read in bam:
    s = read.reference_start
    if read.is_reverse:
        s = read.reference_end
    e = s + 1
    overlaps = anno.findOverlaps(read.reference_name s, e, matchType=4)
    # For a GTF, parse the results to see if this is the correct type.
    # For a BED3/BED6 file there will be no exons, so anything other than
    # None indicates an appropriate overlap

    do something

deeptoolsintervals in the GTF/BED parsing and custom interval tree library used by deepTools. anno holds the interval tree and the associated information (gene/transcript name or ID, exons, strand, etc.). the matchType=4 is the pivotal part since that specifies that only overlaps having the exact same start position should be included. Assuming you have a stranded library, you might instead do something like switch between matchType=4 (the same start position) and matchType=5 (the same end position) depending on the orientation of the alignment and the the strand= of the gene.
Note that you'll need to make a set() to hold read names to find the mates (that or jump around the file to get them, which is probably slower).

Hey @devon-ryan, is this documented somewhere? I searched the DeepTools docs but did not find it. I am looking for other, more tuned, ways to overlap reads with features.

@fridaymeetssunday This is a semi-separate module, you can find the documentation for it here. It's shipped with deepTools but I never included it as part of the standard API documentation, since I figured it'd be internal use only.

Perfect. I was looking to see if something like the option GTF_MATCH_WITHIN existed (it does :)), though in the meantime I figure that good ol' bedtools intersect will solve my problem for now.

You can also do that with findIntersect() in R, I believe.


text:
How to convert fastq to fast5
fast5 is a variant of HDF5 the native format in which raw data from Oxford Nanopore MinION are provided. You can easily extract the reads in fast5 format into a standard fastq format, using for example poretools.
Say I have aligned these reads in fastq format to an external reference genome, resulting in a SAM file. Say I have then taken a subset of the  SAM file, according to the bitwise flag, to include only the reads that map to the reference. With the read ID, I can then grep them out from the file containing the reads in fastq format, generating a subset file in fastq format containing only the IDs that have mapped to the reference.
Now my question is, can we subset reads from the fast5 archive according to the list of mapping reads as taken from the file with reads in fastq format? This is for educational purposes, so that we have a smaller starting archive, and the fast5 -> fastq extraction takes less cpu time.

This seems a nice question to play around if we could have sample data. Could you provide (or a link to) a fast5 and the list of mappings ?

There are some sample reads from my TEDx talk last year here, reads from runs I did during the AGTA conference last year here, and reads from a chimeric read investigation here. The ONT FAST5 "standard" change so frequently that it's difficult to find something that represents what's currently available.

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 created a python program called fast5seek to do what (I think) you're after.
As you have mentioned this is for educational purposes I have added a tonne of comments to the code too so I think you shouldn't have any issues following it.
The docs on the GitHub repo have all the info, but for those reading along at home
pip3 install fast5seek
fast5seek -i /path/to/fast5s -r in.fastq in.bam in.sam -o out.txt

What it does is read in <in.fastq|in.bam|in.sam> and extract the read id from each header. It then goes through all the fast5 files under /path/to/fast5s and checks whether their read id is in the set of read ids from <in.fastq|in.bam|in.sam>. If it is, the path to the file is written to it's own line in out.txt.
If no output (-o) is given, it will write the output to stdout.
So if you wanted to pipe these paths into another program, you could do something like
mkdir subset_dir/
fast5seek -i /path/to/fast5s/ -r in.fastq | xargs cp -t subset_dir/

The above example would copy the fast5 files that are found in your fastq/BAM/SAM to subset_dir/.

As a caution, the run id should be checked (if that information is available). There's a small possibility that two FASTQ files from different runs will have the same read id.

I have now added this check for fastq files. Although, run_id is not contained within BAM/SAM files (at least not for the ones I have) so this check will only occur on fastq.

I do this very frequently, using the read file name to identify FAST5 files associated with particular reads. If a FASTQ record includes the channel number and read number (and preferably the runID as well), then I use that information to find the associated FAST5 read.
If reads are called with Albacore, then the sequencing_summary.txt file has additional information. The first column is the exact name of the file that was called, the second column is the read ID of the FASTQ sequence, and the third column is the run ID associated with that read ID. This is much easier to work with, but does require calling reads with Albacore (which, admittedly, is what seems to be producing the best results at the moment).

yes this is what I did, basically: 1. extract fast5 to fastq; 2. align the fastq to the reference; 3. get the reads ID that map from the SAM of alignment; 4. grep the mapping reads ID from the fastq; 5. extract the name of the fast5 from the fasta header of each read (last field); 6. find the corresponding fast5 in the original fast5 directory. But I was wondering if I did not miss any software that could do this on one run.

Added a note about albacore, which has a direct link between FASTQ read IDs and FAST5 file names

In light of your above comment @AminaEchchiki I have added in compatibility for BAM and SAM files too.


text:
Sleuth: transcripts with beta close to 0 are considered differentially expressed in a likelihood-ratio test
I'm comparing the results that I obtain when doing a DE analysis with the Wald test and the likelihood-ratio test. One the thing that I've noticed is that there are many genes with 'beta' close to zero that are considered differentially expressed between the conditions.
I know that the likelihood-ratio test does't use the beta values to calculate the p-values, but I find it strange that transcripts with similar expression values between the conditions are being considered differentially expressed.
Volcano plot (Wald test q-values):

Volcano plot (likelihood-ratio test q-values):

I’ve never used Sleuth but there are several odd things about this data as far as I’m concerned: the y-axis on the Wald test plot for instance seems way off (inflated). The q-values in the LR plot are a lot more reasonable.

Cross-posted on biostars and the kallisto mailing list.

What does this look like if you plot Beta vs average TPM across all samples, excluding those with qval < 0.05?

What is the actual question? Why are some genes with beta = 0 and qvalue < 0.05?

Yes, that's the thing I can't understand

This question was answered here

Could you post here also the solution to your problem? Just to avoid having to chase the anser from site to site


text:
Can gene co-expression networks be used to help identify differentially expressed genes?
In my RNAseq dataset of differentiated stem-cell lines, some samples have far fewer significantly differentially expressed genes than others. QC shows that this is because there are way fewer reads for one sample than the others.
Can gene co-expression networks be used to infer differential expression? For example, if certain genes just miss the cutoff q-value cutoff  If so, what papers have done so?

Are you using Tuxedo protocol for analysis?

How many samples do you have? It seems strange that because a single sample has fewer reads there are less DEG. Did you check for outliers? What is your normalization process?

While you can use networks to find differentially expressed genes (see the WGCNA package, which does this) in my experience this ends up largely matching what you'd get using a traditional package with a looser threshold for significance. Given the time savings of traditional packages, there's rarely any gain to using networks (it won't hurt to try, just note that it'll take some time). If some very interesting genes are just above your significance threshold then change your threshold. You need to validate your findings in some way anyway, so your p-value threshold is partly just a way to protect you from wasting your time (but if background knowledge suggests that your time wouldn't be wasted...).

Co-expression network will give you an idea about the genes having similar expression patterns and the nodes will be decided on the basis of correlation scores and nothing much you will get from differential gene expression perspective.
I would suggest you to try some other models for differential expression analysis like baySeq , DESeq, edgeR, NOIseq. Depending upon the model you use for analysis you might get some extra significant DEG. Also, try different significance cutoff and see how it's changing the analysis result to get the threshold that suits the data best.


text:
Split FASTQ and matching BAM into matching chunks
I am running a slow downstream analysis on a large set of nanopore reads (approx 3 million), and would like to split them into smaller chunks, run the analysis in massively parallel, and then recombine. Originally I just split the FASTQ into chunks, re-aligned each chunk, and then merged the output, but here I would like to use an existing alignment so I can compare results with existing analyses (i.e. the alignments must be the same).
How can I efficiently split a FASTQ file and a BAM file giving the alignment of the FASTA file into chunks, ensuring that all of the reads in FASTQ chunk 1 are in BAM chunk 1, vice versa and so on?
My FASTQ is approximately 45GB and my BAM is 33GB, so I would prefer to avoid storing one of the two files in memory if possible.
EDIT: Here's some pseudocode of exactly what I'm trying to do:
# input: in.bam, in.fastq, chunk_size
i <- 0
for fastq_read in in.fastq:
    bam_read <- extract fastq_read.read_name from in.bam
    n <- i modulo chunk_size
    write fastq_read to out.n.fastq
    write bam_read to out.n.bam
    i <- i + 1

I could swap the above to iterate through the bam file and fetch from the fastq if that's easier.

Does your BAM file already happen to have the same ordering as the fastq files (some aligners produce that by default)? If that's not the case you're pretty much left with (1) sorting files or (2) hashing read names.

Unfortunately they're not in the same order, and additionally some reads exist multiple times in the BAM (supplementary alignments) and a handful of reads do not exist in the BAM at all.

Yikes, your only option is something like what Michael Hall posted then.

Hmm, it's hard to think of a super efficient way of doing this (assuming the files aren't ordered the same - if they are then this whole answer is basically redundant). And also assuming the read ids for both files aren't a perfect intersection.
Off the top of my head you probably want to build a set of read ids for the fastq file and another for the bam. Some python code to get started on that:
import pysam
import itertools

def get_read_id_fastq(ref_path):
    """Extracts the read ids from a fastq file."""
    read_ids = set()
    with open(ref_path, 'r') as ref:
        for line in ref:
            if line.startswith('@'):  # i.e if line is header
                # split the line on spaces, take the first element, remove @
                read_id = line.split(' ')[0].replace('@', '')
                read_ids.add(read_id)

    return read_ids

def get_read_id_bam(ref_path):
    """Extract the read ids from a BAM file."""

    read_ids = set()
    with pysam.AlignmentFile(ref_path, 'r', ignore_truncation=True) as ref:
        for read in ref:
            # query_name is the query template name
            read_ids.add(read.query_name)

    return read_ids

fastq_ids = get_read_id_fastq(fastq_path)
bam_ids = get_read_id_bam(bam_path)

Then take the intersection of these two sets and that's your common read ids.
common_ids = fastq_ids & bam_ids

The next part will be a bit more involved. You will have to iterate through each file one at a time. I would suggest that for the first one you iterate through, create a running dictionary with the read id that was written as the key and the 'chuck number' it was written to as the value. You could create a cycle for your chunk size to manage this effectively
chunk_cycle = itertools.chain(*zip(range(chunk_size)*len(common_ids)))
write_idx = {}

The next part will probably require you to put in some hard-fought times getting it to work. I'll put in some rough pseudo-code to give you an idea.
for line in file:
    read_id = # get line read_id
    if read_id in common_ids:
        chunk_num = chunk_cycle.next()
        write_idx[read_id] = chunk_num
        file_to_write_to = 'out.{}.bam'.format(chunk_num)
        # open this file or write to it if already open

When you go to do the next file when you find a read id is in the common set, you look up it's value in write_idx and this will give you the chunk number to write to.
The reads wont be in the exact same order between the two files, but you could sort later if you needed it (not sure it would matter?).
Hope this helps. Sorry I couldn't give you more, but this should give you a head start hopefully.

Thanks! In fact, I wouldn't need the first part, since every read in the BAM came from the FASTQ. If there are FASTQ records that don't exist in the BAM, I don't mind.

Update: anyone interested can grab my implementation of this solution at https://github.com/scottgigante/nanopore-scripts/blob/master/split_bam_and_fasta.py


text:
How to count reads in bam per bed interval with bedtools
I recently installed Ubuntu 16.04 (because I was still using 12.04). But it seems my bedtools scripts don't work properly anymore. I can't figure out how to use the new bedtools for my old ways. What I want to do is get the number of reads from a bam file, per interval from a bed file. It was very simple with my old version of Ubuntu and bedtools:
bedtools coverage -abam file.bam -b All_peaks.bed > file.cov.txt

But this doesn't seem to work anymore. I used to get small files as results (~3MB), but now huge files are created (~4GB) including read names??
I am using bedtools v2.25 now.

The order of -a and -b switched at some point. You want:
bedtools coverage -a All_peaks.bed -b file.bam > file.cov.txt

For reference, this is the end of the help output in version 2.25:
Default Output:
     After each entry in A, reports:
       1) The number of features in B that overlapped the A interval.
       2) The number of bases in A that had non-zero coverage.
       3) The length of the entry in A.
       4) The fraction of bases in A that had non-zero coverage.

And this is the equivalent output from version 2.19:
Default Output:
     After each entry in B, reports:
       1) The number of features in A that overlapped the B interval.
       2) The number of bases in B that had non-zero coverage.
       3) The length of the entry in B.
       4) The fraction of bases in B that had non-zero coverage.

Update: The change in behavior happened in version 2.24:

We have changed the behavior of the coverage tool such that it is consistent with the other tools. Specifically, coverage is now
  computed for the intervals in the A file based on the overlaps with the B file, rather than vice versa.

Thanks Devon! I wasn't aware of the 'switch'. Switching -a and -b solves indeed my problem!


text:
Using signal peptide and the expression levels of signal recognition particle in secretome analysis
I have not found any work which investigates assessment of differences in levels of secreted proteins by taking advantage of differential expression of the genes which mediate the secretory pathway.
For example, suppose the gene whose product binds specifically to a certain signal peptide is down-regulated. I would assume that this implies a down-regulation in all the levels of secreted proteins to which it binds.
Is such an analysis possible to date? If so, where?
Thanks in advance.

One big problem is that expression levels don't always end up the same way: there are post-transcriptional changes that interfere with the final function. Finding this regulations is hard to do only using RNA-seq/expression data, and even less easy to understand what do they mean. That's why metabolomic analysis are useful, to understand what compounds are available.

This is something that Qiagen's Ingenuity Pathway Analysis is meant to solve. Unfortunately, the process of solving pathway dependency problems like these is difficult, and frequently leads to the discovery of super-linked proteins that are seen everywhere (and consequently fairly uninteresting from a biological perspective). The challenge is in filtering out the useless noise of "statistical significance", and identifying the biologically-significant signal. That's not something that computers are particularly good at doing at the moment.
My own suggestion is that you would be better off finding someone who has domain knowledge in the area of secreted proteins, and trying to present results to them in such a way that it makes it easy to interpret the results and ask further questions. If you are that person, perhaps there's a more specific question that you could ask to further investigate the problem that you have.


text:
How to confirm exon shuffling in a gene?
Note: this question has also been asked on reddit
I'm trying to confirm that the sequence of a novel gene is derived by exon shuffling between several different genes. I have the promoter sequence, gene sequence, and mRNA (with defined exon/intro boundaries). I've tried performing several database searches with the sequences, but each search produces a set of different unrelated hits. How would one go about using this info to confirm the exon shuffling hypothesis? Any tips are appreciated. Thanks!

The best proof is showing it in the wet lab. Try to clone your shuffled exon gene (cDNA made from mRNA), and sequence it the old fashion way (Sanger).

These suggestions are based on your comment that you have the gene sequence.
For a statistical approach, you could formulate a null model and show that the sequence similarity you observe is much higher than expected under the null model, but reasonable under your exon shuffling model. Ideally, you would show this individually for each exon, to convince the reader about the exons' disparate origins. I'll propose a reasonable null model, though suggestions are gladly accepted here. Pick a random gene, assume it's built using exon shuffling, and search the rest of the transcriptome for sequences that best substantiate that argument.
For an evolutionary approach, you could search for orthologs (homologues? I never remember the difference) of your gene in other species. Hopefully, the results would show that your gene formed more recently than the exons it's built out of. If your gene actually has some other origin, maybe you'll find out it's older than some of its putative components. This can't exactly "confirm" your hypothesis, but it will test your hypothesis and just be generally a cool thing to try.


text:
How to calculate p-values for fold changes?
I'm currently working with data from a Luminex multiplex assay. In this assay, the concentrations of 17 different analyte proteins were identified in 12 groups in triplicate. One of these 17 groups was used as the control, and the log2 fold changes were calculated for the analyte concentration of each sample in each group using the average control concentration for that analyte.
However, now I would like to calculate a p-value for the identified fold changes if possible.
My current preliminary idea is to perform the t test for each group compared with the control group for each analyte, but I don't think that my group sizes are large enough for this. My statistical knowledge is lacking, so if there is a better way of calculating these p-values I don't know about it.
Is there a better way of obtaining these p-values?

For future reference, basic statistics questions like this should probably go on cross-validated. This one is so basic, though, that it was faster to simply answer it here than to move it.

The key point here is whether or not the values are approximately normally distributed and whether any transformations can be applied to make them so.
For a t-test, the most important thing is that there is no relationship between the mean and the variance.
Its also correct that you don't want to divide by the mean of the controls. You want to retain all the information.
So here is how i'd start, first calculate the mean and standard deviation of each of your 18 groups (17 treatments + control), and plot this on a scatter plot (i.e. mean on X, SD on Y).
Is there a relationship between the two?
Repeat the exercise with log transformed values. Is the relationship more or less.
Pick the set of values that has the weakest mean-SD relationship. I'm guessing the log transformed might well be better.
If $X = (x_1, x_2, x_3)$ is your controls and $Y_i = (y_{i,1}, y_{i,2}, y_{i,3})$ are the treated samples from treatment i, to test the difference between treatment and control, do a t-test on X and Y (not T/mean(c))
for example in R
X = c(x_1, x_2, x_3)
Y_1 = c(y_11, y_12, y_13)
t.test (X, Y_1)

or in excel is if the three values for X are in column A and the three values for $Y_1$ in column B
=T.TEST(A1:A3,B1:B3,2,3)

If you are going to do 17 of these, remember to do a multiple testing correction by multiplying the resulting p-values by 17.
If you want to compare treatment 1 ($Y_1$) vs treatment 2 ($Y_2$) do the comparison directly, not taking the control into account:
t.test(Y_1, Y_2)

Remember that if you use log transformed values, you calculate log fold changes using mean(Y)-mean(X) rather than mean(Y)/mean(X)

Your null hypothesis would be that the fold-changes are 0, so you can either do the T-test accordingly or simply do away with the fold changes and perform the T-test between the raw values in the group (this is preferable to performing a T-test of one group vs. 0, since it allows you to assess the expected variability around 0). Note, however, that 3 samples is pretty much the bare minimum needed for any statistics, so your power will be terrible and your estimation of variance will likely be pretty inaccurate. Consequently, take any p-values with an appropriate grain of salt.
To my knowledge there aren't any good alternatives for your situation without having a good expected background distribution.

You can't calculate a p-value on the fold-change values, you need to use the  concentrations in triplicate thus giving a measure of the variance for the t-test to use.
t-test assumes your data are normally distributed, if they aren't you're going to get spurious p-values. If you aren't sure a non-parametric test like Wilcoxon is better. It will be less sensitive although with only 3 replicates your experiment is low powered anyway.

With N=3 per group a Wilcoxon won't give a p-value below 0.1 for a two-sided test or 0.05 for a 1-sided test.

@DevonRyan ah yes. Good point!

So to get a p value you must have one null hypothesis like the fold changes are significant or not. Then you test the hypothesis against the data and calculate the p value (likeliness or unlikeliness) for a given threshold.
In your case, t test is good enough to calculate the p-value.
http://willett.ece.wisc.edu/wp-uploads/2016/01/05b-TandP.pdf


text:
Which measure should be used in a PCA or RNA-seq data? TPM or counts?
I'm trying to understand the magnitude of batch effects in my RNA-seq samples, and I was wondering which expression units are more suitable to draw a PCA. I'm thinking of either counts or TPM, but things like rlog or vst could work too.
Additionally, I'm wondering whether either of this units should be log-transformed first, to avoid high-abundance transcripts driving the PCA.

tldr: log transform counts and TPMs, but rlog/vst are preferred
TPM should be log transformed to get more useful results. If you're using DESeq2 already (given the reference to rlog and vst, this seems likely), then please go ahead and use rlog or vst. That will give you more reasonable results than raw counts. If you're stuck with counts for some reason, then first use normalized counts so they're at least a bit more comparable and then log transform them so your highly expressed genes aren't driving everything.
Edit: As an aside, if you know what the batch effect is (e.g., library prep. date), it's sometimes convenient to include that in your model. You can then assess the genes that are actually changed due to that, which is sometimes handy to know (e.g., which genes might be more/less prone to degradation).

by "normalized counts" you mean something like quantile-normalization or batch-effect removal?

@mgalardini Quantile normalized, or whatever method you prefer (e.g., the ones use by default be DESeq2/edgeR/limma).

PCA works best when the input data is approximately normally distributed on each dimension. It would be a good idea to do some initial data quality checks to verify that this is the case (and transform the data appropriately if not), or at least verify that the data is approximately normally distributed in the aggregate.
For looking at Illumina RNASeq data, what worked best for me (i.e. produced the most normal-looking data) was the following steps:

Removing genes that had low raw counts in all samples
Using DESeq's variance-stabilized transform (which transforms counts into a log-like distribution)
Further normalising the VST values by dividing by the longest transcript length within each gene (which I call VSTPk)

These steps are stated in a bit more detail in our Th2 paper that was published at the end of last year:
http://jem.rupress.org/content/early/2016/12/01/jem.20160470#materials-methods

Given that I'm working with bacteria I will probably need to skip step 3. Thanks for the insight!

Do you have a reference for the claim that PCA assumes normality, because its not something I've come across before, and indeed, many sources say it is not. e.g. https://stats.stackexchange.com/questions/32105/pca-of-non-gaussian-data

I'll preface this with "I am not a statistician", and am basing this on my memories of conversations I had with a biostatistician I worked with. The operations carried out for a PCA assume that things like mean and variance and euclidean distance work normally and predictably; a grossly non-normal distribution can affect this. PCA is fairly robust to non-normal distributions, but not completely immune.


text:
Analyzing Illumina Counts
I'm pretty new to all of this--forgive me if this is a simple question.
When I download illumina counts from GEO (like the supplementary file in
GSE89225). Can I do comparisons directly on that file? Is there some normalization procedure I should go through?
Thanks!

What type of comparison?

The counts files for GSE89225 is the output of HTSeq-count as a large matrix. Unless you are developing a differential expression package yourself you should not attempt to directly use this. Rather, you should load it into R and use packages such as DESeq2, edgeR, or limma (those are the most popular ones).
For convenience, in DESeq2 you would want the DESeqDatasetFromMatrix() function after loading this.


text:
Is there a way to quickly verify the presence of some SNPs in Fastq files?
Is there a tool that can scan fastq files without assembling them for a custom list of user defined snps?

Where do fastq come from? How does the list look like? Do you have a reference genome? Please clarify.

They're bacteria sequences from a MiSeq, paired end fastq. They are all the same species and I do have a good reference genome. I'd like to provide a list of genomic coordinates I. The reference genome and find of those snps are present in the new sequences. I'm looking for an easy way just to find the small list I'd like instead of all variants from the reference.

could you edit the question and add a sample of your vcf file and the reference genome coordinates file? would be cool to play with, maybe some sed/awk one-liner will do the job

Try this pipeline:

Clean raw fastq files (trimmomatic)
Align clean fastq files to the reference genome (dna-to-dna aligner of your choice)
Convert the alignment sam file to bam (samtools)
Call SNPs from bam file (samtools)

Here's a start tutorial on this: http://ged.msu.edu/angus/tutorials-2013/snp_tutorial.html

Thanks, that's a great tutorial. I'm able to call snps and generate vcf files, but I'm looking to quickly query the vcf files for a specific list of snps, preferably listed by reference genome coordinates. I'm not sure on the wording so I haven't been able to find a tool.

You could use ARIBA, which will map you reads to the regions of your interest (more on that in a second) and then assemble/call variants.

The tool is mostly used to find known variants related to AMR (antimicrobial resistance) from a set of databases, but according to the project's wiki you can build your own database with your region of interest.
Additionally, reference sequences can be either of:

    Presence/absence sequences. ARIBA will look for these sequences in the input reads and report any that it finds, and also any variants between these sequences and the input reads.

    Variants only sequences. These should have known variant details specified in the metadata file (see below). ARIBA reports only when it finds at least one of the given variants in each of these these sequences. If a sample has one of these sequences, but does not have one of the supplied variants, then it is not reported. If you supply a variants only sequence, but no variant, then the sequence will be removed during sanity checks.

@DanK, please consider accepting one of the answers if they helped you solve your problem, thanks


text:
VerifyBamID freemix score
We all have used freemix score at some stage to check contamination or swaps in our sequencing experiment, but can anyone once for all explain how the score is calculated?

Have you read the paper? As an aside, we haven't all used VerifyBamID, it's mostly only useful in human experiments.

@DevonRyan Yes I have read the paper but could not grasp how the freemix score is estimated? Posted the same question on Biostars about a year ago but did not get any answer!!

Here's my current understanding of the freemix metric output by verifyBamID. The freemix score is a measure of contamination that is calculated without knowing the genotypes of the individuals. In the paper it is equation (2), and the authors refer to it as the "sequence-only" method. However, in addition to the sequencing reads, it requires some knowledge of SNPs and their allele frequency from the population. They demonstrate in the paper that using the SNPs and allele frequency estimates from a different population actually reduces the estimated level of contamination.
As for intuition for what equation (2) is measuring, imagine what would happen to the number of heterozygous sites detected if DNA from two different individuals was sequenced together. The total number of heterozygous sites would increase because it would count all sites where the two individuals are homozygous for opposite alleles or at least one of the two individuals is heterozygous. This would lead to an increased fraction of heterozygous sites over that expected from Hardy-Weinberg Equilibrium (i.e. 2pq).
Here's the explanation of freemix directly from the author, Hyun Min Kang:

The key idea of FREEMIX estimate is to use excessive heterozygosity to estimate the level of contamination. Especially for common SNPs, you will observe higher fraction of heterozygous alleles than 2*p*(1-p), and it turns out that you can quantify the contamination very well if you know the population allele frequency already. If you do not have accurate population allele frequency information, than it would be harder to estimate FREEMIX parameters using verifyBamID.

Additional sources:

Verifying Sample Identities - Implementation
Interpreting output files
verifyBamID Google Group


text:
Differences between TM and GDT-TS scores for structure comparison
What are the relative advantages / disadvantages of the TM (template modelling) and GDT (global distance test, total score) scores, used in protein structure prediction competitions, such as CASP? When should I choose one over the other?

Hi. This question seems to need more information to get a good answer as it looks (to me) too general. Could you please [edit] adding what have you found/read about each tool, or any paper you read on the topic to help me answer you?