text:
How to efficiently get human gene names from NCBI based on a large list of SNPs
I found a good answer related to my question here:
How to get a list of genes corresponding to the list of SNPs (rs ids)?
But it says about small number of SNPs. I want to get gene names based on approximately 12000 to 19000 SNPs. NCBI says about 10 read per second with API key. Is there any efficient way to get gene names for thousands of SNPs? Is that database available to store locally? Any advice please?

Yes, you can download the NCBI data locally, see this page.

Modify SNPs and gene annotations for your genome assembly of interest. This example shows how one might do this for hg19, for instance:

Download SNPs and convert to BED:
$ wget -qO- ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/common_all_20180423.vcf.gz \
    | gunzip -c \
    | convert2bed --input=vcf --output=bed --sort-tmpdir=${PWD} - \
    > hg19.snp151.bed

Filter this BED file for entries of interest:
$ grep -wFf snps-of-interest.txt hg19.snp151.bed | cut -f1-6 > hg19.snp151.filtered.bed6

Download gene annotations and convert to BED:
$ wget -O - ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz \
    | gunzip -c \
    | awk '($3 == "gene")' \
    | gtf2bed \
    > gencode.v19.genes.bed

Note: In the future, newer versions of Gencode annotations may be available for your reference assembly. Always check their site to confirm what is most recent.
Map genes to SNPs:
$ bedmap --echo --echo-map-id hg19.snp151.filtered.bed6 gencode.v19.genes.bed > snps_with_associated_gene_names.bed

A list of genes that overlap a SNP of interest will be in the seventh column of the output.
If this must be done within Python, you could use something like subprocess.call(...), where ... includes commands wrapped in bash -c.
I don't recommend this because of the need to escape quote marks and deal with Python's quirks in running command-line tasks, which makes a Python-based approach very difficult to set up and maintain. I'd suggest learning some basic bash scripting to make this easy and fast.

I have SNP list. Can your suggested technique be applied in Python?

Maybe look into subprocess.call() to run command line tasks inside Python.

Run them through the Ensembl VEP. You'll get the gene name and all the information about how the genes are affected by the variants.


text:
Identification of differential genes across 8 groups
I need to identify differential genes across 8 groups. I know this can be done using DESeq2 or EdgeR as these are better suited and equipped.
However, I was thinking if something like this piece of code could be tweaked. Here it works for two groups each with 3 replicates and generates a heatmap. The code is in R
pvals <- apply(dm, 1, function(x) wilcox.test(x~c(rep(0,3), rep(1,9)))$p.value)
heatmap(dm[pvals<0.05,])

Thank you for your inputs and suggestions.

What is your reason not to use DESeq2 or EdgeR? Their models will work much better for RNAseq data than a simple Wilcoxon test.

Welcome to Bioinformatics.SE! At the moment, your question simply asks "what about this?" while ignoring a tremendous amount of scholarship on differential expression analysis. Unless you can more clearly present a specific question or concern, there is a good chance this question will be closed as too broad or unclear.

For a single measurement of some fixed, observable quantity like height or grades, using a Wilcoxon test isn't a bad idea. But RNA-seq is a very complicated system and can't simply be analyzed by something like a Wilcoxon test. Your results will be full of mistakes if you do this. I would recommend that you read some introduction to RNA-seq papers or textbooks and watch some presentations by Wolfgang Huber, Gordon Smyth, and Lior Pachter to see why this is the case

source: Stack Exchange

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

sort: votes


text:
Delete all 4 lines of a fastq read from a fastq file using read ID
I have the following error when running bowtie2:

Error: Read HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG
  has more read characters than quality values.

I now want to remove all 4 lines of this specific read from the fastq file.
How can I use awk or sed to do this?

Run tail on your FASTQ file. If a job aborted while writing to disk, this read will be the last line. This error message likely means the file is truncated and missing reads. That read is not as much of a concern as the rest of the data. It’s best to run your data processing again to make sure the file is complete.

Don't do it: Your FASTQ file is either malformed or your FASTQ record spans more than four lines, which is allowed in FASTQ.  For a detailed description of what can go wrong in FASTQ parsing see for example http://biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO-module.html#FastqGeneralIterator.
If the FASTQ is malformed, then you should really ask yourself how this happened in the first place and fix the source of the problem. If the record is valid FASTQ, then I suggest parsing the read with for example FastqGeneralIterator and dumping the parsed result back to FASTQ in a 4-line-per-record form.

If you are 100% sure the read only has 4 lines (they can have more), you can use this sed command:
sed -i.bak '/^@HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG/,+3d'

The -i.bak makes sed modify the original file and create a backup copy with the same name and the extension .bak. The command just means "delete the line matching the pattern and the next three lines".

Maybe the "@" should be added to the pattern. What would happen if the read name was also on the "+" line?

When I can I like to do file processing line-by-line in the spirit of UNIX tools. You can read 4 lines from a Fastq file into 4 tab-separated values on a single line using paste, and then use grep to filter out the record in question. (Then you just have to turn the tabs back into line breaks.)
paste - - - - < reads.fastq \
    | grep -v 'HWI-D00466:116:CC62WANXX:3:1102:7363:63646' \
    | tr '\t' '\n' \
    > reads-fixed.fastq

However, it's likely that other reads in your Fastq file are corrupted as well, in which case it's probably better to write a script in Python or some other language that discards all reads whose read length doesn't match the quality string length.
But of course the most important question is: how did the data become corrupted in the first place? Where did this file come from, and what pre-processing steps were performed on it? Is getting rid of this read (or all reads with mismatched sequence/quality length) going to solve the problem for you? These aren't questions that we can answer for you, but probably deserve good answers. :-)

Here's an awk one-liner that works as intended, unlike the bioawk answer.
This removes all instances of fastq records in which the length of the seq and qual fields do not match. Again, append | gzip > filtered.fastq.gz to gzip the output.
zcat bad_file.fq.gz | awk '{ pos = NR%4; if ( pos == 1) {h=$0} else if (pos == 2) \
  {s=$0} else if (pos==3){c=$0} else if ( pos==0 ) \
  {q=$0; if ( length(q) == length(s) ){printf("%s\n%s\n%s\n%s\n", h, s, c, q)}}}'

Delete the \ characters at the end of the first two lines to make it easy to paste into your terminal.
Here is a shorter version of the same thing thanks to @terdon
zcat file.fastq.gz | awk '{ record[++k]=$0; if(NR%4==2){slen=length($0)} else if(NR%4==0){if(length($0) == slen){ for(i in record){print record[i]}} k=0; }}'

Heh, I was just writing one too :) Here's what I came up with, it's a bit shorter than yours: zcat file.fastq.gz | awk '{ record[++k]=$0; if(NR%4==2){slen=length($0)} else if(NR%4==0){if(length($0) == slen){ for(i in record){print record[i]}} k=0; }}' feel free to add it to your answer if you want to.

To generalize this to remove all reads that do not have the same length of read and sequence quality, you can use this bioawk one-liner:
bioawk -cfastx '{ if (length($seq) == length($qual) ){printf("@%s %s\n%s\n+\n%s\n", $name, $comment, $seq, $qual)} }' my_fastq.gz

This strips anything in the "fastq comment field", +, but that almost never has anything there anyway.
Append | gzip > filtered.fastq.gz to the above command to gzip the output.

This doesn't seem to work for me (Arch Linux, bioawk version 20110810). Even after correcting the syntax error (corrected in your answer), it seems to exit after the first record it finds where the lengths are not the same. Can you confirm it works on your machine?

Thanks for the edit! Yep, I just confirmed that it works for me on all records. I'm using Ubuntu 16.04. >>bioawk --version returns awk version 20110810. It looks like the last time that I compiled bioawk was from the github repo in 2016.

That's really weird. I have the same version. Can you try a case where the record with different lengths is the second one in the file? Does it not just print the first record there?

I got the same error as you - I made the second entry in the fastq file have a different length of seq and qual fields and bioawk just exited before printing that record. I wonder if this is intended behavior? I have to add that in my above comment, I just tried the bioawk one-liner on a normal fastq file without any length mismatches. Oops!

I don't see how it can be intended. Sounds like a bug in bioawk.

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


text:
Hisat2 : which option should mention for strand specific library read
I was trying hisat2 I get confused by the strand options...
The question has been already asked on github here but didn't get any satisfactory answer.
I'm surprise how much this strand information could be confusing in different tools.
Looking at the hisat2 manual I get confuse by the two different options:
--rna-strandness

Specify strand-specific information: the default is unstranded.
  For single-end reads, use F or R. 'F' means a read corresponds to a transcript. 'R' means a read corresponds to the reverse complemented counterpart of a transcript. For paired-end reads, use either FR or RF.
  With this option being used, every read alignment will have an XS attribute tag: '+' means a read belongs to a transcript on '+' strand of genome. '-' means a read belongs to a transcript on '-' strand of genome.

and --fr/--rf/--ff

