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

For example, from my BLAST search:

except I would like the graph to contain information including information for the distribution of all 6139 hits across my query sequence.

Since this is a large number of sequences to visualize a distribution, I am considering using a score of the number of hits on a certain region, and having an output more similar to:

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

Thanks.

# Answer
You probably want to include query start (qstart) and query end (qend) in your blast output.

Something like this:

```
blastn -query your.fasta -out blast.out.txt -db your.db -outfmt '6 qseqid sseqid qstart qend length evalue'

```

In R you can take the "qstart:qend" from each line for density plot.

There are many ways in R to plot the densities of these start and end amino acids.

Let me show an example with a small data frame:

```
qstart <- c(200, 300, 250, 400, 500)
qend <- c(300, 450, 400, 600, 650)

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

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

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

```
> 4 votes

---
Tags: r, blast, visualization
---


title: Formula for k-mer coverage

text:
# Question
Title: Formula for k-mer coverage
Let $C$ be base coverage, $R$ is the length of reads and $K$ is the length of $k$-mer. Then $k$-mer coverage $C\_k$ can be computed as $C\_k = C\cdot(R - K + 1)/R$.

Could someone please explain why is this equation valid (I'm mostly confused as why it is divided by $R$)?

Source: Velvet manual Section 5.1

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

# Answer
> 5 votes
I was still puzzled from the answers, so I tried to calculate with all the steps. I take this definition "$C\_k$ is the number of reads containing a k-mer." and corresponding definition for coverage ($C$): "$C$ is the number of reads covering a base".

Coverage is $C = \frac{T \cdot R}{L}$, where $T$ is total number of reads, $R$ is read length and $L$ is length of genome. Given the $C\_k$ definition, $C\_k = \frac{T (R - K + 1)}{L-K+1}$, where $R - K + 1$ is just number of kmers in a read, and $L-K+1$ is number of kmers in a genome. Then,

$$C\_k = \frac{T (R - K + 1)}{L-K+1} = \frac{T (R - K + 1)}{L-K+1} \cdot \frac{R}{R} = \frac{R - K + 1}{R} \cdot \frac{T \cdot R}{L - K + 1}$$

since $L \>\> K$, we can approximate $L - K + 1 \approx L$, then we reduce the expression to

$$\frac{R - K + 1}{R} \cdot \frac{T \cdot R}{L} = \frac{R - K + 1}{R} \cdot C$$

which is the formula for $C\_k$.

---
Tags: k-mer
---


# Question
Title: What database do gene IDs starting with sp_v2 refer to?
I am a student trying to analyze GEO2r datas for one of my courses. The IDs given in the output are different for different series. I need to convert all of them to a similar format.

In this process I encountered the following type of ID which I don't know its origin:

```
”ID"    "logFC"
"SP_v2 4634"    "-0.9897758"
"SP_v2 3382"    "-0.8391782"
"SP_v2 4210"    "-1.1693583"
"SP_v2 2117"    "-1.0504727"
"SP_v2 3488"    "-0.9756444"
"SP_v2 1128"    "-0.8289103"
"SP_v2 2735"    "-0.8629999"
...

```

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

# Answer
This is possibly a wild goose chase, but a lot of searching led me to this sample:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM2536852

Which appears to exhibit accessions of the right format. This in turn leads to the platform:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL22166

Which is a custom spotted cDNA array, which the metadata seems to suggest is a human platform.

The full annotation table gives gene symbols and Entrez IDs for most of the probes on the array:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc=GPL22166&id=53508&db=GeoDb\_blob144

To take the example of the first line of your results:

```
ID          ORF     Entrez gene SEQUENCE
SP_v2 4634  ETNK1   55500       AAAGCAGCTTCATCTTTCAAAATTGATTTGCTCTGGTTTT

```

The Entrez Gene record for that ID matches up:

https://www.ncbi.nlm.nih.gov/gene/?term=55500
> 5 votes

---
Tags: identifiers
---


# Question
Title: Documentation for Transfrag class codes (cuffcompare)
Is there any extensive documentation or description for each class of Transfrag class codes as reported by the cuffcompare tool in the Cufflinks package?

Official doc might not be the best.

E.g. what does *contained* mean (class code: c)? Or, what is a *generic exonic overlap* (class code: o)?

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

---
Tags: transcriptome
---


# Question
Title: How GFF3 attributes (9th column) varies from one gene prediction algorithm to another
GFF3 files are in tabular format with 9 fields per line, separated by tabs. The first 8 fields share almost same data structure, but the 9th field varies a lot depending on feature type and gene prediction algorithm.

Presently I am trying to build my own GFF module that will parse a GFF file and return one structured data object (a **python dictionary structure**). For that, I need a good understanding of the 9th field of the GFF format.

So can any one please help me by providing various GFF files?

# Answer
The first place to start is the GFF3 specification. This is the official word on what is and is not allowed in a GFF3 file. For example, users can define arbitrary attribute keys, so long as they do not begin with an uppercase letter (these are reserved for "official" use).

But your question doesn't seem to be about what is *allowed*, but what is *commonly used*. I have a question for you: which gene predictors are YOU using? Or are you using gene annotations produced by others?

Here are a few examples.

## NCBI RefSeq

If you're using GFF3 from NCBI, chances are it looks something like this.

```
NC_007070.3     Gnomon  gene    329235  331223  .       +       .       ID=gene14;Dbxref=BEEBASE:GB42168,GeneID:551678;Name=LOC551678;gbkey=Gene;gene=LOC551678;gene_biotype=protein_coding
NC_007070.3     Gnomon  mRNA    329235  331223  .       +       .       ID=rna27;Parent=gene14;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;Name=XM_624067.4;gbkey=mRNA;gene=LOC551678;model_evidence=Supporting evidence includes similarity to: 66 ESTs%2C 24 Proteins%2C and 99%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 113 samples with support for all annotated introns;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    329235  329459  .       +       .       ID=id117;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    329850  330082  .       +       .       ID=id118;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    330166  330301  .       +       .       ID=id119;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    330376  331223  .       +       .       ID=id120;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  CDS     329333  329459  .       +       0       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     329850  330082  .       +       2       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     330166  330301  .       +       0       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     330376  330416  .       +       2       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1

```

## MAKER

The MAKER annotation workflow (paper, software) is a pretty commonly used gene annotation tool, and produces GFF3 output like this.

```
scaffold_12     maker   gene    652527  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959;Name=maker-scaffold_12-augustus-gene-0.959
scaffold_12     maker   mRNA    652527  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1;Parent=maker-scaffold_12-augustus-gene-0.959;Name=maker-scaffold_12-augustus-gene-0.959-mRNA-1;_AED=0.24;_eAED=0.18;_QI=0|0|0|0.66|0.5|0.33|3|0|218
scaffold_12     maker   exon    652527  652817  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1203;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   exon    654877  655170  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1204;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   exon    655272  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1205;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     652527  652817  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     654877  655170  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     655272  655343  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   gene    941547  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851;Name=snap-scaffold_12-processed-gene-0.851
scaffold_12     maker   mRNA    941547  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1;Parent=snap-scaffold_12-processed-gene-0.851;Name=snap-scaffold_12-processed-gene-0.851-mRNA-1;_AED=0.22;_eAED=0.22;_QI=0|0|0.25|0.25|1|1|4|661|95
scaffold_12     maker   exon    941547  941631  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1206;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    942343  942367  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1207;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    942780  942920  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1208;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    943200  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1209;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     941547  941631  .       +       0       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     942343  942367  .       +       2       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     942780  942920  .       +       1       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     943200  943236  .       +       1       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   three_prime_UTR 943237  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:three_prime_utr;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1

```

## GeneMark, SNAP, Augustus

If we look at the output of several common *ab initio* gene prediction tools, none of them actually produces GFF3 by default. Here is some output from GeneMark (paper, software).

```
543   1   +  Initial         799051     799097         47          1 2 - -
543   2   +  Terminal        799236     799266         31          3 3 - -

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

```

Here is some output from SNAP (paper, software).

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

```

Both of these formats appear to be tab-delimited, but that's where the similarity with GFF3 ends.

The output of Augustus (paper, software) looks like GTF (a variant of the GFF format) once you lose all the lines beginning with a `#` symbol, but Augustus also has a GFF3 output mode.

```
# Predicted genes for sequence number 1 on both strands
# start gene g1
scaffold_12     AUGUSTUS        gene    2841    3400    1       +       .       g1
scaffold_12     AUGUSTUS        transcript      2841    3400    1       +       .       g1.t1
scaffold_12     AUGUSTUS        start_codon     2841    2843    .       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        intron  3027    3097    1       +       .       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        CDS     2841    3026    1       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        CDS     3098    3400    1       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        stop_codon      3398    3400    .       +       0       transcript_id "g1.t1"; gene_id "g1";
# protein sequence = [MAIKNAEHDLRVIVDAIEGLGLKVAPHKTEAMAFPASALCGRRGAAPPKIRLGGSSILVGSRSRWYISGHSENSSKSP
# RTEGKETTPLQQRDPLDAPLWVSGVVAHCCGGPEGQEGCPGLAAQGSDQGVLRIRDGLLCGYDGCGDHRPRPSDSSAGGGLCRP]
# end gene g1
###
# start gene g2
scaffold_12     AUGUSTUS        gene    4712    15229   0.21    +       .       g2
scaffold_12     AUGUSTUS        transcript      4712    15229   0.21    +       .       g2.t1
scaffold_12     AUGUSTUS        start_codon     4712    4714    .       +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        intron  4858    5591    0.74    +       .       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        intron  5686    15028   0.23    +       .       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     4712    4857    0.72    +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     5592    5685    0.42    +       1       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     15029   15229   0.43    +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        stop_codon      15227   15229   .       +       0       transcript_id "g2.t1"; gene_id "g2";
# protein sequence = [MGRNSHRSCCVVNCKITSAKSDCKFYKFPTAKWKINQRKMWVAAVKRQKYIKDEISHAETQTEITEVTGATKVNYANK
# KYICLLFVRTYVLRMLVDVALSNLRFSLFGIRKSLEIFGQSEKADQTRWRLPSCEMEWIESRKGKMRE]
# end gene g2
###

```

## Parsing attributes

So...how should YOU handle attributes in GFF3's 9th column? That depends a lot on what you want to do with the data.

The most important attributes are `ID` and `Parent`, which are used to define relationships between features and subfeatures. (These relationships implicitly define a directed acyclic graph of features, although most GFF3 parsers don't directly support traversal of this graph.)

But not only do you have to handle the pre-defined attributes discussed in the GFF3 specification, you also have to be able to handle any number of arbitrary attributes whose keys you may not know beforehand.

The simplest way to handle this would be to parse the attribute column into a dictionary of key/value pairs. Once it's in this form, it's trivial to see what attributes are there and how to access them.

For example, if we have a feature with the following attributes...

```
ID=mRNA42;Parent=gene19;integrity=0.95;foo=bar

```

...we would want it in a dictionary like so.

```
attributes = {
    'ID': 'mrna42',
    'Parent': 'gene19',
    'integrity': '0.95',
    'foo': 'bar',
}

```

Python code to parse that might look something like this.

```
attributes = dict()
for keyvaluepair in attributestring.split(';'):
    for key, value in keyvaluepair.split('='):
        attributes[key] = value

```

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

---
Tags: file-formats, python, gff3, sequence-annotation
---


# Question
Title: Bioperl - how can i print first result of search sequence per iteration?
```
use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;

$query = "LEGK";
$query_obj = Bio::DB::Query::GenBank->new(-db => 'protein',
                                          -query => $query );

$gb_obj = Bio::DB::GenBank->new;

$stream_obj = $gb_obj->get_Stream_by_query($query_obj);

while ($seq_obj = $stream_obj->next_seq) {
    # do something with the sequence object
    print ">$query",' ', $seq_obj->display_id, ' ', $seq_obj->desc,"\n", $seq_obj->seq[,'\n';

```

How can I print first occurrence of protein sequence?

# Answer
> 4 votes
In the simplest case, if you just want to stop after the first record was printed, you can just add `exit` (I also corrected the syntax errors you had and added `use strict` and `use warnings`; I suggest you get into the habit of using those two, they save you from a lot of grief in the long run):

```
#!/usr/bin/env perl
use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;
use strict;
use warnings;
my $query = "LEGK";
my $query_obj = Bio::DB::Query::GenBank->new(-db => 'protein',
                                          -query => $query );
my $gb_obj = Bio::DB::GenBank->new;
my $stream_obj = $gb_obj->get_Stream_by_query($query_obj);
while (my $seq_obj = $stream_obj->next_seq) {
 print ">$query",' ', $seq_obj->display_id, ' ', $seq_obj->desc,"\n", $seq_obj->seq, "\n";
    exit;
}

```

If instead, you want to only process the first result for a list of multiple queries:

```
#!/usr/bin/env perl
use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;
use strict;
use warnings;
my @queries = qw(LEGK TP53 ACT1);
query:foreach my $query (@queries) {
 my $query_obj = Bio::DB::Query::GenBank->new(-db => 'protein',
                                          -query => $query );
 my $gb_obj = Bio::DB::GenBank->new;
 my $stream_obj = $gb_obj->get_Stream_by_query($query_obj);
 while (my $seq_obj = $stream_obj->next_seq) {
  print ">$query",' ', $seq_obj->display_id, ' ', $seq_obj->desc,"\n", $seq_obj->seq, "\n";
        ## Move to the next query
        next query;
    }
}

```

Or, to print only the first sequence if multiple quries have identical sequences:

```
#!/usr/bin/env perl
use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;
use strict;
use warnings;
my @queries = qw(LEGK TP53 ACT1);
foreach my $query (@queries) {
 my %seen;
 my $query_obj = Bio::DB::Query::GenBank->new(-db => 'protein',
                                          -query => $query );
 my $gb_obj = Bio::DB::GenBank->new;
 my $stream_obj = $gb_obj->get_Stream_by_query($query_obj);
 while (my $seq_obj = $stream_obj->next_seq) {
  if ($seen{$seq_obj->seq}) {
   ## Here, since we're not using an argument (next query), this will move
   ## to the next iteration of the current loop and not the `foreach`.
   next;
  }
  print ">$query",' ', $seq_obj->display_id, ' ', $seq_obj->desc,"\n", $seq_obj->seq, "\n";
    }
}

```

---
Tags: database, public-databases, perl
---


title: Chip-Seq data description

text:
# Question
Title: Chip-Seq data description
I have the following Chip-Seq data and could not found a description of it. Can you help me with it? More, I would like to find out the order of the nucleotides in binding regions; that is, what number they are.

```
@GWZHISEQ01:319:HA1NPADXX:1:1101:1182:2196 1:N:0:ACTGAT
GGAGGAAGTGTGTTGCTGCAGACAGTAGCAACCAGACCCACACTGCGCGTA
+
CCCFFFFFHHFHHIJJJJIJJIIIIHGJJJJJJJCGHIJJIJJFHIIJJII
@GWZHISEQ01:319:HA1NPADXX:1:1101:1272:2109 1:N:0:ACTGAT
NCCACGCTAGGCTCAGCTTGTCGGCCTGGCTAAGCAGTTGCGAAAGTGCGC
+
#1=DDDFFHHHHHJJIJJJJJJIGEHGGJIGHJFHGJIJJJJJJJ@CHCBG
@GWZHISEQ01:319:HA1NPADXX:1:1101:1418:2141 1:N:0:ACTGAT
GCACGCACTACCCAGAGATCATCCAAAGCCTGAAGCCACAGGGCGCACTCG
+
CCCFFFFFGHGHHJJJIIJJJJIIHGGGIJIIJJIIIJJJIJJJIJIJIJH
@GWZHISEQ01:319:HA1NPADXX:1:1101:1302:2160 1:N:0:ACTGAT
ACACTCTTTCCCTACACGACGCTCTTCCGAGATCGGAAGAGCACACGTCTG
+
CCCFFFFFHDHHHJJJJIJJJJIIIIIJIHJIJJJIIJHIIJJJJJJJJHH

```

# Answer
> 4 votes
That's a fastq file, you will want to align it to the genome, call peaks, and then use something like MEME to determine binding motifs.

---
Tags: chip-seq
---


# Question
Title: RNAseq: Z score, Intensity, and Resources
I'm very new to bioinformatics in general, and I'm trying to understand some basic concepts.

I have RNAseq data, and bioinformatics people tell me that intensities cannot be compared across patients. So there are all of these pipelines to compare intensities to Z scores--are those as simple as just plugging data into a bioconductor package?

It would be great to get some overview description of why Z values are important/why you can't compare intensities across patients, and/or some pointers towards any resources that I can read to learn more.

Thanks much for your time!

# Answer
> 3 votes
It depends on what test or analysis you want to do, whether you need intensities (expression values) or z-scores.

If you want to do statistical analysis, such as finding differentially expressed genes between groups of patients (e.g., with limma), you don't want to use z-scores. But you use normalized intensities (for microarrays) or start with raw counts for RNAseq. Then follow the user guide (limma has a great user guide, also for starters).

For visualization in heatmaps or for other clustering (e.g., k-means, fuzzy) it is useful to use z-scores. Z-scores are a form of transformation (scaling), where every genes is sort of "reset" to the mean of all samples, using also the standard deviation. If you want to know exactly what a z-score is, a simple google search can tell you the details.

In R you can use the scale function for z-score transformation. Be aware that the function works on columns though. Which means you have to transpose your matrix first if genes are in the rows, and then transpose them back after transformation.

# Answer
> 2 votes
In cancer genomics experiments often involve sequencing a large number of tumour samples, but few or no matched normals. This is partly due to financial constraints, but also due to ethical ones: it is easy to get ethical approval for tumour tissue, which is removed as a matter of course in a patient's care anyway, whereas it is hard to get approval for normal tissue, which would not normally be removed and would require an additional procedure.

However, this presents a problem. Most of the procedures we use in standard gene expression analysis involve comparing one group to another after suitable normalisations of the data. But here we have no suitable comparison group. To get around this problem in cancer RNA-seq analysis we do something called outlier detection. That is, we assume that for a given gene, the majority of the patients have "normal" expression, and look for the few patients that don't look like the rest. We do this with the Z-score.

To calculate the Z-score for an observation we subtract the mean of all observations and divide by the standard deviation. Thus the Z score of an observation is how many standard deviations an observation is from the mean of all observations - or how unusual it is.

Thus for Gene A in Patient 1, we calculate how many standard deviations Gene A is from the mean of Gene A across all patients. A very big (or small/negative) value tells us that the expression of Gene A is unusual in patient 1 compared to the other patients.

(NB: in some cancer studies they use the median and median absolute deviation rather than the mean and standard deviation to compute a "robust Z-score" as the data is unlikely to be normally distributed).

---
Tags: rna-seq, bioconductor, python, normalization
---


# Question
Title: Normalization methods with RNA-Seq ERCC spike in?
ERCC spike-in is a set of synthetic controls developed for RNA-Seq. I'm interested in using it to normalize my RNA-Seq samples. In particular, I'd like to use the spike-ins to remove technical bias and any variation that should not be part of my analysis.

The site doesn't give any details on how I can do that.

**Q:** What are the possible normalization strategies? Can you briefly describe them?

# Answer
> 9 votes
You may consider using RUVSeq. Here is an excerpt from the 2013 Nature Biotechnology publication:

> We evaluate the performance of the External RNA Control Consortium (ERCC) spike-in controls and investigate the possibility of using them directly for normalization. We show that the spike-ins are not reliable enough to be used in standard global-scaling or regression-based normalization procedures. We propose a normalization strategy, called remove unwanted variation (RUV), that adjusts for nuisance technical effects by performing factor analysis on suitable sets of control genes (e.g., ERCC spike-ins) or samples (e.g., replicate libraries).

RUVSeq essentially fits a generalized linear model (GLM) to the expression data, where your expression matrix $Y$ is a $m$ by $n$ matrix, where $m$ is the number of samples and $n$ the number of genes. The model boils down to

$Y = X\*\beta + Z\*\gamma + W\*\alpha + \epsilon$

where $X$ describes the conditions of interest (e.g., treatment vs. control), $Z$ describes observed covariates (e.g., gender) and $W$ describes unobserved covariates (e.g., batch, temperature, lab). $\beta$, $\gamma$ and $\alpha$ are parameter matrices which record the contribution of $X$, $Z$ and $W$, and $\epsilon$ is random noise. For subset of carefully selected genes (e.g., ERCC spike-ins, housekeeping genes, or technical replicates) we can assume that $X$ and $Z$ are zero, and find $W$ - the "unwanted variation" in your sample.

# Answer
> 4 votes
We have added ERCC spike-ins to all our RNASeq data, just in case other people might find it useful in the future. However, I have never used it in my own analyses because I can't think of a reasonable way that it could be used.

The typical recommendation for ERCC is to add it in proportion to the input RNA amount, but that makes an assumption that total cell RNA counts are similar across different cells (which is demonstrably false by looking at single cell RNASeq results).

I have yet to think of a situation where ERCC would provide better results than a "housekeeping" gene set sampled from the original reads.

---
Tags: rna-seq, transcriptome, normalization
---


# Question
Title: Extract nanopore read ID & start times from fastq file
I have a fastq file from minION (albacore) that contains information on the read ID and the start time of the read. I want to extract these two bits of information into a single csv file.

I've been trying to figure out a grep/awk/sed solution, but without success.

E.g.

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

```

Should produce

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

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

```

The `awk` command gets the first of every 4 lines, printing the first and fifth "word". `sed` is then used to strip the initial `@` and replace `start_time=` with `,`. The output on your example file is:

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

```
> 9 votes

# Answer
Since the string `start_time` will only appear on the header line, or else you don't have a valid fastq file, you can simply do:

```
$ perl -ne '/^@(\S+).*start_time=(.*)/ && print "$1, $2\n"' file.fastq
93a12f52-95e5-40c7-8c3e-70bf94ed0720,2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c,2017-07-04T06:56:41Z

```

Alternatively, since you mentioned `awk` and `sed`:

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

```

Or, doing the whole thing in `awk`:

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

```

And if the `@` annoy you:

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

```

And in `sed`:

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

```

Or, if your sed supports it:

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

```

Finally, since `grep` can't do replacements, to use it you would have to do something like:

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

```

And to get the commas:

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

```
> 6 votes

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

* filename
* **read\_id**
* run\_id
* channel
* **start\_time**
* duration
* num\_events
* template\_start
* num\_events\_template
* template\_duration
* num\_called\_template
* sequence\_length\_template
* mean\_qscore\_template
* strand\_score\_template
* complement\_start
* num\_events\_complement
* complement\_duration
* num\_called\_complement
* sequence\_length\_complement
* mean\_qscore\_complement
* strand\_score\_complement
* sequence\_length\_2d
* mean\_qscore\_2d
> 3 votes

---
Tags: fastq, nanopore, awk
---


title: Multiple sequence alignment with R

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

Any recommendations and feedback are really appreciated :)

# Answer
The ‹msa› package on Bioconductor does exactly that.

It doesn’t hand the results on a silver platter, though: you’ll need to read the vignette carefully to learn how to use it. But after that it’s pretty powerful.
> 4 votes

---
Tags: r, sequence-alignment
---


# Question
Title: Biopython: resseq doesn't match pdb file
I have a PDB file, and I need to extract its residue sequence numbers (resSeq's). Based on manual inspection of the first few lines of the PDB file (pasted below), I would think that resSeq's should begin with 22, 23. However, Biopython's PDB module suggests otherwise (output attached below as well). I wonder if it's a Biopython bug or if I have problems understanding the PDB format.

```
ATOM      1  N   GLY A  22      78.171  89.858  59.231  1.00 21.24           N
ATOM      2  CA  GLY A  22      79.174  88.827  58.999  1.00 20.87           C
ATOM      3  C   GLY A  22      80.438  89.415  58.391  1.00 21.89           C
ATOM      4  O   GLY A  22      80.362  90.202  57.440  1.00 23.18           O
ATOM      5  N   LEU A  23      81.588  89.069  58.972  1.00 21.51           N
ATOM      6  CA  LEU A  23      82.895  89.555  58.527  1.00 20.80           C
ATOM      7  C   LEU A  23      83.288  89.020  57.162  1.00 22.41           C
ATOM      8  O   LEU A  23      82.889  87.923  56.788  1.00 22.93           O
ATOM      9  CB  LEU A  23      83.973  89.232  59.560  1.00 20.97           C
ATOM     10  CG  LEU A  23      84.225  87.818  60.062  1.00 13.32           C
ATOM     11  CD1 LEU A  23      85.448  87.888  60.939  1.00 15.24           C
ATOM     12  CD2 LEU A  23      83.035  87.258  60.829  1.00 12.21           C

```

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

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

```

# Answer
Let's check it. I copied your file fragment into `a.pdb`.

```
>>> from Bio import PDB
>>> for res in PDB.PDBParser().get_structure('pdb', 'a.pdb')[0].get_residues():
...   print(res.get_full_id())
...
('pdb', 0, 'A', (' ', 22, ' '))
('pdb', 0, 'A', (' ', 23, ' '))

```

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

---
Tags: pdb, python
---


# Question
Title: Best way to detect long insertions in bisulfite sequencing data?
I am interested in identifying indels in whole genome bisulfite sequencing data (76bp paired end). Currently, I do this by setting the `-rfg` and `-rdg` affine gap penalty scores for bowtie2 to more permissive values than the default 5+3N and mapping using the bisulfite sequencing alignment wrapper, bismark.

My question is, what values of `-rfg` and `-rdg` will allow me to identify longest possible indels without sacrificing alignment quality? Is it better to set the affine penalty to zero with a high penalty for initially opening an indel (ex. 8+0N)? Or is it better to keep the initial penalty low and having a nonzero penalty for extension (ex. 2+1N)?

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

---
Tags: ngs, sequence-alignment, read-mapping, methylation, indel
---


# Question
Title: FASTQC overrepresented sequences after trimming
I have a set of RNA-seq samples from different experiments (Single and Paired End, depending on the experiment). I ran FASTQC in all the samples and found overrepresented adapter sequences:

I removed the adaptors (TruSeq adaptors) using Cutadapt (in addition, I removed low quality and N bases from the 3' end of the reads). After that, I ran again FASTQC and the output is the following (representative example) :

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

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

---
Tags: ngs, fastq, trimming
---


# Question
Title: Improve a reference genome with sequencing data
I have a DNA sample which I know doesn't *quite* match my reference genome - my culture comes from a subpopulation which has undergone significant mutation since the reference was created.

The example I have in mind is E.coli. We've tried assembly using a couple of different tools and the de-novo assembly isn't as high quality as we would like, despite having tonnes of data. Approaching this from a Bayesian point of view, the reference genome provides a very good prior if we could use it wisely

From visual inspection with IGV, a significant number of both SNPs and SVs appear to be present, but an assembly built entirely from my own sequencing data is not high enough quality for my purposes.

How can I modify this reference genome to match my sample with new sequencing data (preferably with Oxford Nanopore Technologies long reads, but I can also use these to scaffold short reads if necessary), taking advantage of my knowledge that the existing reference is *mostly* very good, without having to access the reads which were originally used to construct the reference genome?

The goal of the project isn't to determine where the SVs are, I just need a reference that accurately represents my sample in order to use the data for downstream analysis (as the training set for machine learning.) So by a high quality reference, I mean one which represents as well as possible the sample that was sequenced. To make matters worse, this may not be the one which has the highest alignment identity if there are systematic sequencing errors, as in nanopore sequencing!

# Answer
Depending on the coverage of your data and the complexity of the genome, you could either reassemble the genome *de novo* or run a reference guided (or reference assisted) assembly. It sounds like you're leaning more towards the latter.

There are a couple of reference-guided assembly tools available: AlignGraph and Ragout. These may or may not be appropriate depending on the organism of interest and your data types. For example, these tools are very unlikely to work well on Oxford Nanopore reads which have not been error corrected using Nanopolish or Canu -correct.
> 4 votes

# Answer
One approach to this is to use whatever data you have to iteratively update the reference genome. You can keep chain files along the way so you can convert coordinates (e.g. in gff files) from the original reference to your new pseudoreference.

A simple approach might be:

1. Align new data to existing reference
2. Call variants (e.g. samtools mpileup, GATK, or whatever is best for you)
3. Create new reference incorporating variants from 2
4. Rinse and repeat (i.e. go to 1)

You can track some simple stats as you do this - e.g. the number of new variants should decrease, the number of reads mapped should increase, and the mismatch rate should decrease, with every iteration of the above loop. Once the pseudoreference stabilises, you know you can't do much more.
> 4 votes

# Answer
You can use nanopolish using the illumina reads. Also have a look at pilon.
> 3 votes

# Answer
If it is a short-read draft assembly and you have long-reads (ONT or Pacbio) run links to scaffold the genome and then run Pilon iteratively to try to polish and fill gaps using the short-reads.
> 3 votes

# Answer
Kindel (which I wrote) can infer consensus from low quality alignments of short reads to viral references, and extending it to work with single molecule reads and larger genomes is on my to-do list, though I imagine this will require some redesign.
Presumably you're dealing with a bacterial or fungal genome in this case ? I do also have a basic C++ version, but it's a long way away from being user friendly. Anyway, it may be worth a look – feel free to get in touch with any issues you encounter. I'd use Pilon as mentioned above
> 0 votes

---
Tags: genome, phylogenetics, assembly
---


title: Annovar index out of date

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

> WARNING: Your index file annovar\_latest/humandb/hg19\_cadd.txt.idx is out of date and will not be used. ANNOVAR can still generate correct results without index file.

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

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

---
Tags: sequence-annotation
---


# Question
Title: What are the available cloud computing services for bioinformatics?
I'm looking for cloud computing services that can be used for doing bioinformatics. An example I found is InsideDNA and there is Amazon of course. A little description of these would be appreciated.

# Answer
> 4 votes
It really depends on what you are trying to do, but here are a few services that I know of.

* GATK on Google Genomics Cloud: Google and the Broad offer a cloud instance tailored to GATK pipelines.
* Genomics on Amazon Web Services: I don't think there is anything that makes this unique, but Amazon offers some resources to help get started with genomics/life sciences-centered cloud solutions.
* Illumina Bioinformatics: Illumina is working on a whole suite of bioinformatics software for the cloud.
* Cancer Genomics Cloud: This is specific to cancer genomics, but I believe Seven Bridges allows you to push all sorts of data into the tool and analyze it.

# Answer
> 9 votes
I have trialed the free version of InsideDNA, and these were my notes:

* Cost: \\$225/month for a team of 5 with 50TB storage or \\$45/month with 10TB storage for individuals (assuming 6 month package: https://insidedna.me/pricing).
* Software installed: Around 600 bioinformatic tools available and standard command line tools; some popular tools missing (like CD-HIT), but should be possible to install on request.
* Jobs: Maximum of 32 CPUs and 208 RAM per job submission. Test jobs generally worked, although a larger test job failed.
* Other points: Command line was sometimes slow, `wget` queries were slow, and `scp` was blocked. However, these may be resolvable issues.

Overall, I felt InsideDNA could be useful for groups without their own computational infrastructure and could be used for easily sharing resources between groups. The packages on offer seem not expensive, but I had a few issues, and I don't know how good their sys admin support would be.

I have not used the Amazon service, so can't comment beyond the details on their website. Also there are a few alternative companies, such as Genestack and DNAnexus, but I haven't directly tested them either.

# Answer
> 5 votes
I'm not sure what kinds of bioinformatics tasks you would like to perform, therefore it is difficult to give a good recommendation.

If you're specifically working on statistical genetics, I can recommend Hail \[1\]. Hail is an open-source tool for analyzing genetics data at the tens of terabyte scale. Most of Hail's users do their science in Jupyter notebooks that are backed by Google Cloud Platform Dataproc clusters. Hail permits you to perform a variety of statistical genetics tasks including:

* filtering and aggregation for quality control
* subsetting, linear regression, linear mixed model regression, and linear burden testing
* utilities for computing various measures of relatedness
* principal components analysis
* variant splitting
* import/export from a variety of formats including PLINK, VCF, and BGEN, and
* a python API which enables the use of libraries like matplotlib for plotting analysis results

To learn specifically about using Hail with the Google Cloud Platform and Jupyter notebooks, I strongly recommend Liam's Hail forum post about his cloud-tools repository.

Here's an example, from the Hail tutorial, of using Hail to perform some quality control and display a scatter plot of the first two principal components of the individuals:

```
from hail import *
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

hc = HailContext()

table = hc.import_table('data/1kg_annotations.txt', impute=True).key_by('Sample')
common_vds = (hc.read('data/1kg.vds')
              .annotate_samples_table(table, root='sa')
              .sample_qc()
              .filter_samples_expr('sa.qc.dpMean >= 4 && sa.qc.callRate >= 0.97')
              .filter_genotypes('''let ab = g.ad[1] / g.ad.sum() in
                         ((g.isHomRef && ab <= 0.1) ||
                          (g.isHet && ab >= 0.25 && ab <= 0.75) ||
                          (g.isHomVar && ab >= 0.9))''')
              .variant_qc()
              .filter_variants_expr('va.qc.AF > 0.01')
              .ld_prune(memory_per_core=512, num_cores=4))

pca = common_vds.pca('sa.pca', k=5, eigenvalues='global.eigen')
pca_table = pca.samples_table().to_pandas()

colors = {'AFR': 'green', 'AMR': 'red', 'EAS': 'black', 'EUR': 'blue', 'SAS': 'cyan'}
plt.scatter(pca_table["sa.pca.PC1"], pca_table["sa.pca.PC2"],
            c = pca_table["sa.SuperPopulation"].map(colors),
            alpha = .5)
plt.xlim(-0.6, 0.6)
plt.xlabel("PC1")
plt.ylabel("PC2")
legend_entries = [mpatches.Patch(color=c, label=pheno) for pheno, c in colors.items()]
plt.legend(handles=legend_entries, loc=2)
plt.show()

```

\[1\] Disclaimer: I work on Hail

# Answer
> 4 votes
Depending on your applications and uses, you might be interested in checking out CyVerse. It is an NSF funded initiative that provides you with data storage, high performance computing resources, and easy access to commonly used tools. As far as I know, it is free to use once you have an account. I also usually encounter it being used with plant and microbial genomics, so not sure how it will work with something like human genomics projects. But might be worth checking out at least. :)

More information: http://www.cyverse.org/about

# Answer
> 1 votes
### Google Genomics

Google has an API called Google Genomics.

### SNPedia

"SNPedia is a wiki investigating human genetics." snpedia.com

### Promethease

"Promethease is a literature retrieval system that builds a personal DNA report" promethease.com

### DNA Land

"Compare DNA with reference data from different populations" dna.land

### The CyDAS Project

And, there's the CyDAS project which has an API that can analyze ISCN formulae. Per their web site: their API "lets you analyze a Karyotype for virtually all information which can be extracted from karyotypes and the rearrangements therein: gains and losses of chromosomal material, break points, junctions..." It's a **free** service, but I don't know how up to date it is.

# Answer
> 1 votes
DNAnexus -- http://dnanexus.com

BaseSpace -- http://basespace.illumina.com

Seven Bridges Genomics -- http://www.sbgenomics.com

Curoverse http://curoverse.com

InsideDNA http://insidedna.me/

# Answer
> 0 votes
There are tonnes of them. On top of the excellent ones everybody mentioned

1. iRods
2. Arvados
3. Galaxy

---
Tags: linux
---


# Question
Title: How can I edit a specific FASTQ read in place, given the read ID?
I am introducing SNVs into specific samples in order to estimate false negative rates for a variant calling pipeline. I know reads can be simulated but I would actually prefer to use the real data so as to keep everything else equal.

Given a read ID, I want to edit a single basecall (e.g. the 12th base) for just that read from within a large FASTQ file containing millions of reads.

example, i want to change the 12th base ('C') in read 31027 to a 'T':

```
@70630 1:N:0:ATCACG
GAAGGTCCATGGATAATACTCAATTTTCCACAACAGCTTTTGTACTCTAGATCATTGATATTTACCAAAAGTCACTCAAACTCATCCTATGCATAATTCTAGTCCACCAATCATGATATGATGGAGAACATGGTTGTAATCAGGAAGACAG
+
DDDDDHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIGHHIIIIIIIIIIIIIIIIIIIIIIHIIHIHIIIHHIGIIIIHHHHIIH@
@31027 1:N:0:ATCACG
CAAAAGTCACTCAAACTCATCCTATGCATAATTCTAGTCCACCAATCATGATATGATGGAGAACATGGTTGTAATCAGGAAGACAGATAAAGCAGCAGACCAAAAGTAATCTGAGAAATTATATTTGAATCACTCAGATATACATCAAATA
+
DDBDDHC@GHFHHIIE@CEHHIHEFFCGFCH?HHEGFCCEHHCGH@HFHHHCHHIIIHEHIIHII@CECHIHIIIEECCDGEHFFHHHHEHHFHHGGIHHHDGHFHIIIGHHHHHHEHHIIHIIIGGHHHCHHHCHIIIHHIEH@GHIHIC
@87319 1:N:0:ATCACG
CAATTAAGCTTTGGCAACGGTGGTCAAGATGAGATGCATATGGAGATAATAACTAAAAGTCAATCGAGACTCATCGTATGCATATTTCTAGTCCATCGATCATGAAATGATAGGATAGCTAGAATGAAAAGTAAATTTCCAGAAGGTCCAT
+
DDDDDIIIIIIIIIIIIIHIHHIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIHHIIIIIIIIIIIIIIIIIHIIHIIIIIIIIIIIIHIIIIIIIIIHIHIH

```

Of course I could stream through the entire FASTQ file rewriting everything to a new file until i get to the read of interest, do the edits to that read, write it out to the new file, and then keep going until the end. BioPython would work fine for this as would other approaches.

However, is there an efficient way to do this without the full read/write stream? Could this be done with sed/awk ? What about an indexed fastq file?

cheers David

# Answer
> 5 votes
You'll first need to determine the appropriate line number, which you can do with `grep -m1 -nw "@31027" foo.fastq`. After that, note that you can provide a line number to `sed`:

```
sed -i '123456s/CAAAAGTCACTCA/CAAAAGTCACTTA/' foo.fastq

```

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

# Answer
> 2 votes
If your fastq file is simple, if each of your reads only has a single line of DNA, you could do:

```
awk '{
        if($1 == "@31027"){
            a=NR+1
        }
        if(NR==a){
            split($0,s,"");
            s[12]="T";
            for(i in s){
                seq = sprintf("%s%s",seq,s[i]);
            }
            $0=seq
        }
    }1;' file.fastq

```

---
Tags: python, ngs, fastq, awk
---


title: Subset FASTA file by species name

text:
# Question
Title: Subset FASTA file by species name
I have a problem: I've managed to download a massive fasta file of 1500 sequences, but now I want to split them into separate fasta files based on the genus.

EDIT
The fasta file looks like this:

```
terminase_large.fasta
>YP_009300697.1 terminase large subunit [Arthrobacter phage Mudcat]
MGLSNTATPLYYGQF...

>YP_009208724.1 hypothetical protein ADP65_00072 [Achromobacter phage phiAxp-3]
MSNVLLKQELDEWL...

>YP_009148449.1 large terminase subunit [Delftia phage RG-2014]
MSEPRKLVKKTLD...
...

```

I would like to end up with something like this:

```
mycobacterium_terminase.fasta
>Mycobacterium phage JAMaL
MVRKKPPPELE...

>Mycobacterium phage Bruin
MEVCGYTLDDI...

>Mycobacterium phage Zaka
MSLDNHLPELA...

salmonella_terminase.fasta
>Salmonella virus SETP7
MNVDITATEPQ...

>Salmonella virus SE2
MEGGSDSLIAM...
and so on...

```

So that the label on the alignment would have the name of the phage and not the protein's ID.

I've managed to extract the genus names of my organism with this:

```
outfile = open('species.txt', 'w')
with open('terminase_large.fasta') as fd:
    for line in fd:
        if line.startswith('>'):
            if '[' in line:
                name=line.split('[')[-1]
                name=name.split(' ', 1)[0]
                outfile.write(name[:] + "\n")
outfile.close()

```

And got to extract only the unique names with this:

```
lines_seen = set()
outfile = open('species2.txt', "w")
for line in open("species.txt", "r"):
    if line not in lines_seen:  # not a duplicate
        outfile.write(line)
        lines_seen.add(line)
outfile.close()

```

(Can I merge those two scripts together?)

Now, my genus names look like this:

```
Arthrobacter
Achromobacter
Delftia
....

```

I tried automating my script to get the Entrez data, but it gives me the 'Supplied id parameter is empty' message

My code looks like this:

```
from Bio import Entrez
Entrez.email = "example@example.org"

for line in open("species2.txt", "r"):
    searchterm = "(terminase large subunit AND viruses[Organism]) AND" +line+ "AND refseq[Filter]"
    searchResultHandle = Entrez.esearch(db="protein", term=searchterm, retmax=1000)
    searchResult = Entrez.read(searchResultHandle)
    ids = searchResult["IdList"]

    handle = Entrez.efetch(db="protein", id=ids, rettype="fasta", retmode="text")
    record = handle.read()

    out_handle = open('terminase_large_'+str(line[:-1])+'.fasta', 'w')
    out_handle.write(record.rstrip('\n'))

```

Can someone help me with it?

# Answer
Splitting into multiple files and changing the IDs can be easily done:

```
perl -pe 'if(/>/){/\[(.*?)\]\s*$/; $_="> $1\n"}' file.fa |
    awk '(/^>/){name=$2} {print >> name".fa"}'

```

That assumes all your FASTA headers have `[foo bar baz]` as the last element of a line. It will create a file called `foo.fa` (the bacterium's name) with all sequences saved there.

### Explanation

* `perl -pe` : run the script given by `-e` on each line of the input file, and print the resulting line.
* `if(/>/)` : if this line starts with a `>`. \
* `/\[(.*?)\]\s*$/` : match an opening bracket (`\[`), then capture (that's what the parentheses do, they capture a pattern so we can refer to it as `$1`) everything until the first `]` (`.*?\]`)
* `$_="> $1\n"` : the `$_` special variable in Perl is (in this case) the current line. So, `$_=foo` means "make the current line read `foo`. Since the `-p` prints each input line, changing the value of `$_` means the changed value will be printed. So here, we are printing `>`, whatever was in the square brackets (`$1`) and a newline character.

The output of the perl command alone on your example input file is:

```
$ perl -pe 'if(/>/){/\[(.*)\]\s*$/; $_="> $1\n"}' file.fa
> Arthrobacter phage Mudcat
MGLSNTATPLYYGQF...

> Achromobacter phage phiAxp-3
MSNVLLKQELDEWL...

> Delftia phage RG-2014
MSEPRKLVKKTLD...

```

So, we now pass it through `awk` which does:

* `(/^>/){a=$2}` : if this line starts with an `>`, save the second field (the bacterial species) as the variable `name`.
* `{print >> name".fa"}` : print each line into a file whose name is the current value of the variable `name` with a `.fa.` extension.

---

If you prefer python scripts to the one-liner approach, you can do the same thing with:

```
#!/usr/bin/env python
import re
outFile = None
with open("file.fa", "r") as inFile, open(species + ".fa", 'a') as outFile:
    for line in inFile:
        line = line.rstrip()
        if line.startswith('>'):
            regex = re.compile('.*\[((.+?)\s+.*?)\].*')
            matches = regex.search(line)
            species = matches[2]
            outFile.write('>%s\n' % matches[1])
        else:
            outFile.write("%s\n" % line)

```

---

As for your script, you've got the right idea, but have a small bug. You forgot to remove the `\n` from your input file, so it looks for `Arthrobacter\n` instead of `Arthrobacter`. The golden rule of debugging is "print all the things". If you add `print("Searchterm: ",searchterm)` you will see:

```
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDArthrobacter
AND refseq[Filter]
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDAchromobacter
AND refseq[Filter]
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDDelftia
AND refseq[Filter]