The upstream/downstream mate orientations for a valid paired-end alignment against the forward reference strand. E.g., if --fr is specified and there is a candidate paired-end alignment where mate 1 appears upstream of the reverse complement of mate 2 and the fragment length constraints (-I and -X) are met, that alignment is valid. Also, if mate 2 appears upstream of the reverse complement of mate 1 and all other constraints are met, that too is valid. --rf likewise requires that an upstream mate1 be reverse-complemented and a downstream mate2 be forward-oriented. --ff requires both an upstream mate 1 and a downstream mate 2 to be forward-oriented. Default: --fr (appropriate for Illumina's Paired-end Sequencing Assay).

Should both need to be specified ? Are they redundant ?
I didn't find so many person struggling with that, so I guess there is something obvious that I don't get. I hope someone may enlighten me about it !

--fr/--rf/--ff should rarely be set, since they refer to the relative orientation of reads and basically everything these days is --fr. Mate-pairs were --rf, but you don't see those much any more.
--rna-strandedness is a very different option, since it sets how reads are expected to align against genes. Most stranded protocols in use these days follow the dUTP-method, where read #2 in a pair has the same orientation as the transcript from which it arose. So either R or RF would typically be appropriate, unless the library is unstranded. In practice, I expect this is more useful if you plan to run stringTie downstream, since then the XS auxiliary tag is set appropriately.

I found this resource t visualise the relative orientation of reads. The --fr/--rf/--ff option is kind of disturbing... you say Mate-pairs were --rf, but you don't see those much any more. but In a "old" tool like tophat rf was not possible, when in a much recent tool they have introduced it...

There are some other sequencing types (i.e., not RNAseq) where different library types are now being produced. Presumably they wanted to allow some flexibility so hisat2 could be used with them as well.

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


text:
Filtering step for read counts data
I have around 1200 samples as columns and 60,000 genes with Htseq-Counts data. Before normalization with voom function I want to do filtering step.
I want to remove genes whose expression is == 0 in at least 10 samples.
Can I do this with read counts itself or should I transform counts to cpm?
If filtering with counts data is fine May I know how to do this.

If you have your counts in a data.frame called counts, something like this might work:
filtered.counts <- counts[rowSums(counts==0)<10, ]

For example lets assume the following data frame.
> A <- c(0,0,0,0,0)
> B <- c(0,1,0,0,0)
> C <- c(0,2,0,2,0)
> D <- c(0,5,1,1,2)
>
> counts <- data.frame(A=A, B=B, C=C, D=D)
> counts
  A B C D
1 0 0 0 0
2 0 1 2 5
3 0 0 0 1
4 0 0 2 1
5 0 0 0 2

And now you would want to remove the genes with at least 3 zeroes or more:
> filtered.counts <- counts[rowSums(counts==0)<3, ]
>
> filtered.counts
  A B C D
2 0 1 2 5
4 0 0 2 1

Every gene with 3 or more zeroes is now removed.

In at least 10 samples. So I guess it should be this filtered.counts <- counts[rowSums(counts==0)>=10, ]

No, you keep all rows with less than 10 zeroes with my example.

In your way I will be having only 1/4 th of genes remaining. I want to remove genes whose expression is == 0 in at least 10 samples. With my command now I have more than half genes remaining after filtering.

If you use it your way, you keep all the genes with more than 10 zeroes, check it for your self if you don't believe me.

So, the filtered.counts output that u gave above should be used for normalization?

If this filtering approach is what you want yes. However, in edgeR or limma manual they also give other options, such as at least 1 cpm in more than 3 samples. Usually 3 samples is then based on groups of 3 samples. If you have only 1/4 of your genes left with this filtering it might be too stringent. For limma voom, you need to have cpm values btw.

Ok. As you said I did it in edgeR way.
keep <- rowSums(cpm(counts) > 0.5) >= 10; y <- counts[keep,] So, now genes in "y" can be used for "voom" function to normalize the data. Do you think this is right?

Yes I think it is right, good luck!

10 out of 1200 samples is a very strict cut off. I'm not surprised you don't have many genes remaining.

Yes, but I'm interested in looking for a specific gene. So, to get that in my results I took that cutoff.

While this answers explains how to do it I want to address when and why and which thresholds to do it.
Filtering the genes with low counts is usually done because the counts are not reliable it would be noise, specially when there are low number of samples these genes disturb power of the analysis. However with 1200 samples having 10 or more samples without any count is valuable information because you have enough samples to estimate if this is due to a latent or known variable or if it is due to noise.
The CPM are used because the size of the library per sample affects the power of the experiment. The higher the library size the easier is to find if a gene is expressed or not. CPM are calculated to take into account that effect. Filtering by counts or CPM will result in different genes filtered if you use other thresholds than 0 (for instance I have filtered those genes that the average CPM where below 1, which would be different to filtering by those that the average counts is below 1.)

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


text:
Extracting expression data from GSE dataset downloaded from GEO
I have downloaded GSE16146 dataset from GEO using GEOquery R package.
I would like to extract "Data table" from downloaded GSE16146.
>library("GEOquery")
>GSE16146 <- getGEO("GSE16146")
>Table(GSE16146)[1:5,]

This returns the following error:
> Table(GSE16146)[1:5,]
Error in (function (classes, fdef, mtable)  :
unable to find an inherited method for function ‘Table’ for signature ‘"list"’

Can anyone help?

Have you tried the function getGSEDataTables to access GSE Data Tables from GEO, as reported in this tutorial? Also, what version of R/GEOquery are you using?

GEOquery version 2.42.0, R version 3.4.0  > getGSEDataTables(gset)
Error in as.character.default(<S4 object of class "GSE">) :
no method for coercing this S4 class to a vector

The error is telling you that you can't calculate a Table for a list. Also you might have tried to use table instead of Table.

According to the manual, all you need to do is:
library('GEOquery')
gseGSE16146 <- getGEO('GSE16146', GSEMatrix=FALSE)

As explanation, getGEO() outputs by default to GSEMatrix=TRUE and returns a list of ExpressionSet objects. You should get what you were looking for with:
Table(GSMList(gseGSE16146)[[1]])[1:5,]

The manual has also a paragraph about this:

The GSE has  a  metadata section, just like the other classes.
  However, it doesn’t have a GEODataTable.  Instead, it contains two
  lists, accessible using GPLList and GSMList, that are each lists of
  GPL and GSM objects.

I tested the code I posted on R 3.3.2 and I get access to the Data Table just fine.
update:
Answering how to extract the actual expression data (not what asked in the original question):
library(GEOquery)
data = getGEO("GSE16146")
datExpr = exprs (data[[1]])

Thanks. It works. Would you suggest how can I extract all the expression matrix for 80 samples (list) from the gseGSE16146?

Hi thanks for your answer! What I didn't get is how do I access the actual expression data?

@Tapper I updated my answer above with how to extract the expression data from the experiment. hope it helps?

@aechchiki Thank you! The snippet does not contain the data but a list in my case though: GSM4446615 GSM4446616 GSM4446617 GSM4446618. The file getGEO() got was anyway to small to contain the dataset imho. I finally got the data by downloading the big data file myself and processing it with: gse = read.table("GSE123456.txt", sep = "\t"); data=CreateSeuratObject(raw.data = gse, ..) -- What do you think?


text:
Difference between BWA-backtrack and BWA-MEM
Many of my colleagues recommend I use BWA-MEM instead of regular old BWA. The problem is I don't understand why and reading the BWA man page doesn't seem to help the matter.
What is the difference between BWA and BWA-MEM? And, in which instances would you employ one over the other?

You should read the BWA and BWA mem publications.

This paper: https://arxiv.org/abs/1303.3997

Possible duplicate: https://bioinformatics.stackexchange.com/questions/418/why-is-bwa-mem-the-standard-algorithm-when-using-bwa

TL;DR:
BWA-backtrack is based on backtracking.
This approach is appropriate only when the dissimilarity between the reads and the reference is low,
or when you want to find all best hits or enumerate all possible alignments up to a specified
number of errors.
In all other situations, BWA-MEM is preferable as it can,
thanks to its sophisticated strategy based on maximum exact matches,
deal with errors better and also automatically switch between local and
global alignment modes.
Long description:
I would like to provide some algorithmic insight since I believe that it
can be very useful in this case.
Both BWA-backtrack and BWA-MEM use the same indexing strategy
(heavily relying on BWT-index), but the actual algorithms are quite different.
BWT-index (and also other full-text indexes such as
suffix arrays or suffix trees)
can easily find exact matches (imagine Ctrl+F-like search in a text editor),
but any differences between the read and the reference,
such as sequencing errors or genomic variants,
make the situation complicated. One then needs to
somehow transform inexact matching to exact matching,
and all three BWA mappers (note that there exists also BWA-SW, but it is deprecated)
use quite different strategies.
BWA-backtrack looks for substrings of the reference, which would be similar
to the entire read (end-to-end) using an algorithm
called backtracking.
First, it searches occurences of the read without any "corrections".
If nothing found, it consideres all possible single edits; then two edits, etc.
To make mapping efficient, one usually wants to stop with the first
found alignment as this would be the best one.
When required, it is also possible to find the other equally good alignments
or to enumerate all alignments up to some edit distance or withing some divergance rate
(see the -N option of BWA-backtrack).
It turns out that
the time required for finding an alignment
can be exponential in the number of errors, which is probably the main problem of
backtracking-based approaches.
To prevent huge overheads due to dissimilar reads,
one needs to limit the number of allowed errors to some reasonable number
(see maxDiff in the BWA man page)
and consider the other reads unaligned.
In the case of BWA-backtrack,
the minimum required identity level is ~97% with the default options (see the -n option).
In fact, the algorithm is more complicated and uses various heuristics such as
seed-and-extend or Z-dropoff in order to make the computation fast enough (at the price of lower accuracy). If you are interested
in more details, all these tricks are well described in the paper.
BWA-MEM uses quite a different strategy. It detects long exact matches
between the read and the reference, and then chains them into local or global alignments,
based on what is more appropriate in that specific case.
Such an automatic local-global switching can be very powerful and
BWA-MEM works well with various types of data
(short reads, long reads, low error rates, high error rates, etc.).

To quote the Introduction to BWA on sourceforge:

BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome. It consists of three algorithms: BWA-backtrack, BWA-SW and BWA-MEM. The first algorithm is designed for Illumina sequence reads up to 100bp, while the rest two for longer sequences ranged from 70bp to 1Mbp. BWA-MEM and BWA-SW share similar features such as long-read support and split alignment, but BWA-MEM, which is the latest, is generally recommended for high-quality queries as it is faster and more accurate. BWA-MEM also has better performance than BWA-backtrack for 70-100bp Illumina reads.

In short, for anything where you have read lengths over 70bp BWA-MEM is faster, and more accurate.

One important addition is that BWA-MEM will split reads more aggressively. Whereas BWA-ALN requires (if I remember correctly) one end-to-end alignment for anchoring, BWA-MEM allows for partial alignments.
Practicable, use the -M flag for bwa mem if you are processing the alignments with older or less mainstream programs if they don't understand the BAM supplementary read flag yet.

@Manuel: What specific programs do you have in mind that need -M?  That flag was added for use with Picard < 1.96, which is now ancient — Picard itself hasn't benefited from or needed this flag for four years now.  See also https://github.com/lh3/bwa/pull/26


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

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

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

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

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

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

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

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

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

@IanSudbery thanks for the tip, will do that

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

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

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

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

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

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

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

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

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

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

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

for alignment in bamfile.fetch():

    if alignment.is_unmapped:
        continue

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


text:
PASA pipeline: compare experimental transcripts to the reference annotation
I would like to ask if anyone has experience in running a subset of the PASA pipeline, in particular for the reconciliation of some experimental 'transcripts' with the reference annotation.
In more detail, I am working with RNA-seq data from D. melanogaster. I have reconstructed the 'transcripts' using Trinity. I have aligned these 'transcripts' to the reference genome using GMAP. Now I would like to match these 'experimental transcripts' with the reference annotation, to see how they compare.
I was wondering if I can run just a subset of the whole PASA pipeline, so basically skipping the step of alignment to the reference and using the alignment file that I generated externally.
Also, I would be glad to receive more general redirecting to any released software allowing the comparison of experimental transcripts to the reference annotation.
Cross-posted on the PASA Google Group pasapipeline-users
Cross-posted on biostars

Go easy on the cross-posting. You've already directly asked the authors, just be patient for their reply.

Hi Amina Echhiki, thanks for your question and welcome to Bioinformatics Stack Exchange. As you've included this (or a similar question) on other sites, it would be appreciated if you could report on the answers you have received as your own answer to this question. Stack Exchange encourages people to answer their own questions as they gain knowledge about the subject.

Answer from the authors:

PASA ends up using gff3 format instead of sam format for uploading the
  alignments.  If somebody has the gmap gff3 output, then can be done to
  just upload that directly using the custom alignment importer, but not
  with sam. They suggest just having PASA rerun GMAP as part of its
  regular routine.


text:
tools to reconcile experimental transcripts with reference annotation
Looking for tools to reconcile alignment file of experimental transcripts mapped to genome (SAM/BAM) with the reference transcriptome annotation (GTF) from Ensembl (organism: D. melanogaster).
The aim would be to check which transcripts reported in the alignment file are also reported in the reference annotation, generating a new annotation (e.g., GTF) including the new transcripts as seen in the alignment file (experimental transcripts to genome).
Any idea?
-- Edit:
Up to now, I tried an in-house script that does the following:

generate a BED file from the BAM file with the exon coordinates (start/end) per transcript - using BEDOPS v2p4p27
convert the reference GTF to BED to retrieve the coordinates per annotated feature
compare the experimental vs the reference BED intervals
infer the suite of annotated exons per experimental transcript, thus, ideally, the isoform it corresponds to

However, with this simple approach, I only get around 20% of the experimental transcripts that can be identified in all their length using reference annotation. This might be due either to a too simple approach (e.g. if the start/end of the exon in the experimental transcript is shifted by 1 base as compared to the reference, it won't be detected), or to a real biological signal (since I am also using data from long-read sequencing technologies).

Can you convert the GTF file (from the BAM file via BED) to GTF? Then you can use cuffcompare. I've never tried making a GTF out of BAM files, so I haven't a clue if that works.

I will try that, and post an answer here if I get to a solution

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

I think that's exactly the solution we are converging to ;)

As per my answer to @_julien_roux on twitter:
Trying to find novel transcripts within the context of an existing annotation is much less straightforward. You probably need to do a "genome-guided assembly" with Trinity and PASA:
http://pasapipeline.github.io/#A_ComprehensiveTranscriptome
We did something similar in much simpler organisms in our recent paper:
http://dx.doi.org/10.1186/s12864-017-3505-0
Your biggest problems are going to be alternative spliced gene models, the large amount of fragmentation inherent in de novo transcriptome assembly and largely underreported non-coding RNAs in most genomes. It was easier for us as the dictyostelids have few and small introns - fly is going to be a different ball-game. Be prepared to spend a lot of time optimising parameters before finding a result you're happy with (maybe) ;)
You may think the fly genome pretty good and complete, but be prepared to find lots of errors, corrections and novelty. We were very surprised at the amount of new annotation we were able to find in D. discoideum - an organism where every single gene had been manually curated!
Good luck!

Sure, we know this is not a simple project ;) We now have the reconstructed transcripts, from different methods and types of data and would like to evaluate them by comparing to reference annotation (one out many ways we are trying to evaluate these)

@jul635 ah ok. Then I'd strongly recommend PASA, it gives an interactive report of all the changes made to an existing annotation by your transcripts. I'd also recommend Transrate as great tool for quality assessing your assembled transcripts.

Thanks to @Devon_Ryan's suggestion, I had a go at using GMAP with our own de-novo assembly. Here's what I started with:

Genome assembled from nanopore long reads, corrected using Illumina cDNA reads
genome-guided Trinity transcriptome assembly, with reads grouped based on mapping to the assembled genome, then assembled by Trinity

The transcriptome assembly was filtered to only include transcripts that had reasonable coverage from cDNA reads, and further filtered to cluster similar transcripts and only use the longest isoform. This process is detailed in our paper. If I were to do that again, I'd probably create a supertranscript containing all exons, rather than the longest isoform.
I first created an index for the genome using gmap_build:
gmap_build -k 15 -j 8 -d 'MINBRA' \
  Nb_ONTCFED_65bpTrim_t1.contigs.hpcleaned.fasta -D MINBRA_GMAP

I then used gmap to map the transcriptome to the genome, generating gene models as a GFF3 file:
gmap -t 10 -d MINBRA -D ./MINBRA_GMAP \
  -A locusTagged_CDLI_CD98LMOHC50_TBNOCFED.tran.fasta \
  -f 2 > Nb_TrinTran_vs_Canu.gff3

Here's an example from the annotation file that was produced by GMAP:
...
tig00022147 MINBRA  gene    413327  414778  .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.path1;Name=MINBRA_TGG_483_c1_g1_i1
tig00022147 MINBRA  mRNA    413327  414778  .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.path1;coverage=100.0;identity=99.9;matches=1452;mismatches=0;indels=2;unknowns=0
tig00022147 MINBRA  exon    413327  414778  99  +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna1.exon1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna1;Target=MINBRA_TGG_483_c1_g1_i1 1454 1 .
tig00022147 MINBRA  CDS 413417  413557  100 +   0   ID=MINBRA_TGG_483_c1_g1_i1.mrna1.cds1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna1;Target=MINBRA_TGG_483_c1_g1_i1 231 91 .
###
tig00001947 MINBRA  gene    37526   38961   .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.path2;Name=MINBRA_TGG_483_c1_g1_i1
tig00001947 MINBRA  mRNA    37526   38961   .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna2;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.path2;coverage=100.0;identity=98.7;matches=1435;mismatches=1;indels=18;unknowns=0
tig00001947 MINBRA  exon    37526   38961   98  +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna2.exon1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna2;Target=MINBRA_TGG_483_c1_g1_i1 1454 1 .
tig00001947 MINBRA  CDS 37614   37750   97  +   0   ID=MINBRA_TGG_483_c1_g1_i1.mrna2.cds1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna2;Target=MINBRA_TGG_483_c1_g1_i1 231 91 .
###
tig00022147 MINBRA  gene    464550  465456  .   +   .   ID=MINBRA_TGG_489_c13_g1_i2.path1;Name=MINBRA_TGG_489_c13_g1_i2
tig00022147 MINBRA  mRNA    464550  465456  .   +   .   ID=MINBRA_TGG_489_c13_g1_i2.mrna1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.path1;coverage=100.0;identity=99.3;matches=907;mismatches=0;indels=6;unknowns=0
tig00022147 MINBRA  exon    464550  465456  99  +   .   ID=MINBRA_TGG_489_c13_g1_i2.mrna1.exon1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.mrna1;Target=MINBRA_TGG_489_c13_g1_i2 913 1 .
tig00022147 MINBRA  CDS 464909  465053  98  +   0   ID=MINBRA_TGG_489_c13_g1_i2.mrna1.cds1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.mrna1;Target=MINBRA_TGG_489_c13_g1_i2 507 361 .
###
...

GMAP told me on standard output what transcripts were not mappable to the genome (136 of 36830, or about 0.4%).
I've now used STAR to map these cDNA reads to the transcripts (using the GFF3 annotation file for initial evidence) in order to confirm that the mapping carried out by GMAP is reasonable. I ended up throwing out another 7.4% of transcripts. First, I generated a genome index using the generated gff3 file:
~/install/star/STAR-2.5.4b/bin/Linux_x86_64/STAR --runThreadN 10 \
   --runMode genomeGenerate --genomeDir ./star_genome \
   --genomeFastaFiles Nb_ONTCFED_65bpTrim_t1.contigs.hpcleaned.fasta \
   --sjdbGTFfile Nb_TrinTran_vs_Canu.gff3 \
   --sjdbGTFtagExonParentTranscript Parent

Followed by mapping my RNASeq experiments using STAR. STAR can sort the genome mapping, but not the transcriptome mapping, so an additional sort/index was done on that:
for x in $(seq 1 9); do
  ~/install/star/STAR-2.5.4b/bin/Linux_x86_64/STAR --runThreadN 10 \
     --readFilesCommand zcat --genomeDir ./star_genome \
     --readFilesIn *-1563-0${x}-*.fastq.gz \
     --outFileNamePrefix mapped/1563-0${x}_vs_BNOCFED \
     --outSAMtype BAM SortedByCoordinate --quantMode TranscriptomeSAM
  samtools sort mapped/1563-0${x}.toTranscriptome.*.bam \
    > mapped/1563-0${x}.toTranscriptome.sortedByCoord.out.bam
  samtools index mapped/1563-0${x}.toTranscriptome.sortedByCoord.out.bam
done

BAM hits were collated into a single file:
(for x in $(seq 1 9);
    do samtools idxstats 1563-0${x}.toTranscriptome.sortedByCoord.out.bam | \
      perl -pe "s/^/1563-0${x}\t/";
 done) | gzip  > readCounts_all_vs_transcriptome.tsv.gz

Which was then processed using R to identify unmapped transcript annotations:
#!/usr/bin/Rscript
data.df <-
    read.delim("readCounts_all_vs_transcriptome.tsv.gz", header=FALSE,
               col.names=c("Sample","Transcript",
                           "Length","Mapped","Unmapped"));
gene.counts <- tapply(data.df$Mapped + data.df$Unmapped,
                      list(data.df$Transcript, data.df$Sample), sum);
sum(rowSums(gene.counts)==0) / nrow(gene.counts);
TranscriptsToRemove <- rownames(gene.counts)[rowSums(gene.counts) == 0];
cat(TranscriptsToRemove, file="transcripts_to_remove.txt", sep="\n");

The unmapped annotations were then excluded from the gff3 file:
grep -v -F -f mapped/transcripts_to_remove.txt Nb_TrinTran_vs_Canu.gff3 \
  > mappable_Nb_TrinTran_vs_Canu.gff3


text:
How to convert GFF3 to GTF2
I would like to convert a file in gff3 format to a gtf2.2 format.
The reason why I would like to do this is: I have a set of transcripts assembled by a bunch of different software (and using RNA-seq data from different sequencing technologies) and I would like to compare each of those sets with the reference transcriptome annotation for that species (D. melanogaster). I already asked for the community advice about how to proceed on that, but to run the suggested software (TACO and cuffmerge) I need to have a GTF file containing the experimental transcripts to compare to the reference transcripts (GTF).
Up to now, I unsuccessfully tried:

gffread utility in the Cufflinks package (gffread input.gff3 -T -o output.gtf): this results in an empty output.gtf file and an empty log file (used Cufflinks v.2.2.1) - I contacted the authors via their Google group but haven't heard of them yet
gff3_to_gtf utility in the GenomeTools (gt) package (gt gff3_to_gtf input.gff3 -o output.gtf): the output is not created and the log file is not informative - I contacted the authors via their mailing list but haven't heard of them yet
GFF3_to_GTF utility in the FML package (./gff3_to_gtf_converter.pl input.gff3 output.gtf): the output just contains a header (##gff-version 2.5) and the log is empty

The gff3 file was created as output of GMAP, and contains the transcripts as found by alignment to the reference (specifying option -f gff3_match_cdna).
[Edit: what I found out a poteriori, is that this format is not a standard gff3, thus the conversion is not trivial...]
Here is the beginning of the gff3 file I tried to convert:
$ head r9_gmap_match-cdna.gff
##gff-version   3
# Generated by GMAP version 2017-04-24 using call:  /home/aechchik/software/gmap-2017-04-24/bin/gmap.sse42 -d gmapidx -D /scratch/beegfs/monthly/aechchik/isoforms/ref/chromosomes/ -f gff3_match_cdna -n 0 -t 20 r9_2d.fasta
2L    gmapidx    cDNA_match    18442664    18443024    79    -    . ID=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1;Name=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d;Target=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d 141 489;Gap=M13 I10 M8 D3 M8 D3 M32 D1 M4 D2 M5 D3 M34 I3 M6 D1 M7 D1 M5 D1 M30 D2 M16 I1 M8 D3 M10 D1 M5 D1 M6 I1 M8 I1 M8 D1 M33 D3 M33 I1 M28 D3 M25
###
3R    gmapidx    cDNA_match    15853880    15855465    96    +    . ID=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1;Name=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d;Target=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d 80 1645;Gap=M46 D2 M12 D2 M157 D3 M4 D1 M50 I2 M3 I1 M12 I1 M68 I1 M66 I1 M53 D2 M47 D1 M16 D1 M35 D1 M155 D1 M28 D1 M166 D2 M47 D1 M8 D4 M69 D1 M28 D1 M5 D1 M12 D1 M202 I1 M115 D1 M61 I2 M7 D1 M7 I1 M36 D2 M41
3R    gmapidx    cDNA_match    15855529    15855742    97    +    . ID=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1;Name=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d;Target=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d 1646 1856;Gap=M21 D1 M68 D1 M85 D1 M37
###
X    gmapidx    cDNA_match    14837810    14838142    93    -    . ID=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1;Name=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d;Target=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d 74 406;Gap=M13 D1 M182 I1 M14 I2 M56 I1 M30 D2 M21 D1 M13
X    gmapidx    cDNA_match    14837470    14837753    92    -    . ID=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1;Name=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d;Target=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d 407 688;Gap=I2 M5 I1 M64 I1 M44 D1 M9 D5 M31 D3 M23 I1 M19 I1 M25 I1 M26 D1 M19 I1 M9
###

I can't get cufflinks to work with this either. Strange. Could you explain what you need, exactly? I mean, GTF is basically GFF2. Would it be OK to just remove the attributes? What information do you need to keep? Basically, what will you use the GTF for? Do you need the gene_id and transcript_id fields and, if so, where would that information come from? It might be helpful if you came into our chat room. If you do, ping me (@terdon) to let me know.

I tried to convert this with rtracklayer in R and also failed (rtracklayer doesn't want to output proper GTF format) :( I wonder if a bit of python/perl would suffice.

The following bit of python code should work:
#!/usr/bin/env python
import sys

lastTranscript = [None, None, None, []]  # ID, chrom, strand, [(start, end, score), ...]

def getID(s):
    """Parse out the ID attribute"""
    s = s.split(";")
    for k in s:
        if k.startswith("ID="):
            return k[3:]
    return None

def dumpLastTranscript():
    """Print the last transcript"""
    bounds = sorted(lastTranscript[3])
    print("{}\tgmapidx\tgene\t{}\t{}\t.\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], bounds[0][0], bounds[-1][1], lastTranscript[2], lastTranscript[0], lastTranscript[0]))
    print("{}\tgmapidx\ttranscript\t{}\t{}\t.\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], bounds[0][0], bounds[-1][1], lastTranscript[2], lastTranscript[0], lastTranscript[0]))
    for start, end, score in bounds:
        print("{}\tgmapidx\texon\t{}\t{}\t{}\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], start, end, score, lastTranscript[2], lastTranscript[0], lastTranscript[0]))

def handleLine(cols):
    """Handle a single line, appending the exon bounds to the previous if relevant"""
    ID = getID(cols[8])
    assert(ID is not None)
    if lastTranscript[0] is not None and lastTranscript[0] != ID:
        dumpLastTranscript()
        lastTranscript[3] = []
    lastTranscript[0] = ID
    lastTranscript[1] = cols[0]
    lastTranscript[2] = cols[6]
    lastTranscript[3].append((int(cols[3]), int(cols[4]), cols[5]))

f = open(sys.argv[1])
for line in f:
    if line.startswith("#"):
        continue
    cols = line.strip().split("\t")
    handleLine(cols)

dumpLastTranscript()
f.close()

If you saved that as gff2gtf.py and made it executable, the usage would be ./gff2gtf.py foo.gff3 > foo.gtf. With the example that you provided in your post, the result is:
2L  gmapidx gene    18442664    18443024    .   -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
2L  gmapidx transcript  18442664    18443024    .   -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
2L  gmapidx exon    18442664    18443024    79  -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
3R  gmapidx gene    15853880    15855742    .   +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx transcript  15853880    15855742    .   +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx exon    15853880    15855465    96  +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx exon    15855529    15855742    97  +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
X   gmapidx gene    14837470    14838142    .   -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx transcript  14837470    14838142    .   -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx exon    14837470    14837753    92  -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx exon    14837810    14838142    93  -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";

That at least looks like correct gtf 2.2 format to me.

Close to be GTF2.2 but is not. Feature type from 3rd column should be filtered to match the GTF2.2 specification. Only 9 feature types are accepted. Gene is not as example. It became accepted in version 2.5. See here https://agat.readthedocs.io/en/latest/gxf.html#gtf2-2

@juke34 As there is no coherent GFF specification, you're always going to have to tweak something or other.

There are clear specifications (I would say excepted for GFF0 and GTF1) but it is easy to not follow them ^^. Have you look at the review I made? I tried to gather all the information in one place:  https://agat.readthedocs.io/en/latest/gxf.html

The GenomeTools package has lots of small applications that let you massage GFF3. It may have a converter, or it may let you tweak your existing GFF3 so that gffread works

Oh... I should read past the first item on the list... Try using gt gff3 to transform and output GFF3 files?

yes @Dan, I tried gt (Genome Tools package) - see point 2 of my question, and it is not working as it is.

Please edit your question and include an example of how the OP can use your suggested tool to do what they ask. As it stands, this is more of a comment than an answer. Please also see https://meta.stackexchange.com/q/225370/203101

@terdon Do you mean the answer?

@KamilSJaron d'oh! Yes, sorry, I was asking Dan to edit this answer, not the question.

I'm quite sure it will work with agat_convert_sp_gff2gtf.pl from AGAT.
See the review of the different solution I made here:
https://agat.readthedocs.io/en/latest/gff_to_gtf.html
AGAT it the only one you can decide wich version of GTF you require as output.


text:
STAR-long parameters for aligning RNA ONT reads to genome
Are there any suggested parameters to align ONT reads to the reference genome using STAR-long? For now, I used the parameters suggested here, but I noticed a weird behaviour.
I have RNA reads (D. melanogaster) from R7 and R9 flowcells, separately. I only selected to analyze 2D reads in pass category.
I have, respectively, 113249 reads for R7 and 40318 reads for R9. I aligned those reads and get (only!) 150 uniquely mapped reads for R7 data and 8017 uniquely mapped reads for R9 data. I tried to run again the same command on a different server with fresh compilation, but the output file is consistent with these 150 reads.
However, if I align the same with GMAP, I get 78016 uniquely mapped reads for R7 and 33523 uniquely mapped reads for R9, so I suspect that something went wrong in the alignment run.
I am aware that the two mappers behave very differently, STAR-long being  more precise and preferring to report mappings of fewer reads but at better loci, and GMAP being overall less precise, trying to map the most of the reads but at not-so-good loci.
I was wondering if some of you had experience with this and could suggest me the best parameters for RNA reads from ONT?

I guess my reply would be more appropriate as a comment - but I can't. I'm not sure if that setting of SE makes sense, but whatever. This is not necessarily a solution for the behaviour you see, but I would suggest trying https://github.com/lh3/minimap2 for cDNA-seq alignment.

At my hand, STAR-long doesn't work well with noisy reads. It worked for iso-seq mostly because the error rate of iso-seq reads is much lower.

I've had great results using minimap2, particularly when combined with a pre-treatment of Canu for error correction (using minimap2 for the read-to-read mapping):
# correct reads
~/install/canu/canu-1.6/Linux-amd64/bin/canu overlapper=minimap \
   genomeSize=100M minReadLength=100 minOverlapLength=30 -correct \
   -p 4T1_BC06 -d 4T1_BC06 \
   -nanopore-raw workspace/pass/barcode06/fastq_runid_*.fastq
# align reads
~/install/minimap2/minimap2 -p 10 \
   -a ~/db/fasta/mmus/ucsc/mmus_ucsc_all_cdna.idx \
   -x splice <(pv all/4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  samtools sort > 4T1_BC06_all_vs_mmusAll.bam

Update: the most recent nanopore basecaller, combined with the most recent minimap2, seems to now do a decent job with mapping, and pre-correction of reads no longer seems necessary. More recently I've been using LAST to map to the transcriptome, and minimap2 to map to the genome. Minimap2 has the ability to use a homopolymer-compressed genome index, which means that the most common consistent error for nanopore reads (i.e. misjudgement of homopolymer length) will not influence the mapping rate.

Yep I also ended up using minimap2, thanks. At the time of this question, I was just trying to check if I was doing some nonsense with parameters for STAR, or if its behavior was simply due to its algorithm. For reads correction, aren't you afraid of altering the indel composition of your reads if you use Canu?

No... why would that be a problem? Canu is correcting errors based on the consensus of overlapping reads; I don't see why/how that correction would cause more problems than it solves.


text:
Exon-exon junctions: compare experimental transcripts to reference annotation
My aim is to parse an experimental transcript set (obtained by RNAseq) to check which splice junctions are already reported in a reference annotation and which ones are new.
I tried Regtools:

feed the experimental alignment file (BAM: aligned transcripts to reference genome) to Regtools junctions extract
use Regtools junctions annotate to compare with the reference annotation (GTF)

Exact commands I used for this are:

[..]/regtools/build/regtools junctions extract aln_sorted.bam -o regtools.bed - where the sorted BAM file has index in the same folder
[..]/regtools/build/regtools junctions annotate regtools.bed [..]/ref/genome.fa [..]/ref/annotation.gtf -o regtools_annotate_out

The alignments have been generated with Minimap2.
It did not behave as expected: this gives me as output a file with only junctions annotated as 'new' ($11==N), while other annotations should be possible according to the manual:

DA - This exon-exon junction is present in the transcriptome provided by the user(GTF) The ends of this junction are hence known donor and known acceptor sites according to the GTF file.
NDA - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The ends of this junction are
   known donor and known acceptor sites according to the GTF file.
D - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The donor of this junction is a known donor but the acceptor is novel.
A - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The acceptor of this junction is a known acceptor but the donor is novel.
N - This exon-exon junction is not present(novel) in the transcriptome provided by the user(GTF) The ends of this junction are hence not known donor/acceptor sites according to the GTF file.

My questions are:

am I doing something wrong?
do you have any experience with Regtools and can help debugging this behavior?
do you have other suggestion to achieve the output I am seeking?

You can use the script that comes with minimap2:
# install the k8 javascript shell
curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` $HOME/bin/k8    # assuming $HOME/bin in your $PATH

# acquire the latest minimap2
git clone https://github.com/lh3/minimap2.git

# print novel/erroneous junctions
minimap2/misc/paftools.js junceval -e anno.gtf rna-seq.sam > output.txt

Some example output:
P       read1    4       chr2    227332588       227340301       [(227332588,227340291)]
P       read2    1       chr9    19376379        19376489        [(19376388,19376493)]
P       read2    3       chr9    19378506        19378707        [(19378514,19378707)]
N       read3    1       chr1    160647335       160647472
N       read3    2       chr1    160647581       160659672

Here column 4–6 gives junctions not in the annotation. Column 7, if not empty, gives the list of annotated junctions/introns that overlap the junction in read alignment. Note that a big fraction of outputted records are errors in read alignment.

Hi @user172818, I have two questions on this: (1) could you please explain what the value in $3 represent exactly - and (2) I had cases where the value in $7 is a list of more than one tuple: do you have an explanation ?

@aechchiki 1) the exon number of the read. 2) an intron on a read may overlap multiple distinct introns in annotations due to alt splicing. The list gives all.

For future googlers, Ding et al., 2017 published a review: "Comparison of Alternative Splicing Junction Detection Tools Using RNA-Seq Data".
Table 1 has the list of software in this review.


text:
full visualisation of draft genomes alignment
I have two draft genomes (aprox. 500Mb each), one with ~10'000 scaffolds (genome1.fasta) and one with ~1'000 contigs (genome2.fasta). I would like to visualize their pairwise alignment, which a while ago I found it could easily be done by mummerplot in MUMmer.
In mummer, pairwise draft genomes alignment is supported:

Many times it is necessary to align two genomes that have not yet been completed, or two genomes with multiple chromosomes. This can make things a little more complicated, since a separate alignment would have to be generated for each possible pairing of the sequences. However, both NUCmer and PROmer automate this process and accept multi-FastA inputs, thus simplifying the process of aligning two sets of contigs, scaffolds or chromosomes.

It is still ok if I do it aligning small portions (e.g. contig vs scaffold), as I could get here. However, if I do it all-vs-all, I end up with a "confetti"-plot (self-coined term) instead of a dot-plot :

I am not sure how to interpret this, and what I would like to see is more of a straight line as the one. Probably it's just a matter of filtering, so I'd like to ask for advice.
I used nucmer through delta-filter as follows:
# mummer version //4.0.0beta1

# nucmer
nucmer --prefix=ref_qry genome1.fasta genome2.fasta

# ref_qry.delta generated from previous step
delta-filter -r -q ref_qry.delta > ref_qry.filter

# plot
mummerplot ref_qry.filter -R genome1.fasta -Q genome2.fasta --filter --layout

Reasons for options -r -q in delta-filter step :

To determine a one-to-one mapping of each reference and query, combine the options and use delta-filter -r -q.

I'll also accept other methods for visualisation (doesn't need to be withh MUMmer). Any input is greatly appreciated!

EDIT
I found some solutions which I reported as answers. I will accept anyone giving answer on how to tune MUMmer parameters to achieve the goal I describe in the original question.

http://dgenies.toulouse.inra.fr/

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

Thanks to @user172818.

Also, Mashmap has an embedded visualisation tool generating dot-plots to visualize mappings, ideal for:

highly similar species or to compare different assemblies of the same species.
mapping genome assembly or long reads (PacBio/ONT) to reference genome

An example from their website:


text:
swap genotypes for specific fields of a vcf file
Is there a tool already established to "swap" the ref/alt information in the GT, AD and PL fields of a vcf file?
For the GT field, I need to transform: 0/0 -> 1/1 ;  1/1 -> 0/0 ;  0/1 -> 1/0 ; and finally 1/0 -> 0/1
For the AD field, I need to swap the columns: x,y -> y,x
For the PL field, I need to swap the columns: x,y,z -> z,y,x
Sample input:
Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  0/0:17,0:17:48:0,48,720 0/1:14,8:22:99:167,0,347 .....

Sample output:
    Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  1/1:0,17:17:48:720,48,0 1/0:8,14:22:99:347,0,167 .....

note: the same operations will need to be carried out for all samples if the VCF file has multiple samples (the ... stand for other potential samples).
Before I dive into scripting ( basically playing with regex), I would like to know if something to process this is already available. I looked into bcftools and fixref but they don't seem to help for this.

What do you mean by "note: the ... stand for other columns to be swapped as well,"? You mean carry out the same operations on the GT, AD and PL of any other samples in the VCF file?

yes, better wording here! thanks

How generic does this need to be? I mean, is it enough to give you a solution that works specifically for the order of fields you have (GT:AD:DP:GQ:PL) or do we also need to be able to cope with different orders (so need to infer the position of the AD field from the INFO field)?

for now, a solution for the order of fields I have is okay. I don't have a lot of experience with vcf file, and I thought that the order of fields was "fixed" for all vcfs :9

If we can always assume your fields will be GT:AD:DP:GQ:PL, you can use this relatively simple, brute approach awk script:
awk -vOFS="\t" -F"\t" '{
  ## just print the header lines
  if(/^#/){
    print;
  }
  else{
    ## The number of samples will be the extra number of fields
    ## in addition to the default VCF 9
    sampleNum = NF - 9
    ## Fix each sample
    for(i=1; i<=sampleNum; i++){
      ## get the fields
      num=split($(9+i), fields,":");
      ## genotype
      if(fields[1] ~ "0/0"){ fields[1] = "1/1"}
      else if(fields[1] == "1/1"){ fields[1] = "0/0" }
      else if(fields[1] == "0/1"){ fields[1] = "1/0" }
      else if(fields[1] == "1/0"){ fields[1] = "0/1" }
      ## AD
      split(fields[2],values,",")
      fields[2] = values[2]","values[1]
      ## PL
      split(fields[2],values,",")
      ## Now set these as the sample data
      $(9+i) = fields[1];
      for(k=2;k<=num;k++){
        $(9+i) = $(9+i)":"fields[k]
      }
    }
    print;
  }
}' file.vcf

On your example input, this returns:
Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  1/1:0,17:17:48:0,48,720 1/0:8,14:22:99:167,0,347    .....

thanks @terdon, I guess this means we don't really have a software able to handle it and I'll need to dive into scripting :) your solution works so i'll accept it, although I'd prefer to use some available software in this case.

@aechchiki yeah, I doubt there will be software for this since you are doing something that is very much not standard: you are swapping data around between different alleles.


text:
Get start and end coordinates per chromosome
I will use GATK for SNP calling (HaplotypeCaller). I need to feed the interval file in the command, otherwise I get errors (even though I want to use the whole genome, not a subset  - btw it's not the point of the question, so will not elaborate).
The interval file for GATK can be for example a bed file, with 3 columns: chr_name  chr_start  chr_end. I don't have this file, but have the genome and the reads.
To obtain the intervals per chr/scaf, I proceeded this way:

converted the alignment to interval file, using bam2bed    (BEDOPS)
extract the columns $1,2,3 (using cut)
get the min & max coordinates per chr/scaff with a small R script:

bd <- read.table('all.bed', h=FALSE)
bd_min <- aggregate( bd$V3 ~ bd$V1 , bd, function(x) min(x))
bd_max <- aggregate( bd$V2 ~ bd$V1 , bd, function(x) max(x))
bd_coord <- merge(bd_min, bd_max, by='bd$V1')
write.table(bd_coord, file='coords.bed', sep='\t', quote=FALSE, col.names=FALSE, row.names=FALSE)

These should thus represent the mapped intervals, conceptually same as not feeding the intervals at all (= considering all the sequences).
The ideas for this:

in the bed file: $1 is be the feature, $2 the start position, $3 the end position
so, I need to get the min of $2 for each unique item in $1(for the start coord) and the max of $3 for each item in $1 (for the end coord)

The output looks plausible:
$ cat coords.bed | head -3
Bla_chrm1   678 43860826
Bla_chrm10  181 20381540
Bla_chrm11  343 20367560

My question here:

is this a correct way to proceed?
is there a standardized way to perform this?

My main concern is that I will have to re-build this "coordinate" files for every GATK run, because I think the coordinates per chr/scaf will shift, even slightly, at each GATK run with different datasets.
My other minor concern, is that this approach is a bit slow (given the 160M lines in the current bed file), so an unix-tools solution will also be accepted. I tried to compute myself the max/min using awk (e.g. awk '$3>max[$1]{max[$1]=$3; row[$1]=$0} END{for (i in row) print row[i]}'), but I get different results than with the R approach...

You can try this, since you have the bam file, all of them are aligned to the same reference genome (I hope)(and apologies for the awful awk.. been a while since I used it):
samtools view -H some.bam | awk '{if($1~/@SQ/){print $2"\t1\t"$3}}' | sed 's/[SL]N://g' > coords.bed

For the human genome I got this:
chr1    1   195471971
chr2    1   182113224
chr3    1   160039680
chr4    1   156508116
chr5    1   151834684
chr6    1   149736546
chr7    1   145441459
chr8    1   129401213
chr9    1   124595110

What I did is to get the chromosome length, and make every line start at 1.
If you use R, you can do:
library(Rsamtools)
CL = scanBamHeader(<bam file>)[[1]]$targets
write.table(data.frame(names(CL),start=1,end=CL),
"coords.bed",row.names=FALSE,quote=F,col.names=FALSE,sep="\t")

thanks for your input @StupidWolf. yes I'm using the same reference, to which i also aligned all the reads. it helps the question, on how to get chromosome coordiantes - although not solving the GATK error,  still getting "Input files reference and reads have incompatible contigs: No overlapping contigs found.")

Weird... are you providing the fasta file with the same header? grep ">" <your.fasta> and samtools view -H  | grep "@SQ" ... you should see the same chromosome names

yes, weird, agreed. and yes, i verified that they contain the same fasta IDs, which I also see in the bam header... I opened an issue on GATK forum, as self-comment to my question: https://gatk.broadinstitute.org/hc/en-us/community/posts/360057671431-GATK-HaplotypeCaller-errors-with-or-without-the-L-interval-file-

I upvoted your answer since it's what I was asking for in the original question. I'll deal with the GATK issue on my own. :up:

Thanks dude.. Can you make a small bam file and provide the fasta via ftp or something? Most likely the GATK people might ask for it. I cannot use GATK right now, but I can try with bcftools..

Not sure if I completely understood your question, but if you want simply an interval of your entire genome  you could do instead:
samtools faidx MyGenome.fasta

Obviously that means that you need access to that, not sure if you have.
The resulting fai file contains a format similar to that:
000000F 33203223        94      60      61
000001F 28828106        33756799        60      61
000002F 27810542        63065468        60      61

Where the first column is your chromosome and the second one your length.
I guess you could then just do a awk '{OFS=FS="\t"}{print $1,"0",$2-1}' to get your file.

yes, I do have the genome, and have indexed it. as reported in the question itself, I could obtain (what I understood being) an interval file (just: chrID, start, end). and with @StupidWolf's answer I just tried to expand this interval to start from pos1 for each chr/scaf feed it in GATK. but I still get the same GATK errors (full link in my self-comment to my question) I got for the other interval files. weird, I agree. 100% sure I used one and only ref to map the reads to. I only have that one.

I think the reported error from GATK is clear, the annotation which should look as following:
Bla_scaf1190  830  763
Is wrong as the start is higher then the stop. The question is why you would make an interval file which covers only mapped regions. Isnt the point of intervals to reduce the mapped space ? Otherwise speed should be identical for full space and your intervals

that's one test I tried, and that one error is clear. in the two other options (1.no intervals at all, technically taking the whole genome;  2.intervals for mapped regions of the chr; and I could add what I just tried, so 3.intervals for the whole chr length), I have the other error still persisting "Input files reference and reads have incompatible contigs: No overlapping contigs found"

This probably comes too late, but for the record GATK supports multiple kinds of interval file format. BED is the most common, and there's the Picard format (.interval_list) which is BED-like with a few tweaks. Both require coordinates. However, GATK also supports a simpler format (.list or .intervals) that allows you to just give the contig names (eg chr1) without any coordinates. Sounds like that would satisfy your need. Incidentally, GATK also supports using VCF files as interval file if you want to re-genotype specific regions based on variants located there (ideally paired with the interval padding arguments).
See the GATK doc article for more details:
https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists
(The doc has a few formatting issues but I just notified the team so they'll fix it)


text:
estimate genome size: kmer-based approach from PacBio reads
Can anyone suggest a software/method for kmer analysis using PacBio reads (RSII)?
Something similar to Jellyfish, that I saw in a nice tutorial - but must be suitable for long, noisy reads. kmercounterexact from BBMapTools might also be an option, but again tricky with long, noisy reads.
The closest solution I could find was to estimate the genome size from overlap statistics (for the raw reads, the error corrected preassembled reads or by mapping the raw reads to the assembled contigs) by:

running an initial assembly using a small seed read length
then plot the pre-assembled read overlap histogram.

What do you think of this approach?
I would like to use this kmer analysis to estimate genome size, but for now I could not find much. I do not have companion Illumina reads - if so, I could try the approach proposed by Carvalho et al. 2016
I have reads & in-house genome assemblies (Canu) from different isolates of the same species of fungus to compare.

I don't think there is a method that would estimate a genome size using raw long reads.
The genome size estimates based on raw reads are done by fitting a model to kmer spectra (for instance Genomescope). The kmer spectra built from long reads are really messy due to the high error rate of long reads. That makes fitting of a model quite difficult. These methods assume small error rates. If the error rate is as high as 0.1 (which is normal for long reads), then the probability of having one kmer of 23 bases correct is ~8%. Meaning that if the coverage is 100x, just ~8x of homozygous kmers would be the real genomic kmers and it's very likely that kmer peaks of the sequencing errors, heterozygous kmers and homozygous kmers will just blend into one.
I can imagine, that if your coverage would crazy high (like 600x or more) the kmer spectra could separate errors and genomic kmers. You could also try to decrease kmer size to increase the proportion of correct kmers.
Another option would be to use a trick that they used to generate solid kmers for genome assembly using generalized de brujin graph using noisy reads. Maybe it would be better to take a look on the newer version of the assembler, that is called Flye.
But to be honest, I think for a long reads, a genome assembly is much more reliable way how to predict a genome size than any of the kmer tricks. I would simply assemble the genome, remap back reads, check the uniformity of coverage and conclude the genome size out it.
-- edit --
I tried to find some comprehensive evaluation of the assembly sizes compared to genome sizes, but I could not, so I just checked two examples, one big, one small and they both match quite well. Ammopiptanthus nanus was assembled to 823.74 Mb while kmer estimate was 889 Mb (ref). Plasmodium falciparum was assembled to 23.5 Mb while the original genome size estimate is 22.8 Mb. I know it's just two, but there is also a good reason to think that the genome size will be reasonable, feel free to add here more examples.
If you would be worried about contamination in the assembly, you can first run it though one of the methods to detect a contamination (like Blobology).

ok, I can get genome size from assembly, as discussed. thanks for the input

Using a genome assembly as a way to estimate the genome size is not reliable due to potential contamination causing assembly errors, possible failure to collapse some haplotypes, or assembly errors from inconsistent coverage.

@conchoecia I do agree, but the contamination problem is in kmer spectra analysis as well. But in my experience, long read assemblies I have seen so far were quite close to flow cytometry estimates and consistent with kmer spectra analysis. I will do some more research. Do you have any example where the assembly size was very biased?

@conchoecia better?

Hi @KamilSJaron, the only caveat that I had in mind is that both canu and FALCON do not handle haploid assemblies well if the animal is diploid and highly heteryzygous, like over 5% in the case of FALCON. For example, I am assembling ctenophore genomes which regularly have uncollapsed haplotypes due to their high heterozygosity. As a result, the genome assembly size is approximately twice what is expected based on a kmer spectrum histogram from Illumina reads. Most people probably won't be working with high-heterozygosity animals, but it is important to point out. That's all!

Ah, got you, I ve been there. I was assembling amphioxus genome that IS extremely heterozygous. I ended up setting the parameters so both haplotypes got assembled separately and it worked very well (assembled precisely the diploid genome size)

I have recently come across this submission from Wang et al. (2020), which is a method using the tool kmerfreq on error-corrected long reads (e.g. generated using Canu, FLAS, or LoRMA). It looks to be a promising approach.


text:
calling diploid SNVs from long reads
I'd like to call diploid SNV variants from long-reads data (~80SMRTcells PacBio).
I have generated a draft reference genome for an indivudual from a heterozygous (~4%) species (Canu+Haplomerger2).
I can use this reference for some applications:

RNAseq mapping;
functional annotation;
duplicate detection;
synteny analysis

but not others:

analysis of compound heterozygotes;
allele-specific expression;
variant linkage

This is mainly why I need phasing. FYI, I tried Falcon + FalconUnzip too, but the daligner phase was just taking too much time/resources and I need to find a strategy to make it work in reasonable time (working on it).
Meanwhile, for phasing, It was suggested to use HapCUT2 - from Canu's FAQ

(1) "Avoid collapsing the genome" in Canu to get raw diploid assembly as complete & separate as possible (2) HaploMerger2 to get high-quality reference and alternative haploid assemblies (3) HapCUT2 or other phasing tools to get the high-quality haplotype assembly based on the reference haploid assembly.

However, HapCUT2 requires the input of a diploid SNV file. Which could be done if I had short-reads data, but yet haven't found a method for doing so with in my hands just long-reads data.
Strategies I tried to get the SNVs:

marginPhase : method paper

but I had a issue in running it

Clairvoyante : method paper

but I currently don't have access to suitable GPU machines to build the model - and I can't use the models provided by the authors because my organism is too different from human

Any input would be greatly appreciated!

At 4% heterozygosity, falcon-unzip is probably the best hope. Trio binning would be better but you don't have the right data.

@user172818 is it possible in your knowledge to run Falcon-Unzip on non-Falcon output? From here it looks like it's a no.

You may want to try the longshot tool (https://github.com/pjedge/longshot) developed for calling variants in diploid genomes from long read data.

You could try aligning your reads to the draft reference genome with for example minimap2 and calling variants with freebayes. It appears there is a protocol for long reads: https://github.com/ekg/freebayes#re-genotyping-known-variants-and-calling-long-haplotypes
Disclaimer: I have not run freebayes myself, and it may be that this is a terrible idea. I would expect that calling SNVs (Single-nucleotide variants) from PacBio reads is very challenging.  However, with enough coverage, I think it could be possible.


text:
Are duplicate variants against the VCF standard?
I have been trying to understand an error that the EBI's vcf_validator gives when run on my vcf file. Consider this minimal file:
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  .   AACAC   AATAC   0.00    .   .   GT  0/1
chr12   529516  .   C   T   .   PASS    .   GT  0/1

If I then run vcf_validator on it, it complains about:
$ cat $(LC_ALL="C" ~/vcf_validator_linux  -i test.vcf |& awk '/written to/{print $NF}' )
According to the VCF specification, the input file is not valid
Error: Duplicated variant chr12:529516:C>T found. This occurs 2 time(s), first time in line 6.

This is because the two variants are indeed the same. Changing
AACAC to AATAC at position 529514 just means change C to T at position 529516.
OK, but how is that against the VCF standard? What makes that an invalid VCF file? I am trying to automatically validate input VCF files and stop the analysis if they're invalid, but this seems like a perfectly valid file. Why is  vcf_validator not letting it through? The vcf-validator perl script from vcftools doesn't complain about this. So is it just a bug in  EBI's vcf_validator or is it really somehow illegal to have duplicate variants in a VCF file according to the specs?

What is the use case for having duplicate variant records in a VCF file?

@JohnMarshall none that I know of. But someone sent me a file with such duplicates and I was surprised that they caused the vcf_validator to consider the file invalid, so I asked if there was anything in the specs that would actually make such files invalid.

No. The only point the specification makes is that:

It is permitted to have multiple records with the same POS

I suspect this is intended to let the file list different variants occurring at the same position, but this incidentally allows duplicate variants.
This said, you probably don't want duplicate variants, because it will mess with any statistics you apply over the variants, and will unnecessarily clutter the file and slow down processing, if this happens a lot. You might also get variant annotations on one variant but not another, which can become a problem depending on how you query your variants.

Duplicate variants also point to an issue with upstream analysis. If the variants have different statistics (e.g. in this example one has qual 0.00, and the other has qual "."), then there is a bias in the counting. If the variants have the same statistics, then the lines are redundant.


text:
What is the "primer masterfile" format?
I am working with an amplicon kit from Swift Bioscience and they provide the primer locations in what they call a "masterfile". The format is like this:
$ head -n2 Acel-Amplicon_56g_masterfile.txt
1       43814951        43815064        MPL_S505N-A519T_F/2     43814932        43814951        MPL_S505N-A519T_F/2_primer      43815064  43815081 MPL_S505N-A519T_R/1_primer      GCCTGGATCTCCTTGGTGA     CTTCGGCTCCACCTGGT
1       115256480       115256576       NRAS_G60-S65_new_R/1    115256457       115256480       NRAS_G60-S65_new_R/1_primer     115256576 115256602        NRAS_G60-S65_new_F/2_primer     TGATGGCAAATACACAGAGGAAG CAGGATTCTTACAGAAAACAAGTGGT

This is the format which is understood by Swift's primer clipping tool, but in order to use other tools, I would like to convert it to BEDPE so I can then use standard tools to clip the primers from my bam files. However, I can't do that since I am not sure what the fields in this file are.
This is what I think thy are:

Chromosome
Target region start
Target region end
5' Primer name
5' Primer start
5' Primer end
3' Primer name
3' Primer start
3' Primer end
5' Primer sequence
3' primer sequence

Is my interpretation correct? Also, is this a standard file format or something specific to the Swift tool? I couldn't find any documentation about a "masterfile" format online.

Did you tried to contact with the tech support? Looking at the source code it seems that it accepts some "panel bedfiles"

@Llopis yes, I am talking to their tech support because I have other issues with the file. My question here is about this file format and whether it is some sort of standard format I am not aware of.

this is the swift standard primer file. You are a bit off in your interpretation of the columns (seems like you missed one as there are 12 columns, not 11 which threw everything off):
1) Chromosome
2) Target Start
3) Target Stop
4) Target Name (change which affected everything downstream in counting)
5) 5' Primer Start
6) 5' Primer Stop
7) 5' Primer Name
8) 3' Primer Start
9) 3' Primer Stop
10) 3' Primer Name
11) 5' Primer Sequence
12) 3' Primer Sequence
To convert to a BEDPE describing the primer locations, you only will need columns 1, 5, 6, 8, and 9 (and optionally 4). The last letter of the Region name for the denotes the strand orientation of the region. For example Region 1 has _F denoting forward or (+) strand for BEDPE usage, while region 2 has _R in name denoting reverse or (-) strand for BEDPE. Putting all this together, you can convert to BEDPE with:
$ awk -v OFS="\t" '{
            o=".";
            if($4~/_F\//){
                o="+"
            }
            else if($4~/_R\//){
                o="-"
            }
            print $1,$5,$6,$1,$8,$9,$4,".",o}' foo.masterfile
1   43814932    43814951    1   43815064    43815081    MPL_S505N-A519T_F/2 .   +
1   115256457   115256480   1   115256576   115256602   NRAS_G60-S65_new_R/1    .   -

Thanks, I meant 4 was the target name, I just wrote it wrong in the question. But why would I use the first 4 columns? I want the positions of the primers, not of the target. So, for a bedpe file, wouldn't I need columns 1,5,6,1,8,9 and then - or + depending on the "F" or "R"? Aren't the 1st 4 columns describing the target and not the primer?

I see, it was not clear to me that you wanted the location of primers and not the target. The primer orientation is also included in the primer names. The columns still remain the same, so if you would like primer info just use the information in the corresponding primer columns. The orientation of each primer is in its name as well. You would use columns 1, 5,6,7 for 5' primer and 1, 8,9,10 for 3' primer. Using name columns only for primer orientation (but BEDPE can also accept a 4th column for the name of the region if you'd like).

Yes, sorry, my question wasn't very clear on that (I edited now). Do you have any idea what the numbers at the end of the target names are? I mean the 2 at the end of MPL_S505N-A519T_F/2 and the 1 at the end of NRAS_G60-S65_new_R/1. Could they be the phase? And, if so, with respect to what?

By the way, I took the liberty of adding an awk command to make the conversion since I think that might be useful to future visitors. I hope you don't mind.

The awk command adds a nice touch! I am not super sure what the 1 and 2 mean with respect to the genes or orientations as with IDT orders I know they have a different meaning. Specifically, I think they refer to exon number but that could not be the case here. I don’t want to give you misinformation, but that is something you should check first with respect to your target set.


text:
cn.MOPS fails with 'missing value where TRUE/FALSE needed'
I am trying to use cn.MOPS to call CNVs on a set of whole exome sequencing data bam files. My script (up to the point where it fails) is:
source("https://bioconductor.org/biocLite.R")
library(cn.mops)
BAMFiles <- c('sample1.bam', 'sample2.bam', 'sample3.bam', 'sample4.bam', 'sample5.bam', 'sample6.bam', 'sample7.bam')
segments <- read.table("region.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)
resCNMOPS <- exomecn.mops(X)

The last command, resCNMOPS <- exomecn.mops(X) fails with:
Error in if (all(segMedianT == 0)) { :
  missing value where TRUE/FALSE needed
In addition: Warning message:
In normalizeChromosomes(X, chr = chr, normType = normType, qu = normQu,  :
  Normalization for reference sequence  chrUn_gl000228 not applicable, because of low number of segments

If I run traceback immediately afterwards, I get:
> traceback()
4: FUN(newX[, i], ...)
3: apply(sINI[chrIdx, , drop = FALSE], 2, segment, minSeg = minWidth,
       ...)
2: cn.mops(input = input, I = I, classes = classes, priorImpact = priorImpact,
       cyc = cyc, parallel = parallel, normType = normType, normQu = normQu,
       norm = norm, sizeFactor = sizeFactor, quSizeFactor = quSizeFactor,
       lowerThreshold = lowerThreshold, upperThreshold = upperThreshold,
       minWidth = minWidth, segAlgorithm = segAlgorithm, minReadCount = minReadCount,
       useMedian = useMedian, returnPosterior = returnPosterior,
       ...)
1: exomecn.mops(X)

I am using R version 3.4.4 on an Ubuntu 16.04.5 LTS machine.
My bed file only has one region which is on chrUn_gl000228:
$ grep chrUn_gl000228 region.bed
chrUn_gl000228  112604  113879

If I remove that region and rerun the same code, it works as expected. So, my question is, how can I detect this and fix it automatically? Given a collection of bam files and the bed file describing the capture regions used to sequence the samples and generate the bam files, how can I check if cn.MOPS will complain about any of the segments?
Presumably, once I have built these two objects:
segments <- read.table("region.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)

I should be able to then check them and see if they have the right data and, if not, remove the offending segment from the GRanges object. But how can I do this?

You start your question with biocLite, which indicates that this is not the latest version of R and Bioconductor. Which versions are you using? Have you tried with the most up to date version? Why does it complain that you have low number of segments? Could you post how many segments do you have on X?

@llrs no, I haven't used the latest, but upgrading is not an option (this is a large server). I would be happy to post that, but can you give me a command for it? My R knowledge is extremely limited. That sort of command is exactly the kind of thing I am asking about. I have no idea how to inspect and get statistics about the contents of a GRanges object.

Also, hang on, why are you saying it isn't the latest version? I got that code from the documentation of cn.MOPS. Is there a better way to load biocLite?

Well you will need to use sessionInfo() and paste here the result to see the version of R and Bioconductor. The most up to date documentation of cn.MOPS is on the bioconductor page AFAIK (if packages are not up to the standards they are kicked out of the project). Try with summary(X) and dim(X) or simply print the object and post the first lines

I said you weren't using the latest version because in the latest version the recommended way to install packages changed to BiocManager::install("cn.Mops") you need a package in CRAN to install packages of Bioconductor

@llrs yes, I added the R version (3.4.4). I had already tried summary(X) but that only prints "GRanges object with 191743 ranges and 7 metadata columns".  dim(X) gives NULL.

Let us continue this discussion in chat.

I managed to track down what was producing this error. The problem was that my bed file only had one segment for chromosome chrUn_gl000228. However, the cn.MOPS method relies on normalizing read counts across multiple segments for each chromosome. Since this chromosome only had a single segment, the exomecn.mops(X) was failing.
So, either remove the single segment from the bed file, or use a bed file that has at least 2 segments on each chromosome. I created a fake segment on chrUn_gl000228 and confirmed that the program worked as expected and didn't throw any errors.


text:
Why do SAM and BAM use different coordinate systems?
BAM files are, at least as far as I know, simply binary compressed versions of SAM files. They have the exact same information and are used in the same way. Why then does the SAM format use a 1-based coordinate system while BAM uses a 0-based one?
This is the relevant section of the SAM format specification (emphasis mine):

1-based coordinate system A coordinate system where the first base of
a sequence is one. In this coordinate system, a region is specified by
a closed interval. For example, the region between the 3rd and the 7th
bases inclusive is [3, 7]. The SAM, VCF, GFF and Wiggle formats are
using the 1-based coordinate system.

0-based coordinate system A
coordinate system where the first base of a sequence is zero. In this
coordinate system, a region is specified by a half-closed-half-open
interval. For example, the region between the 3rd and the 7th bases
inclusive is [2, 7). The BAM, BCFv2, BED, and PSL formats are using
the 0-based coordinate system.

The two coordinate systems are a very common source of error in bioinformatics, so there must be some reason why these two largely equivalent formats use different systems. For that matter, why do the VCF and BCF formats, whose relationship is equivalent to that between SAM and BAM, also do this? I am guessing there is a good reason and it is equally applicable to both pairs of sister formats.
Note that the samtools manual states that "all position coordinates are 1-based" when giving regions, so it must handle the conversion internally. Indeed, I can confirm that this command returns identical md5sums for both inputs:
md5sum <(samtools view -b sample.bam chr1:36931633-36931871 | samtools view) \
       <(samtools view sample.bam chr1:36931633-36931871)

Arithmetic on a zero-based coordinate system is less complicated than that on a one-based system, so it appears zero-based is often (not exclusively) used for binary data formats, like BAM or bigBed, or text formats like BED, where computers are used to more efficiently calculate lengths of or set operations on intervals.
A more complete answer on the benefits of zero-based accounting is available here, which references the original note (EWD831) on this subject from Dijkstra.
When reading other textual data formats like SAM or VCF, a one-based index is perhaps more "natural" for humans used to counting digits on their hands. If you visually inspect an annotation on the closed interval [3,7] on a line that denotes some abstract chromosome and starts at position 1, it is clear to you where those data points are positioned. Colloquially, you understand that the interval starts on the third position and ends on the seventh. Or, perhaps, people would typically label the initial item in this set the first item, and not the zeroth item.
A one-based index therefore seems to get used often with formats used for visual inspection or presentation of genomic intervals or sequences, such as in genome browsers — or via CLI tools like samtools, when exporting reads to human-readable SAM.
The language used for working with and organizing data can determine choice of indexing, as well. For programs written in C and offshoot languages, including an API like htslib, binary data is typically worked with as chunks of memory accessed via pointers or arrays, and it is convention to iterate over that data in increments that use a zero-based index. Contrariwise, R's roots in FORTRAN and focus on linear algebraic operations (and matrix indexing) may perhaps explain why that language uses one-based indices for its lists, matrices, and vectors.
The designer of the BAM and SAM formats, Heng Li, comments below(*):

I think 1-based is easier for biologists and
mathematicians while 0-based is easier for C/C++/Java programmers.
SAM/VCF is for biologists to read. BAM/BCF is for a small number of
programmers. I wrote the initial specs of sam, bam and bcf2. PS: when
sam and vcf are stored in memory by htslib, everything is 0-based in
accordance with C conventions. bam and bcf are largely direct dumps of
their in-memory representations.

Note: Original comments were unfortunately deleted from answer by someone, but the quotes were taken verbatim at the time they were written.

That's a fair point about human readable formats, but calling SAM human readable is a bit of a stretch, really. Yes, it's text, but it is designed to be consumed by machines, not humans.