```

So, you need to remove the newline characters and add a space like so (I also made it a bit more "pythonic" and conforming to the Python syntax guidelines):

```
#!/usr/bin/env python
from Bio import Entrez
Entrez.email = "example@example.org"

with open("species2.txt", "r") as in_handle:
    for line in in_handle:
        line = line.rstrip()
        searchterm = ("(terminase large subunit AND viruses[Organism]) " +
                      "AND %s AND refseq[Filter]" % line)
        print("Searchterm: ", searchterm)
        searchResultHandle = Entrez.esearch(db="protein",
                                            term=searchterm, retmax=1000)
        searchResult = Entrez.read(searchResultHandle)
        ids = searchResult["IdList"]

        handle = Entrez.efetch(db="protein", id=ids,
                               rettype="fasta", retmode="text")
        record = handle.read()

        with open('terminase_large_' + str(line[:-1]) + '.fasta', 'w') as out_handle:
            out_handle.write(record.rstrip('\n'))

```
> 7 votes

# Answer
## A (bio)awk-based solution

### Generate the list of genus

```
bioawk -c fastx '{print $comment}' terminase_large.fasta \
    | sed -r 's/.*\[(\w+).*/\1/' \
    | sort -u > genus_names.txt

```

Here is how it works:

1. `bioawk` is like awk but with extra parsing capabilities. `-c fastx` parses the fields in fasta or fastq records. `{print $comment}` will print the `$comment` field which contain the genus names you want to extract.
2. The `sed` command captures (with parentheses) the word (`\w+`) that is just after the opening square bracket (which should be the genus name, if there is only one opening square bracket on the line), and substitutes the whole line (that is, everything before the bracket (`.*`), the bracket (`\[`), the word, and everything after the word (`.*` again)) with just the captured word (`\1`).
3. All this is sorted, and the `-u` option of `sort` ensures there is only one occurrence of each genus name in the output, which we put in the `genus_names.txt` file.

### Extract the fasta records for each genus

```
for genus in $(cat genus_names.txt)
do
    bioawk -c fastx '{print}' terminase_large.fasta \
        | grep ${genus} \
        | awk -F "\t" '{print ">"$4"\n"$2}' \
        | sed -r 's/^>[^\[]+\[(.*)\]/>\1/' > ${genus}.fasta
done

```

Here is how it works:

1. We loop on the content of the previously established list of genus names.
2. For each value of `${genus}`, we parse again the fasta file with `bioawk`. This time, for a given fasta record, we print all the elements that have been parsed. They appear on one line, in the following order (according to `bioawk -c help`): `1:name 2:seq 3:qual 4:comment` (we actually only need `seq` and `comment` for the following steps, but there are less risks of making errors at the later `awk` step if we keep the four fields).
3. We use `grep` to select the resulting lines that contain the genus name we're currently dealing with.
4. The lines that have been successfully selected are then parsed by awk and reformatted into fasta format using only the comment part in the header. The `-F "\t"` is to indicate that the field delimiters are tabulations (which is how bioawk has written the records when we did `{print}`).
5. We use `sed` to keep only the part inside the square brackets (again using a pattern capture between parentheses, that we use in the substitution part with `\1`), and the output is written in a file named using the current value of `${genus}`.

More details on the last `sed` command:

We match the whole header lines. They start with "\>" (`^>`), then there are some non-"\[" characters (`[^\[]+`), then a "\[" (`\[`), then something we capture (`(.*)`), then a "\]" (`\]`). We replace them with the captured part. The sequence lines are unaffected because they don't match the above regular expression.
> 2 votes

---
Tags: fasta, python
---


title: Get canonical transcript from UCSC

text:
# Question
Title: Get canonical transcript from UCSC
I am using the following command to get all refseq genes from UCSC:

```
/usr/bin/mysql --user=genomep --password=password --host=genome-mysql.cse.ucsc.edu \
    -A -D hg38 -e 'select concat(t.name, ".", i.version) name, \
    k.locusLinkId as "EntrezId", t.chrom, t.strand, t.txStart, \
    t.txEnd, t.cdsStart, t.cdsEnd, t.exonCount, t.exonStarts, \
    t.exonEnds, t.score, t.name2 from refGene t join hgFixed.gbCdnaInfo i \
    on t.name = i.acc join hgFixed.refLink k on t.name = k.mrnaAcc'

```

That returns data in the following format (showing the 1st 5 lines):

```
+-------------+-----------+-------+--------+---------+-------+----------+--------+-----------+--------------------+--------------------+-------+-----------+
| name        | EntrezId  | chrom | strand | txStart | txEnd | cdsStart | cdsEnd | exonCount | exonStarts         | exonEnds           | score | name2     |
+-------------+-----------+-------+--------+---------+-------+----------+--------+-----------+--------------------+--------------------+-------+-----------+
| NR_046018.2 | 100287102 | chr1  | +      |   11873 | 14409 |    14409 |  14409 |         3 | 11873,12612,13220, | 12227,12721,14409, |     0 | DDX11L1   |
| NR_106918.1 | 102466751 | chr1  | -      |   17368 | 17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-1 |
| NR_107062.1 | 102465909 | chr1  | -      |   17368 | 17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-2 |
| NR_107063.1 | 102465910 | chr1  | -      |   17368 | 17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-3 |
| NR_128720.1 | 103504738 | chr1  | -      |   17368 | 17436 |    17436 |  17436 |         1 | 17368,             | 17436,             |     0 | MIR6859-4 |
+-------------+-----------+-------+--------+---------+-------+----------+--------+-----------+--------------------+--------------------+-------+-----------+

```

I also want to find the accession of the canonical transcript of each of the genes returned by the command above. Those seem to be stored in the `knownCanonical` table:

```
$ /usr/bin/mysql --user=genomep --password=password --host=genome-mysql.cse.ucsc.edu -A -D hg38 -e 'select * from knownCanonical limit 5'
+-------+------------+-----------+-----------+------------+--------------------+
| chrom | chromStart | chromEnd  | clusterId | transcript | protein            |
+-------+------------+-----------+-----------+------------+--------------------+
| chrX  |  100628669 | 100636806 |         1 | uc004ega.3 | ENSG00000000003.14 |
| chrX  |  100584801 | 100599885 |         2 | uc004efy.5 | ENSG00000000005.5  |
| chr20 |   50934866 |  50958550 |         3 | uc002xvw.2 | ENSG00000000419.12 |
| chr1  |  169853073 | 169893959 |         4 | uc001ggs.5 | ENSG00000000457.13 |
| chr1  |  169795048 | 169854080 |         5 | uc001ggp.4 | ENSG00000000460.16 |
+-------+------------+-----------+-----------+------------+--------------------+

```

However, there seems to be no obvious way to link the `knownCanonical` table to the `refGene`, `hgFixed.gbCdnaInfo` and `hgFixed.refLink` tables used above.

So, how can I modify my 1st query (or, if necessary write a new one) so that my results also include the accession of the gene's canonical transcript?

# Answer
You can link them with the `kgXref` table, since `kgXref.refseq == refGene.name` and `kgXref.kgID == knownCanonical.transcript`. Since it seems that `knownCanonical.transcript` is what you want anyway, you don't even need to join on it:

```
mysql --user=genomep --password=password --host=genome-mysql.cse.ucsc.edu \
-A -D hg38 -e 'select concat(t.name, ".", i.version) name, x.kgID, \
k.locusLinkId as "EntrezId", t.chrom, t.strand, t.txStart, \
t.txEnd, t.cdsStart, t.cdsEnd, t.exonCount, t.exonStarts, \
t.exonEnds, t.score, t.name2 from refGene t join (hgFixed.gbCdnaInfo i, hgFixed.refLink k, kgXref x) \
on (t.name = i.acc and t.name = k.mrnaAcc and t.name = x.refseq) limit 10'

```

The output is then:

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

```

If you want the protein ID then join `knownCanonical` on `x.kgID`.
> 3 votes

---
Tags: public-databases, identifiers, ucsc
---


# Question
Title: BAM to BigWig without intermediary BedGraph
I have a pipeline for generating a BigWig file from a BAM file:

```
BAM -> BedGraph -> BigWig

```

Which uses `bedtools genomecov` for the `BAM -> BedGraph` part and `bedGraphToBigWig` for the `BedGraph -> BigWig` part.

The use of `bedGraphToBigWig` to create the `BigWig` file requires a BedGraph file to reside on disk in uncompressed form as it performs seeks. This is problematic for large genomes and variable coverage BAM files when there are more step changes/lines in the `BedGraph` file. My `BedGraph` files are in the order of 50 Gbytes in size and all that IO for 10-20 BAM files seems unnecessary.

**Are there any tools capable generate `BigWig` without having to use an uncompressed `BedGraph` file on disk? I'd like this conversion to hapen as quickly as possible.**

I have tried the following tools, but they still create/use a `BedGraph` intermediary file:

* `deepTools`

# Some Benchmarks, Ignoring IO

Here are some timings I get for creating a `BigWig` file from a `BAM` file using 3 different pipelines. All files reside on a `tmpfs` i.e. in memory.

## BEDTools and Kent Utils

This is the approach taken by most.

```
time $(bedtools genomecov -bg -ibam test.bam -split -scale 1.0 > test.bedgraph \
  && bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
  && rm test.bedgraph)

real    1m20.015s
user    0m56.608s
sys     0m27.271s

```

## SAMtools and Kent Utils

Replacing `bedtools genomecov` with `samtools depth` and a custom `awk` script (`depth2bedgraph.awk`) to output `bedgraph` format has a significant performance improvement:

```
time $(samtools depth -Q 1 --reference test.fasta test.bam \
  | mawk -f depth2bedgraph.awk \
  > test.bedgraph \
  && bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
  && rm test.bedgraph)

real    0m28.765s
user    0m44.999s
sys     0m1.166s

```

Although it is has less features, we used `mawk` here as it's faster than `gawk` (we don't need those extra features here).

### Parallelising with xargs

If you want to have a `BigWig` file per chromosome, you can easily parallelise this across chromosome/reference sequences. We can use `xargs` to run 5 parallel BAM-\>BedGraph-\>BigWig pipelines, each using the `tmpfs` mounted `/dev/shm` for the intermediary `BedGraph` files.

```
cut -f1 test.fasta.chrom.sizes \
  | xargs -I{} -P 5  bash -c 'mkdir /dev/shm/${1} \
  && samtools depth -Q 1 --reference test.fasta -r "${1}" test.bam \
  | mawk -f scripts/depth2bedgraph.awk \
  > "/dev/shm/${1}/test.bam.bedgraph" \
  && mkdir "./${1}" \
  && bedGraphToBigWig \
    "/dev/shm/${1}/test.bam.bedgraph" \
    test.fasta.chrom.sizes \
    "./${1}/test.bam.bw" \
  && rm "/dev/shm/${1}/test.bam.bedgraph"' -- {}

```

## deepTools

Let's see how `deepTools` performs.

```
time bamCoverage --numberOfProcessors max \
  --minMappingQuality 1 \
  --bam test.bam --binSize 1 --skipNonCoveredRegions \
  --outFileName deeptools.bw

real    0m40.077s
user    3m56.032s
sys     0m9.276s

```

# Answer
> 11 votes
This can be done in R very easily from an indexed .bam file.

Given single-end file for **sample1**.

```
library(GenomicAlignments)
library(rtracklayer)

## read in BAM file (use readGAlignmentPairs for paired-end files)
gr <- readGAlignments('sample1.bam')

## convert to coverages
gr.cov <- coverage(gr)

## export as bigWig
export.bw(gr.cov,'sample1.bigwig')

```

Be aware that this method doesn't include normalization steps (such as normalizing to total coverage). Most of these additional steps can be added if necessary.

# Answer
> 4 votes
Just for reference, there is a not very in depth, but first-hand reasoning (Jim Kent) given here about why bedGraphToBigWig does not support streaming http://genome.soe.ucsc.narkive.com/2S1Z3VpG/bedgraphtobigwig-reading-in-from-stdin

Edit: as noted wigToBigWig allows bedgraph streaming input via stdin but takes (possibly much) more memory than bedGraphToBigWig itself does

```
cat file.bedgraph | wigToBigWig stdin chrom.sizes output.bigwig

```

Note: You can also use `wiggletools bin 10` for example to bin rows and downsample

# Answer
> 2 votes
The QuasR package in R has a useful function called `qExportWig`.

To use it you first need to create a sample sheet and a `qProject` object - this can be done starting with bam files or fastq files (in which case it also performs alignments). You can then create wig or bigWig files without writing a bedGraph file to disk.

An example would look something like this:

```
library(QuasR)

# Create qProject object using qAlign
# Sample file should contain tab-separated file names (BAM or fastq) and sample names
myproj <- qAlign("samples.txt", genome=mygenome, paired="no")

# Export all samples in project as bigWig, scaling to mean number of alignments by default
qExportWig(myproj, binsize=50, createBigWig=T)

```

It has a lot of flexible options, including scaling or log-transforming counts, shifting reads and changing the bin size.

---
Tags: sam, file-formats, format-conversion
---


# Question
Title: Convert a BAM file from one reference to another?
I have a set of BAM files that are aligned using the NCBI GRCh37 human genome reference (with the chromosome names as NC\_000001.10) but I want to analyze it using a BED file that has the UCSC hg19 chromosome names (e.g. chr1). I want to use bedtools to pull out all the on-target and off-target reads.

1. Are NCBI and UCSC directly comparable? Or do I need to re-align the BAM/lift-over the BED to the UCSC reference?
2. Should I convert the BED file or the BAM file? Everyone here uses the UCSC chromosome names/positions so I'll need to convert the eventual files to UCSC anyway.

# Answer
> 24 votes
You're the second person I have ever seen using NCBI "chromosome names" (they're more like supercontig IDs). Normally I would point you to a resource providing mappings between chromosome names, but since no one has added NCBI names (yet, maybe I'll add them now) you're currently out of luck there.

Anyway, the quickest way to do what you want is to `samtools view -H foo.bam > header` to get the BAM header and then change each NCBI "chromosome name" to its corresponding UCSC chromosome name. DO NOT REORDER THE LINES! You can then use `samtools reheader` and be done.

Why, you might ask, would this work? The answer is that chromosome/contig names in BAM files aren't stored in each alignment. Rather, the names are stored in a list in the header and each alignment just contains the integer index into that list (read group IDs are similar, for what it's worth). This also leads to the warning above against reordering entries, since that's a VERY convenient way to start swapping alignments between chromosomes.

As an aside, you'd be well served switching to Gencode or Ensembl chromosome names, they're rather more coherent than the `something_random` mess that's present in hg19 from UCSC.

**Update**: Because I'm nice, here is the conversion between NCBI and UCSC. Note that if you have any alignments to patches that there is simply no UCSC equivalent. One of the many reasons not to use UCSC (avoid their annotations too).

# Answer
> 4 votes
The "right" solution would be realignment, but that's expensive and most of us would not go that route. My preferred solution would be to convert the bed file, as opposed to the bam. Here's why:

1) Reheadering the bam means that you may have reads aligned to contigs without a corresponding entry in UCSC (see Devon's list for the mappings). This is a problem because:

* Some of those reads would likely have been mapped elsewhere if a reference without those contigs was used.
* I'm not even sure what happens to those reads after reheadering - I guess they would need to be marked as unmapped? Lots of potential for screwiness there.

2) It seems cleaner to convert the bed file from UCSC-\>NCBI, where you are guaranteed that every entry has a "home". Then, after you pull your info from the bam, you can always convert chromosome names back if you need to.

# Answer
> 0 votes
Just to point out that if you want to follow @Devon Ryan's answer for a different organism/assembly, that is not in his very useful linked resource, you can download NCBI to UCSC contig to chromosome number mappings from https://www.ncbi.nlm.nih.gov/assembly.

To the site and search for your assembly name. At the bottom of the page is a box called "Global assembly definition" containing a link titled "Download full sequence report".

The downloaded file contains a table with:

* Chromosome numbers in "Sequence-Name"/"Assigned-molecule"
* NCBI names in Refseq-Accn
* UCSC contig name in "UCSC-style-name"

---
Tags: sam, file-formats, bed, format-conversion, reference-genome
---


title: Cannot install chromosomer

text:
# Question
Title: Cannot install chromosomer
I am trying to install chromosomer but I fail. Can anybody help me, please?

```
$ pip install chromosomer

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

```

The output `python`:

```
$ python

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

```

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

```
pip2 install chromosomer

```

---

**But:**

Your Python installation is a bit screwed up: there are Python 2 and Python 3, which are unfortunately incompatible.

On your system `python` and `pip` seem to be aliases for Python 3, which I’d strongly recommend against (lots of tools will break). Instead, they should alias Python 2. Python 3 should be invoked via `python3`/`pip3`.
> 6 votes

---
Tags: python, software-installation
---


# Question
Title: How do you write a .gz fastq file with Biopython?
How do you write a `.gz` (or `.bgz`) `fastq` file using `Biopython`?

I'd rather avoid a separate system call.

The typical way to write an `ASCII` `.fastq` is done as follows:

```
for record in SeqIO.parse(fasta, "fasta"):
    SeqIO.write(record, fastq, "fastq")

```

The `record` is a `SeqRecord` object, `fastq` is the `file handle`, and `"fastq"` is the requested file format. The file format may be `fastq`, `fasta`, etc., but I do not see an option for `.gz`.

Here is the SeqIO API.

# Answer
I'm not sure I'm doing it the best way, but here is an example where I read a compressed gzip fastq file and write the records in block gzip fastq:

```
from Bio import SeqIO, bgzf
# Used to convert the fastq stream into a file handle
from io import StringIO
from gzip import open as gzopen

records = SeqIO.parse(
    # There is actually simpler (thanks @peterjc)
    # StringIO(gzopen("random_10.fastq.gz").read().decode("utf-8")),
    gzopen("random_10.fastq.gz", "rt"),
    format="fastq")

with bgzf.BgzfWriter("test.fastq.bgz", "wb") as outgz:
    SeqIO.write(sequences=records, handle=outgz, format="fastq")

```
> 19 votes

# Answer
Here is how I was able to export to a standard gzipped FASTQ file using Biopython. Basically, instead of using `SeqIO.write`, I directly called the `.format` method of the SeqRecord object. The example code below imports a gzipped FASTQ file, removes reads that do not contain a G in positions 7, 8, and 9, and writes the results to a gzipped FASTQ file.

```
from Bio import SeqIO
import gzip

path_in = "path/to/in.fastq.gz"
path_out = "path/to/out.fastq.gz"
handle_in = gzip.open(path_in, "rt")
handle_out = gzip.open(path_out, "wt")

fq = SeqIO.parse(handle_in, "fastq")
for read in fq:
    # Only export reads that have a G in positions 7, 8,
    # and 9
    if read.seq[6:9] == "GGG":
        handle_out.write(read.format("fastq"))

handle_in.close()
handle_out.close()

```

Using the code example from the original question that converts a FASTA file to a FASTQ file, it would look something like the following:

```
fastq = gzip.open("path/to/out.fastq.gz", "wt")
for record in SeqIO.parse(fasta, "fasta"):
    fastq.write(record.format("fastq"))
fastq.close()

```

**Caveat:** This only works for Python 3. Python 3 differs from Python 2 in how it imports data from binary files (source). In fact, ideally you should be using Python 3.4 or greater since Biopython has compatibility issues with Python 3.3 (source).
> 3 votes

# Answer
As far as I know you can not write compressed fastq on the fly. You will have to write the entire fastq and then compress it.
> 0 votes

---
Tags: file-formats, python
---


# Question
Title: Major variability in speed of BLAST between proteins
I'm running blast+ searches for a list of proteins using specific parameters.

This is my search

```
 blastp -query protein.fasta -remote -db refseq_protein -out
protein_refseq_blast.txt -outfmt '6 qseqid sseqid qstart qend length evalue' -
entrez_query 'txid2 [ORGN] or txid4751 [ORGN]' -evalue 200000 -max_target_seqs 10000

```

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

Thanks.

# Answer
Try changing you `evalue` threshold to something like `0.001`.

If you blasted a sequence against a database containing 10,000 sequences and got a hit with an evalue of `1`, this would mean that, given the size of your database and length of protein, you would expect to find a match as good as (or better than) the hit returned once purely by chance. If the `evalue` threshold in your query was set to `0.01` (for example) then this hit would not be returned.

If you got a hit with an `evalue` of 200,000 this would mean you would expect to see 200,000 matches with a similar score simply by chance. If your `evalue` threshold were set to 200000, then all of the hits with an evalue \<=200000 would be considered significant. If you set it to `0.001` then only those hits with `evalues` \<= `0.001` would be considered significant.

This means that (depending on the size of the data base) almost every sequence in your query set will be considered as significant given your threshold, which will slow down your search hugely.

See the the section on E values here
> 5 votes

# Answer
The evalue associated with a blast HSP (High-scoring Segment Pair) is a measure of how often you would expect to find such an HSP in a database of this size by pure chance. For instance, if you blast the sequence `ACTG` against the entirety of the NCBI's `nr` database, you would expect several hundred thousand hits by pure chance. That is, several hundred thousand sequences that can be aligned to your query but are in no way related to it (not actual homologs). Those hits don't actually mean anything, so you use the evalue to limit the results to those that are unlikely to be found by chance and, therefore, more likely to have biological meaning if hits are found.

So, setting the evalue to 200000 means you are interested even in HSPs whose characteristics (score, length, etc) mean you would expect to find \<=200000 HSPs with such characteristics by pure chance. This evalue is pointless.

Because your e-value is so high, BLAST will report absolutely anything and everything it finds, and the vast majority of your results will be meaningless and irrelevant. If you want a permissive e-value, you could use `1` or even `10` but 200000 is way too much. You usually only want negative e-values, cases where you wouldn't expect to find even one such HSP by chance. For highly dissimilar sequences, I would recommend you use 0.1 or something of that scale.

Changing the evalue will speed things up considerably. That said, you should always expect different times for different queries since the time it takes to find hits will always depend on the number of hits that can be found. If your sequence can be aligned to many of the sequences in the database, it will take longer to fetch all results.
> 3 votes

---
Tags: blast
---


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

# Answer
> 4 votes
There are more R packages available that calculate GC content, for example Ape's `GC.content()` function.

For example:

```
> library(ape)
> data(woodmouse)
> GC.content(woodmouse)
[1] 0.3873347

```

With a sliding window is in `Biostrings` package.

```
> library(Biostrings)
> DNA <- DNAString("ACTGAAACCGTGGCAGTTTGAC")
> letterFrequencyInSlidingView(DNA, view.width=4,letters="CG")

```

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

---
Tags: r, software-recommendation, phylogenetics
---


# Question
Title: How to calculate overall reference coverage with MUMmer?
Is the MUMmer suite capable of calculating reference sequence coverage statistics for all query sequences collectively? It would be possible to achieve by parsing the output of `nucmer` / `show-coords` / `show-tiling` but it seems like there should be a better way. I currently do this using a sensitive read mapper, `samtools depth` and some scripting.

To clarify, I'd like to know the reference coverage achieved using all of the query sequences (i.e. the whole *de novo* assembly, in my case).

# Answer
> 4 votes
I believe all you need to do is to run `dnadiff` from MUMmer. That will run a comparison and output a number of useful metrics.

---
Tags: assembly, coverage
---


title: IGV faster loading remote bigWig files

text:
# Question
Title: IGV faster loading remote bigWig files
I noticed remote bigWig files (http URL) take sometimes about 1 minute to load on IGV. Once they loaded, jumping from one region to another can also take several seconds.

My understanding is that bigWig files have an internal index for different chunks of different sizes. Can the index be modified so that the speed of loading on IGV is improved? Does this need to happen at the point of creation of the bigWig file, or can the index of an existing bigWig be updated to speed up loading?

# Answer
You would need to do this when creating the bigWig file. Note that I'm not aware of anything that allows manually setting the various zoom levels (libBigWig/pyBigWig will allow you to set a max number, but that's it).

This then raises the question of whether IGV is using an appropriate zoom level to begin with. I've never looked at how it chooses which of the zoom levels to use, but I know it uses them since if you don't have any it will fail to load bigWig files.

Honestly, though, if it's taking a minute to jump to a new region then I suspect that the network is just slow or the remote server is overtaxed. Changing the internal indexing isn't going to fix that.
> 1 votes

---
Tags: bigwig, deeptools
---


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

My question has two parts:

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

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

# Answer
The earliest mention of the 30x paradigm I could find is in the original Illumina whole-genome sequencing paper: Bentley, 2008. Specifically, in Figure 5, they show that most SNPs have been found, and that there are few uncovered/uncalled bases by the time you reach 30x:

These days, 30x is still a common standard, but large-scale germline sequencing projects are often pushing down closer to 25x and finding it adequate. Every group doing this seriously has done power calculations based on specifics of their machines and prep (things like error rates and read lengths matter!).

Cancer genomics is going in the other direction. When you have to contend with purity, ploidy, and subclonal populations, much more coverage than 30x is needed. Our group showed in this 2015 paper that even 300x whole-genome coverage of a tumor was likely missing real rare variants in a tumor.

On the whole, the sequence coverage you need really depends on what questions you're asking, and I'd recommend that anyone designing a sequencing experiment consult with both a sequencing expert and a statistician beforehand (and it's even better if those are the same person!)
> 30 votes

# Answer
Solexa Inc. sequenced NA12878 chrX to ~30x in early 2007, which later became part of Bentley (2008). This, I believe, was the first time that 30x showed up. I don't recall they had a particular reason for that. Figure 5 in the published paper was just aftermath. It does not really explain why not 25x or 35x, given that the curves between 25x and 35x in that figure are about linear.

In the abstract of Ajay et al (2011), the authors argued "the current recommendation of ~30x coverage is not adequate". Nonetheless, the discussion section seems to suggest 50–60x would be necessary with GAIIx, but 35x was adequate with HiSeq2000 plus better recent chemistry. Overall, this paper provides a more thorough analysis. The data quality at that time is also closer to data we produce today.

The required coverage is largely determined by two factors: read placement bias (e.g. GC bias) and base/mapping error rate. While GC bias has been reduced with the PCR-free protocol, base error rate has been going the downward since HiSeq2500. I guess 30x coverage would be necessary if you want to achieve the sensitivity with the older 30x data. Illumina, as a sequencing service provider, and our sequencing facility still insist on the 30x threshold.
> 7 votes

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

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

source: Wikipedia: Rule of three (statistics)

Similarly, you can search for related questions like this one:

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

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

---
Tags: genome, ngs, phylogenetics
---


# Question
Title: Tool or script to parse annotated VCF files
I've annotated VCF files using snpEff and looking for a tool or script to parse the VCF file and clean up the file to make it interpretable for a biologist.

# Answer
(edit)you can filter the VCF annotations with **snpsift**, I've also written a **VcfFilterSequenceOntology** http://lindenb.github.io/jvarkit/VcfFilterSequenceOntology.html

I've written vcf2table: http://lindenb.github.io/jvarkit/VcfToTable.html It decodes VEP and SNPeff annotations:

```
>>chr1/10001/T (n 1)
 Variant
 +--------+--------------------+
 | Key    | Value              |
 +--------+--------------------+
 | CHROM  | chr1               |
 (....)
 VEP
 +--------------------------+------+----------------+------------+-----------------+--------+------------------+-----------------------------------------------+-------------+---------+-----------------+----------------------+
 | PolyPhen                 | EXON | SIFT           | ALLELE_NUM | Gene            | SYMBOL | Protein_position | Consequence                                   | Amino_acids | Codons  | Feature         | BIOTYPE              |
 +--------------------------+------+----------------+------------+-----------------+--------+------------------+-----------------------------------------------+-------------+---------+-----------------+----------------------+
 | probably_damaging(0.956) | 8/9  | deleterious(0) | 1          | ENSG00000102967 | DHODH  | 346/395          | missense_variant                              | R/W         | Cgg/Tgg | ENST00000219240 | protein_coding       |
 |                          | 3/4  |                | 1          | ENSG00000102967 | DHODH  |                  | non_coding_exon_variant&nc_transcript_variant |             |         | ENST00000571392 | retained_intron      |
 |                          |      |                | 1          | ENSG00000102967 | DHODH  |                  | downstream_gene_variant                       |             |         | ENST00000572003 | retained_intron      |
 |                          |      |                | 1          | ENSG00000102967 | DHODH  |                  | downstream_gene_variant                       |             |         | ENST00000573843 | retained_intron      |
 |                          |      |                | 1          | ENSG00000102967 | DHODH  |                  | downstream_gene_variant                       |             |         | ENST00000573922 | processed_transcript |
 |                          |      |                | 1          | ENSG00000102967 | DHODH  | -/193            | intron_variant                                |             |         | ENST00000574309 | protein_coding       |
 | probably_damaging(0.946) | 8/9  | deleterious(0) | 1          | ENSG00000102967 | DHODH  | 344/393          | missense_variant                              | R/W         | Cgg/Tgg | ENST00000572887 | protein_coding       |
 +--------------------------+------+----------------+------------+-----------------+--------+------------------+-----------------------------------------------+-------------+---------+-----------------+----------------------+
 Genotypes
 +---------+------+-------+----+----+-----+---------+
 | Sample  | Type | AD    | DP | GQ | GT  | PL      |
 +---------+------+-------+----+----+-----+---------+
 | M10475  | HET  | 10,2  | 15 | 10 | 0/1 | 25,0,10 |
 | M10478  | HET  | 10,4  | 16 | 5  | 0/1 | 40,0,5  |
 | M10500  | HET  | 10,10 | 21 | 7  | 0/1 | 111,0,7 |
 | M128215 | HET  | 15,5  | 24 | 0  | 0/1 | 49,0,0  |
 +---------+------+-------+----+----+-----+---------+

```
> 5 votes

---
Tags: variant-calling, snp
---