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

I tried `Regtools`:

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

Exact commands I used for this are:

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

The alignments have been generated with Minimap2.

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

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

My questions are:

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

# Answer
You can use the script that comes with minimap2:

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

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

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

```

Some example output:

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

```

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

# Answer
For future googlers, Ding et al., 2017 published a review: "Comparison of Alternative Splicing Junction Detection Tools Using RNA-Seq Data".

Table 1 has the list of software in this review.
> 1 votes

---
Tags: rna-seq, software-recommendation, rna-splicing
---


# Question
Title: assigning different characters to present and absent genes from blast output using perl or python
I have file called pangenome.txt which has only one column containing the names of the genes

```
head pangenome.txt
    Ha1_00044161
    Ha11_00031248
    Ha16_00021664
    Ha8_00042762
    Ha4_00025807

```

I have another bunch of files in a tabular blast format (330 in total: individual1.blast.txt,individual2.blast.txt, ...,individual330.blast.txt ) which are the outputs of Diamond protein aligner. Below is part of individual1.blast.txt file

```
head blast.individual1.txt
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00021664   77.4    31  7   0   94  2   49  79  7.3e-159.3
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha8_00042762    74.2    31  8   0   94  2   1   31  9.6e-158.9
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha4_00025807    77.4    31  7   0   94  2   1   31  1.3e-058.5
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00022939   77.4    31  7   0   94  2   295 325 2.1e-057.8
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha12_00032339   71.0    31  9   0   94  2   2   32  2.8e-057.4
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00021553   71.0

```

I want to know which genes from the pangemome.txt is present in my individuals file and which are absent. for the example above this is my desired output:

```
head output.individual1.txt
    Ha1_00044161   0
    Ha11_00031248  0
    Ha16_00021664  1
    Ha8_00042762   1
    Ha4_00025807   1

```

My desired out put is a file with 2 columns (tabular format), the first column the list of the all genes form pangemome.txt pan, and the second column is the presence and absence of these genes in the individuals files if the gene is present present in the blast.individual1.txt, I want to assign 1 and if it is absent I want to assign 0. Can anyone help me to do this? Thanks a lot

# Answer
> 3 votes
I wouldn't do this in Python, myself. This is a very simple text parsing problem and the standard \*nix tools will be able to do it very easily. For example with `awk`:

```
$ awk 'NR==FNR{a[$1]=0; next}
       {
        if($2 in a){
            a[$2]=1
        }
       }END{
        for(i in a){
            print i,a[i]
        }
       }' pangenome.txt blast.individual1.txt
Ha1_00044161 0
Ha4_00025807 1
Ha11_00031248 0
Ha16_00021664 1
Ha8_00042762 1

```

The above is just a one-liner, I expanded it for clarity. You can run it on as many files as you like:

```
awk 'NR==FNR{a[$1]=0; next}{if($2 in a){a[$2]=1}}END{for(i in a){print i,a[i]}}' pangenome.txt blast.individual*.txt

```

For many large files, the fastest approach would probably be to just grep the names (unlike the awk, this will also match if the target gene is in the 1st field, but that's probably not an issue for you):

```
$ grep -wof pangenome.txt blast.individual1.txt.orig
Ha16_00021664
Ha8_00042762
Ha4_00025807

```

That will give you the list of genes found in the target blast files. To get the 1 and 0, you can do:

```
$ sort <(grep -owf pangenome.txt blast.individual1.txt) pangenome.txt |
    uniq -c | awk '{print $2,$1-1}'
Ha1_00044161 0
Ha11_00031248 0
Ha16_00021664 1
Ha4_00025807 1
Ha8_00042762 1

```

For multiple files, just run it on all of them:

```
$ sort <(grep -owf pangenome.txt blast.individual*.txt) pangenome.txt |
    uniq -c | awk '{print $2,$1-1}'

```

# Answer
> 2 votes
Using `python` it should be relatively trivial, using `sets`. The following snippet works with your pangenome file and one blast file; adding more blast files should be trivial, making use of the `union` method of `python`'s `sets`.

```
pangenome = {x.rstrip() for x in open('pangenome.txt')}
blast = {x.rstrip().split()[1] for x in open('blast.txt')}
for x in pangenome:
    if x in blast:
        print('%s\t1' % x)
    else:
        print('%s\t0' % x)

```

# Answer
> 2 votes
Here is a way to do it using a python dictionary:

```
#!/usr/bin/env python3

import sys

gene_filename = sys.argv[1]
blast_filename = sys.argv[2]

# Set to 0 all genes in the gene file
genes = {}
with open(gene_filename, "r") as gene_file:
    for line in gene_file:
        genes[line.strip()] = 0

# Set to 1 those genes found in the blast file
with open(blast_filename, "r") as blast_file:
    for line in blast_file:
        gene = line.strip().split()[1]
        if gene in genes:
            genes[gene] = 1

# Print the results
for (gene, presence) in genes.items():
    print(gene, presence, sep="\t")

exit(0)

```

Executing it:

```
$ ./find_genes_in_blast.py pangenome.txt blast.individual1.txt
Ha1_00044161    0
Ha11_00031248   0
Ha16_00021664   1
Ha8_00042762    1
Ha4_00025807    1

```

---
Tags: python, blast, perl
---


# Question
Title: bed file merge book-end features only if same name in column 4
I have a bed file where some entries have the same name in column 4 but are book-ended in coordinates, e.g. see below:

```
chr1    0       10000   LOC101929784::chr1:13201-13800::E1
chr1    10000   10600   LOC101929784::chr1:13201-13800::XYZ
chr1    10600   10800   LOC101929784::chr1:13201-13800::E1
chr1    10800   12000   DDX11L1::chr1:11868-14409.E1
chr1    12000   12200   DDX11L1::chr1:11868-14409.E1
chr1    12200   13000   DDX11L1::chr1:11868-14409.E1
chr1    13000   13200   DDX11L1::chr1:11868-14409.E1
chr1    13200   15800   WASH7P::chr1:14403-29570.E1
chr1    15800   16000   WASH7P::chr1:14403-29570.E1
chr1    16000   18200   WASH7P::chr1:14403-29570.E1
chr1    18200   18400   WASH7P::chr1:14403-29570.E1
chr1    18400   23000   WASH7P::chr1:14403-29570.E1
chr1    23000   23400   WASH7P::chr1:14403-29570.E1
chr1    23400   24000   WASH7P::chr1:14403-29570.E1
chr1    24000   24200   WASH7P::chr1:14403-29570.E1
chr1    24200   25200   WASH7P::chr1:14403-29570.E1
chr1    25200   25800   WASH7P::chr1:14403-29570.E1
chr1    25800   26200   WASH7P::chr1:14403-29570.E1
chr1    26200   26400   WASH7P::chr1:14403-29570.E1
chr1    26400   27000   WASH7P::chr1:14403-29570.E1
chr1    27000   27600   WASH7P::chr1:14403-29570.E1
chr1    27600   27800   WASH7P::chr1:14403-29570.E1
chr1    27800   29600   WASH7P::chr1:14403-29570.E1
chr1    29600   30800   MIR1302-2HG::chr1:29553-31109.E1
chr1    30800   31400   MIR1302-2HG::chr1:29553-31109.E1

```

I would like to apply an operation that seems to me is a mix of `bedtools merge -d 0` and `bedtools groupby`, but I can't figure out what's the combination.

For example, the 4 entries with `DDX11L1::chr1:11868-14409.E1` I would like merged as they are book-ended, same with the `WASH7P::chr1:14403-29570.E1` entries and the `MIR1302-2HG::chr1:29553-31109.E1` entries. There are two `LOC101929784::chr1:13201-13800::E1` entries, but I don't want to merge them because they are not book-ended, since they are separated by `LOC101929784::chr1:13201-13800::XYZ`.

Any ideas?

# Answer
```
awk 'BEGIN{OFS="\t"}{if($1==lchrom && $4==lname && $2 == lend) {lend = $3}else{if(lchrom) {print lchrom, lstart, lend, lname;}; lchrom=$1; lstart=$2; lend=$3; lname=$4}}END{print lchrom, lstart, lend, lname}' foo.bed > new.bed

```

That produces:

```
chr1    0   10000   LOC101929784::chr1:13201-13800::E1
chr1    10000   10600   LOC101929784::chr1:13201-13800::XYZ
chr1    10600   10800   LOC101929784::chr1:13201-13800::E1
chr1    10800   13200   DDX11L1::chr1:11868-14409.E1
chr1    13200   29600   WASH7P::chr1:14403-29570.E1
chr1    29600   31400   MIR1302-2HG::chr1:29553-31109.E1

```

This explicitly checks for bookends and nothing else and assumes that everything is sorted.
> 5 votes

---
Tags: bed, bedtools
---


# Question
Title: Normalize ATAC-seq/Dnase-seq sequencing reads coverage signals over estimated background
I'm trying to normalize the coverage signals of ATAC-seq reads against its own background using normal distribution, described in this paper

It says:

> Finally, all open chromatin coverage measurements were normalized by standardization to the mean and standard deviation of coverage over a set of 25,000 randomly selected background regions. To select background regions, the set of peak open regions were widened to 20,000 bp, reduced, and subtracted from the genome assembly. Thereafter, 25,000 random positions were selected and widened to reflect the distribution of widths in the set of open peaks. Coverage within these background regions was then calculated, and regions with zero coverage were discarded (~5%). The distribution of counts within background regions approximated a log-normal distribution. Mean and standard deviation of these background regions was calculated and used to transform the coverage measurements for the entire genome.

Does anyone know any existing codes doing that? I know Anshul's pipeline does that but can't find the specific script for that step.

# Answer
> 2 votes
Supplementary File  and  contain the code you're looking for.

---
Tags: normalization, atac-seq
---


# Question
Title: Where can I find gene expression data on one of the cell lines in NCI-60?
I'm trying to retrieve the most expressed genes from a cell line from the NCI-60, lets say cell line NCI-H23, but I'm not sure where to find the gene expression data. My main goal is to create a network using that data and iRefIndex to be loaded into Cytoscape. Does anyone know where I can find the correct data or guide me in the right direction?

# Answer
> 2 votes
CellMiner portal has plenty of data about NCI-60 cell lines, inclusive microarray gene-expression.

Also, NCI-60 is a subset of the ~1000 cell lines covered by the GDSC project for which you have multiple types of measurements (e.g. copy-number, mutation, gene-expression).

# Answer
> 2 votes
Cosmic Cell Lines Project provides gene expression data for 1000+ cell lines (including NCI-60 panel together with NCI-H23 cell line). For each gene in a cell line, you can get Z-score and over/under-expression status calculated across all tumor samples diploid for that gene.

---
Tags: gene, drugs, cell-line
---


title: How to read .bwt file .sa file?

text:
# Question
Title: How to read .bwt file .sa file?
In BWA project in bwt.c file I found bwt\_dump\_bwt method and bwt\_dump\_sa method. I want to use them to read .bwt file and .sa file. I wrote the following program:

```
int main(void)
{
  char *tmp, *prefix;
    bwt_t *bwt;
    prefix = bwa_idx_infer_prefix("hg19bwaidx");

    tmp = calloc(strlen(prefix) + 5, 1);
    strcat(strcpy(tmp, prefix), ".bwt");
    bwt = bwt_restore_bwt(tmp);

  return 0;
}

```

I got the following error:

```
Undefined symbols for architecture x86_64:
  "_bwa_idx_infer_prefix", referenced from:
      _main in readbwt-3022ec.o
  "_bwt_restore_bwt", referenced from:
      _main in readbwt-3022ec.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)

```

Any help ?

Thanks.

# Answer
> 3 votes
In the bwa source code, file bwt.{c,h} implement basic operations on .bwt and .sa files. They depend on utils.{c,h}, but you can create a fake utils.h to cut this dependency as is shown below. Then you only need three files: kvec.h, bwt.h and bwt.c to access .bwt and .sa.

The following example demonstrates how to find super-maximal exact matches (SMEMs). Note that internally bwa concatenates all database sequences into one long string. If you use bwt.{c,h} only, you get coordinates on this long string but not on each chromosome. You have to call APIs in bntseq.{c,h}, which are fairly independent, too. You can go on to implement a mini bwa-mem with this example, but you will need all the bwa-mem source code, except the main function.

## Compiling and running

```
# compile the example
git clone https://github.com/lh3/bwa
mkdir -p tiny-smem && cd tiny-smem
cp ../bwa/{bwt.c,bwt.h,kvec.h} .              # you need at least these three files
### then create utils.h and test.c shown below ###
gcc -g -O2 -Wall bwt.c test.c -o tiny-smem    # compile tiny-smem

# create a bwa index if you don't have it. We use chr20 from GRCh38 as an example:
(cd ../bwa && make)                           # compile bwa for indexing
../bwa/bwa index -p chr20 path/to/human/chr20.fa

# run the example
./tiny-smem chr20.bwt chr20.sa TTTTCTTCTTTTGTAATATAAGCATTgGTTGTTCTAAATTTCTCTCTTGGCACTGCTTTA

```

The last command line above outputs:

```
0       0       26      1       -125734
0       14      27      4       -3657204        +16021887       +11802682       +54692983
0       15      29      1       +8344885
0       17      31      1       -39545379
0       19      37      1       +55060482
0       25      39      2       +1198672        -26409692
0       27      60      1       -125700

```

where column 2 and 3 give the coordinate of SMEM on the query, col4 the number of exact hits and the following columns are the list of hits in the concatenated coordinate.

## File: utils.h

```
#ifndef FAKE_UTILS_H
#define FAKE_UTILS_H

#define xassert(cond, msg) assert(cond)
#define xopen(fn, mode) fopen((fn), (mode))
#define err_fclose(fp) fclose(fp)
#define err_fread_noeof(ptr, size, n, fp) fread((ptr), (size), (n), (fp))
#define err_fwrite(ptr, size, n, fp) fwrite((ptr), (size), (n), (fp))
#define err_fflush(fp) fflush(fp)
#define err_fseek(fp, off, whence) fseek((fp), (off), (whence))
#define err_ftell(fp) ftell(fp)

#endif // ~FAKE_UTILS_H

```

## File: test.c

```
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <stdio.h>
#include "bwt.h"

int main(int argc, char *argv[]) {
  bwt_t *bwt;
  int j, max_intv = 100; // ignore SA intervals larger than this value
  bwtintv_v r = {0,0,0};

  if (argc < 4) { fprintf(stderr, "Usage: ./test <.bwt> <.sa> <seq>\n"); return 1; }
  if ((bwt = bwt_restore_bwt(argv[1])) == NULL) abort(); // read .bwt
  bwt_restore_sa(argv[2], bwt);                          // read .sa
  for (j = 3; j < argc; ++j) {
    int i, k, x, len = strlen(argv[j]);
    uint8_t *seq = (uint8_t*)malloc(len);
    for (i = 0; i < len; ++i) { // convert A/C/G/T/N => 0/1/2/3/4
      int c = tolower(argv[j][i]);
      if (c == 'a') seq[i] = 0;
      else if (c == 'c') seq[i] = 1;
      else if (c == 'g') seq[i] = 2;
      else if (c == 't') seq[i] = 3;
      else seq[i] = 4;
    }
    x = 0, r.n = 0;
    while (x < len) {
      x = bwt_smem1(bwt, len, seq, x, 1, &r, 0); // find SMEMs
      for (k = 0; k < r.n; ++k) { // traverse each SMEM
        bwtintv_t *p = &r.a[k]; // p->info encodes start and end on the query
        printf("%d\t%d\t%d\t%d", j-3, (int)(p->info>>32), (int)p->info, (int)p->x[2]);
        if (r.a[k].x[2] <= max_intv) { // don't retrieve coordinates for huge intervals
          int l_mem = (int)p->info - (int)(p->info>>32);
          for (i = 0; i < p->x[2]; ++i) { // traverse each hit
            uint64_t y = bwt_sa(bwt, p->x[0] + i); // get suffix array coordinate
            if (y < bwt->seq_len>>1) printf("\t+%ld", (long)y);
            else printf("\t-%ld", (long)(bwt->seq_len - y - l_mem));
          }
        }
        putchar('\n');
      }
    }
    free(seq);
  }
  free(r.a);
  bwt_destroy(bwt);
  return 0;
}

```

# Answer
> 4 votes
You need to include `bwa.o` in your call to `clang`, since it holds the definitions of those functions. This will cause a good number of object files, so you might as well just `clang -o foo *.o foo.c -lz -lpthread -lm -lrt` and get it over with. Since this was quite obviously never meant to be done, you'll have to remove at least one `main()` function from `bwa.c`.

Note that further questions regarding this will be considered off-topic and closed, since it's purely about compiling something that you want to write yourself.

---
Tags: sequence-alignment, bwa
---


# Question
Title: Combine Fastq by writeFastq is not working properly
Currently I have 4 fastq file for each patient as WES data. For WES pipeline, I want to combine those fastq files. I am using writeFastq from ShortRead package in Bioconductor.But I think this function is not working properly.

The reason of my doubt is the the size the next step is quality check using FastQC tool. This tool is working on individual file but not on combined file and shows the error that ID line doesn't start with "@" . So I want to know what are the other effective tool for combining multiple fastq files without any loss of information.

Let me know if I am doing something wrong. Thanks

**\[EDIT-1\]** : The codes for combining multiple fastq files into one is as follows:

```
reference_fastq = file.path(paste(getwd(),"refs",sep = "/", collapse = " "), fsep = ".fasta")
fastqDir <- file.path(dataDir, fsep = ".fastq")
outfile_temp <- paste(dirname(dataDir),basename(fastqDir),    paste(basename(fastqDir),"combined","fastq", sep = "."),sep = "/", collapse = " ")
fls <- list.files(fastqDir, "fastq$", full=TRUE)
for (fl in fls) {
  fq <- readFastq(fl)
  writeFastq(fq,outfile_temp, mode="a")
}

```

First four lines of my fastq file (one of 4 indivual file) are as follows :

```
@D00239:185:C9KELANXX:2:2201:1088:2107 1:N:0:TCCTGAGC
GTACAAGAAGCATAACATCAGCATGTGTTTTGGTGAGAACCTTAGGAAGTTTGCAATCATGACAGAAGGTGAAGG
+
/<BBBFFFBFFFFFFFFFFFBF<FFBF<B<FFFFFFFFFFB/B/<////<<<FF<B<FF<FFFFFFFFF/FFFFF

```

First four lines of combines fastq file are as follows:

```
@D00239:185:C9KELANXX:2:2201:1088:2107 1:N:0:TCCTGAGC
GTACAAGAAGCATAACATCAGCATGTGTTTTGGTGAGAACCTTAGGAAGTTTGCAATCATGACAGAAGGTGAAGG
+
/<BBBFFFBFFFFFFFFFFFBF<FFBF<B<FFFFFFFFFFB/B/<////<<<FF<B<FF<FFFFFFFFF/FFFFF

```

It is pretty clear that combined fastq file is also starting wiht proper identifier but why FastQC tool is saying that identifier '@' is missing. Here is the screen shot of the error message:

Thanks.

# Answer
You can drastically simplify your effort because FASTQ files are simple text files; therefore, standard UNIX test file tools work:

```
cat refs/*.fastq > combined.fastq

```

And this even works with `.fastq.gz` files — no unzipping required.

---

Apart from this, your R code contains several errors. For instance, your variable names change throughout, and you’re using the `fsep` parameter wrongly. In fact it’s unclear to me what the code is supposed to do. At an approximation (and untested), the following should work:

```
fastq_dir = file.path(getwd(), "refs")
fastq_files = dir(fastq_dir, pattern = "\\.fastq$", full.names = TRUE)
outfile = file.path(dirname(fastq_dir), "combined.fastq")

for (file in fastq_files) {
    fq = readFastq(file)
    writeFastq(fq, outfile, mode = "a")
}

```

This is equivalent to running the shell command shown above, just in a more convoluted way (parsing and re-generating the FASTQ files).
> 4 votes

---
Tags: r, bioconductor, sequence-analysis
---


title: Combine multiple species trees

text:
# Question
Title: Combine multiple species trees
I have five species trees run using starBeast and Beauti. Here is a link to the paper discussing both, http://journals.plos.org/ploscompbiol/article/file?type=supplementary&id=info:doi/10.1371/journal.pcbi.1003537.s006 Each of these trees has the same individuals, and each locus used (10 loci per tree) has individual models set to it. All the trees used a relaxed log clock normal. I would like to be able to combine these trees into one. I have attempted to find a way to do this within Densitree and starBeast, but to no avail. When all loci are run with all trees, the program simply dies, leading me to believe that running the trees together will all parameters is computationally impossible for my hardware. Is there a way to combine these trees?

# Answer
> 2 votes
Although I've never tried it myself, I think what you need can be achieved through the LogCombiner program from the BEAST package -- assuming the species trees files are comparable (ran with same parameters, same species names, etc.) I am also assuming that you meant 'five species tree files', each file with the posterior sample of species trees.

However, I do believe you should try to run all loci within the same analysis, or at least do some prior check to see if the distributions between data sets are comparable. LogCombiner assumes **independent runs from the same data**.

I would also recommend taking a look (and following) the BEAST google group (which includes starBEAST discussions).

---
Tags: gene, phylogenetics, phylogeny
---


# Question
Title: How do programs like STRUCTURE determine populations from alignments?
If I were to run BAM/SAM files through a program like STRUCTURE, or ADMIXTURE described here how would they determine which individuals belong to what groups?

What part of the files would these programs use to determine which individuals belonged to different populations? I believe they do this via determining allele frequencies for each individual and then putting them into a population based on allele frequency, is the only criteria that they use comparison of allele frequencies between individuals?

# Answer
*structure* works with SNP/genotype information, so won't work directly off BAM/SAM files. You would need to convert to VCF (e.g. using mpileup/vcftools), then process using something like PLINK to convert it into a format that STRUCTURE recognises.

Here's the description I have about *structure* from my PhD thesis. It might be a bit dated now, but I expect most of it is still applicable:

A clustering program called *structure* was created by Falush et al. (2003) that is used to determine population structure using a set of unlinked genotypes. The expected use case for *structure* is in the analysis of population structure, but it can also be useful for determining the utility of a given set of markers for categorising (or quantifying) a given trait.

The main benefit of using *structure* over other methods is that it is able to derive estimated ancestry coefficients from the data, rather than a discrete yes/no for each individual. This allows for the use of cutoff values for the assignment of individuals (useful for tweaking false positive and false negative rates in diagnostic tests), and a quantifiable value where a classification is not useful or does not make sense (as in the case of a continuous trait with high heritability such as height).

The program begins by assigning each individual to a given population. This can be either given in the input data as a ‘popinfo’ flag, or determined randomly by structure. After this, the program constructs a probability distribution for the data, considering the population assignments given, and uses that to estimate population allele frequencies. The program then constructs another probability distribution for these data, and uses that to estimate the population that each individual is likely to have originated from. This process is repeated many times (such a process is known as a random walk), and the underlying theory suggests that the progression of population assignments through the random walk process can be used as an approximation of the actual populations of origin for each individual.

The documentation included with *structure* describes the general algorithm for each iteration (of maybe hundreds of thousands of iterations) as follows:

1. Sample (population allele frequencies for iteration m) from the probability of (allele frequency) given (genotypes of sampled individuals) and (population of origin of individuals for iteration m − 1).
2. Sample (population of origin of individuals for iteration m) from the probability of (population of origin) given (genotypes of sampled individuals) and (allele frequencies for iteration m).

In more simple terms, the steps are similar to the following statements:

1. these populations probably have *this* allele frequency distribution therefore...
2. these individuals probably came from *this* population

The *structure* program is able to accommodate admixture into its model by considering the population of origin as a probability, rather than a single discrete value.

The file that is output from *structure* is a text file with individual identifiers and predicted ancestry coefficients (Q values). While the program does provide means to display these data graphically, the visual representation of data (including sorting) is better controlled using an external program. The most commonly used program is Noah Rosenberg’s distruct (Rosenberg, 2004), a program that allows colour customisation and group re-ordering. For the purpose of this thesis, a custom R script, `snpchip2structure.R`, has been designed which also allows selection of different sorting methods, a scatter plot for two populations, error bars, and a few additional features (see figure 1.7). More details about this script can be found in the programs directory of the supplementary CD for this thesis.
> 2 votes

---
Tags: sam, phylogenetics
---


# Question
Title: Alignment QC differences between HISAT2 and Qualimap
I used hisat2 for aligning reads to to the genome. I have an alignment summary for sample1 as follows:

```
89038751 reads; of these:
  89038751 (100.00%) were paired; of these:
    5641019 (6.34%) aligned concordantly 0 times
    67665552 (76.00%) aligned concordantly exactly 1 time
    15732180 (17.67%) aligned concordantly >1 times
    ----
    5641019 pairs aligned concordantly 0 times; of these:
      651415 (11.55%) aligned discordantly 1 time
    ----
    4989604 pairs aligned 0 times concordantly or discordantly; of these:
      9979208 mates make up the pairs; of these:
        5887822 (59.00%) aligned 0 times
        2656902 (26.62%) aligned exactly 1 time
        1434484 (14.37%) aligned >1 times
96.69% overall alignment rate

```

To check the alignment quality I used Qualimap. In the report I see the following:

```
Reads alignment

Number of mapped reads: 172,189,680
Total number of alignments: 244,302,519
Number of secondary alignments: 72,112,839
Number of non-unique alignments:    22,973,098
Aligned to genes:   2,081
Ambiguous alignments:   173
No feature assigned:    178,975
Not aligned:    5,887,822

Reads genomic origin

Exonic: 2,081 / 1.15%
Intronic:   6,866 / 3.79%
Intergenic: 172,109 / 95.06%
Intronic/intergenic overlapping exon:   302 / 0.17%

Transcript coverage profile

5' bias:    0.69
3' bias:    2.85
5'-3' bias: 0.39

```

Can anyone tell me why the numbers look different? Did I go wrong somewhere? Help me in understanding the results.

# Answer
> 2 votes
89038751 pairs times 2 is 178 million reads. That times the approximately 97% alignment rate is around 172 million reads, which is what qualimap is reporting.

A "non-unique" alignment is one that can map to multiple places (either equally well, or where all of the alignments are of similar enough quality). Such alignments have a single "primary" alignment and one or more "secondary" alignments, which is where the "secondary" metrics come from.

For the qualimap output, there are additional metrics describing how reads overlap genes (or not). Of note are the "ambiguous" alignments, which are those overlapping multiple genes. "No feature" metrics, then, refer to alignments not overlapping any gene (see the qualimap documentation for further details).

My guess is that you have such a high "intergenic" rate because you've either used the wrong GTF file or have told qualimap the wrong strand setting. For unusual organisms one expects many reads to arise from unannotated genes. This is not the case for common model organisms, such as mice and humans, where there are very few remaining commonly expressed unannotated genes.

---
Tags: rna-seq, sam, sequence-alignment, rna
---


# Question
Title: Filter bam using SNP list in bed format with minimum mapping quality and base quality
I have a bam file and a bed file that defines a list of SNPs. I would like to filter the bam file to contain only those reads with a minimum mapping quality that overlap at least one SNP with a minimum base quality.

Samtools seems to almost solve this problem, but not quite:

`samtools depth -b snp_list -q 20 -Q 30 example.bam` allows me to count the number of alignments meeting my desired criteria, but not to produce a bam with these alignments.

`samtools view -L snp_bed -q 30 example.bam` allows me to filter the alignments overlapping any SNP with minimum mapping quality, but not with a minimum base quality at the SNP sites. There is no `-Q` option for samtools view, and the `-q` option is for mapping quality

Can this be achieved in samtools? Or is there another preexisting tool that can?

# Answer
Thanks Devon ;-) using samjdk : http://lindenb.github.io/jvarkit/SamJdk.html

a bed file 'input.bed'

and the following snippet:

```
private IntervalTreeMap<Boolean> treeMap=null;
private final int BASE_QUALITY= 20;
@Override
public Object apply(final SAMRecord record) {
if(record.getReadUnmappedFlag() || record.getCigar()==null) return null;
if(treeMap==null) {
    treeMap = new  IntervalTreeMap<Boolean>();
    final String bedStr;
    try {
        bedStr = IOUtil.slurp(new java.io.File("input.bed"));
        }
    catch(java.io.IOException err)
        {
        throw new RuntimeIOException(err);
        }
    Arrays.stream(bedStr.split("[\n]")).
        filter(L->!(L.startsWith("track") || L.isEmpty() || L.startsWith("browser"))).
        map(L->L.split("\t")).
        map(A->new Interval(A[0],Integer.parseInt(A[1])+1,Integer.parseInt(A[2]))).
        forEach(I->treeMap.put(I,true));
    }

byte quals[]= record.getBaseQualities();
int refpos= record.getStart();
for(CigarElement ce:record.getCigar()) {
    CigarOperator op= ce.getOperator();
    if(op.consumesReferenceBases())
        {
        if(op.consumesReadBases()) {
            for(int i=0;i< ce.getLength();i++)
                {
                Interval t = new Interval(record.getContig(),refpos+i,refpos+i);
                if(!treeMap.containsOverlapping(t)) continue;
                int readpos = record.getReadPositionAtReferencePosition(refpos+i);
                if(readpos==0 ) continue;
                readpos--;
                if(readpos>=quals.length)  continue;
                if(quals[readpos]>BASE_QUALITY) return true;
                }
            }

        refpos+=ce.getLength();
        }

    }
return null;
}

```

usage:

```
samtools view -bu -L input.bed input.bam |\
   java -jar dist/samjdk.jar --body -f snippet.code

```
> 2 votes

# Answer
You're not going to get samtools to do that, though I wouldn't be surprised if something in jvarkit can do it. Baring that, here's a bit of python.

```
import pysam

bed = open("foo.bed")
bam = pysam.AlignmentFile("foo.bam")
of = pysam.AlignmentFile("filtered.bam", mode="wb", template=bam)
seen = set()  # memoize reads seen to not double print them
for line in bed:
    cols = line.strip().split("\t")
    for col in bam.pileup(cols[0], int(cols[1]), int(cols[2])):
        for b in col.pileups:
            if b.alignment.mapping_quality < MIN_MAPQ:
                continue
            pos = b.query_position
            if not pos:
                continue
            if b.alignment.query_qualities[pos] < MIN_PHRED:
                continue
            if b.alignment.query_name not in seen:
                seen.add(b.alignment.query_name)
                of.write(b.alignment)
of.close()
bam.close()
bed.close()

```

Or something along those lines, I haven't debugged that.
> 1 votes

---
Tags: sam, samtools, bed
---


# Question
Title: ENA minion 2D data doesn't actually have 2D data
I am looking at some yeast ONT data from De novo yeast genome assemblies from MinION, PacBio and MiSeq platforms, Giordano et al 2017.

I noticed that the Minion fastq files have both a forward and reverse file. The methods paper says that these reads are from the SQK-MAP006 and SQK-MAP007 protocols which both appear to be 2D. Are these just the 2D forward and reverse files?

Still being unsure, I checked if this was true by counting the number of characters in each read.

```
>> zcat ERR1883385_1.fastq.gz | head -1
@ERR1883385.1 dc5b9de2-ba00-47ab-b698-8ad470704c0a_Basecall_2D_000_template EVO107591_Sanger_FAA52213_CBS432_Phenol_3a_4517_1_ch251_file44_strand/1

>> zcat ERR1883385_2.fastq.gz | head -1
@ERR1883385.1 dc5b9de2-ba00-47ab-b698-8ad470704c0a_Basecall_2D_000_template EVO107591_Sanger_FAA52213_CBS432_Phenol_3a_4517_1_ch251_file44_strand/2

>> zcat ERR1883385_1.fastq.gz | head -2 | awk '{ print length($0); }' -
147
1675

>> zcat ERR1883385_2.fastq.gz | head -2 | awk '{ print length($0); }' -
147
1693

```

As you can see, the read is the same but the forward and reverse strands are different lengths.

Is there some way to just download the `poretools` equivalent to `2D` from ENA without having to download the fast5 files and extract them myself?

# Answer
> 2 votes
Probably not. ENA doesn't have any ability to specifically flag a file as 2D, and the nanopore file formats have changed a lot, so you pretty much need to rely on the submitter to separate out the reads in a meaningful way.

Also, there wasn't really any consistent convention for what ONT used in their read headers. Looking at some of my old files they had a plain sequence ID (i.e. without anything suggesting temp/comp) in different template/complement folders in the FAST5 files for the components of a 2D read and the consensus sequence, but later changed it to \_2D\_template and \_2D\_complement as in the example you have (possibly with an intermediate \_1D\_template and \_1D\_complement format). Parsing/understanding those FAST5 files can be difficult, which is why poretools was created in the first place.

It's still a bit of a mess, with the 1D$^2$ reads ending up in a Frankenstein-like state that's partially 1D and partially 2D.

---
Tags: minion
---


title: removing indels in plink

text:
# Question
Title: removing indels in plink
I have data in Plink format. I removed Indels using --snps-only command in Plink, but still some snps are as below:

```
rs58119544;rs5746975    0   16915280    T   A

```

I wonder, how i can get rid of these snps? the correct should be single snps without this (;) sign. Any idea to fix it please?

# Answer
Use *grep* to get SNPs that have `";"` in the name:

```
>$ grep ";" myMap.txt | cut -f1 -d" " > myBadSNPs.txt

```

Then use plink with --exclude :

```
plink --file myFile
--exclude myBadSNPs.txt
--out myFileFiltered

```
> 1 votes

---
Tags: plink
---


# Question
Title: Optimal design for an RNA-seq experiment with ERCC RNA Spike-In Control Mixes
I am planning an experiment on gene expression. I have 6 conditions and 3 replicates in each condition. I would like to use ERCC spike-in mixes for quality check: lower limit of detection, dynamic range and also check that the fold-changes are accurate. For example, by knowing that spike-in sequence xyz in condition 2 should be a 4-fold higher than in condition 1, does it actually get detected as a 4-fold change in my analysis? To assess the fold-changes, I need to use 2 different mixes. Only one mix should be used per sample. The manual suggests either adding mix 1 in a reference and mix 2 in all other samples or randomly assigning mixes 1 and 2 among samples. The latter is not very clear to me. Should I use the same mix for all replicates within a condition or should I randomise really among all the samples? Or can I randomise among conditions (for example 1, 3, 4 get mix 1 and 2, 5, 6 get mix 2)? So far, I am intending to do the latter, so using mix 1 in conditions 1, 3, 4 and mix 2 in conditions 2, 5, 6. This would also be cheaper because I can only buy equal volumes of mix 1 and 2. Do people have any advice on this?

# Answer
> 2 votes
> I have 6 conditions and 3 replicates in each condition.

If including more biological replicates has a low experimental cost, you should be using 6 replicates per condition (and fewer reads per condition, i.e. maintaining the total number of reads) in order to pick up the majority of differentially-expressed genes. For that, I'd recommend using 3 lanes of a flow cell, with 12 dual-indexed barcodes (with the same index on each end) for each lane. Doing dual indexing will protect a bit against index switching and make the results more robust.

> randomly assigning mixes 1 and 2 among samples. The latter is not very clear to me. Should I use the same mix for all replicates within a condition or should I randomise really among all the samples?

Randomly assign mixes among and within the conditions. Make sure that the mixes are equally represented (as much as possible) across the conditions, but that each condition has at least one of each mix. This can be done using something like the following in R:

```
numReplicates <- 7;
numEachMix <- ceiling(numReplicates / 2);
assignments <- sample(rep(c("mix1","mix2"),numEachMix))[1:numReplicates]

```

For working out the level of fold change this is reliably replicated (or "statistically significant"), you can split the samples into mix 1 and mix 2, then apply the same algorithms for ERCC differential expression as with the non-ERCC genes.

However, as ERCC is typically added in proportion to the amount of RNA (or, if you're lucky and have cell-sorted populations, possibly the number of cells), all it's really testing is the accuracy of sample prep.

---
Tags: rna-seq, rna, quality-control, spike-in, ercc
---


# Question
Title: Why is ribosomal RNA difficult to remove even with Poly(A) selection?
In this answer (actually in a comment), it is stated that:

> As you've noticed from your own analysis, the ribosomal genes have quite variable expression across cells. They're expressed everywhere, and quite **difficult to completely remove from a sample (even with polyA selection)**, with success of removal depending on things like the amount of mechanical damage that the RNA has been exposed to.

Which are the most relevant issues that make Poly(A) selection fail? Looking online, Poly(A) is suggested often but I didn't find an explanation for the different ways in which it can fail.

# Answer
We've found ribosomal RNA to be less of a problem with sequencing that depends on polyA, which suggests the issue might be in the library preparation, rather than the selection.

Many polyA RNA library preparation methods involve amplification, rather than selection, which means that existing transcripts that are present in very high abundance (such as rRNA) will still show up in the mixture.

But polyA-tailed RNA isn't the only RNA present in a sample. For a comprehensive look at RNA, the whole-RNA sample preparation methods can't use *any* form of selection to identify target molecules. In this case, it's typically recommended to use a form of ribosomal RNA removal, where a probe to ribosomal RNA is attached to magnetic beads (or a column, or similar), and used to fish out out ribosomal RNA, leaving the remainder behind. Unfortunately those probes don't work perfectly, either due to the ribosomal RNA being damaged at the probe site, or the probe not matching the ribosomal RNA molecule, or the RNA just not getting into the right physical state to bind with the probe. All these things mean incomplete removal of ribosomal RNA, ending up with a situation where the highly-abundant ribosomal RNA still pokes through in the reads.
> 5 votes

---
Tags: rna-seq, scrnaseq, rna, ribosomal
---


title: Strange per sequence GC content results

text:
# Question
Title: Strange per sequence GC content results
I would be grateful if someone could take a quick look at these FASTQC results. This is rna-seq paired-end data. From the FASTQC manual, an unusual distribution seems to be suggestive of contamination and a shift in the curve is suggestive of a systematic bias. GC content distribution both pre-alignment and post-alignment are strange.

Samples are paired end, strand specific and % of mapped reads is above 95% for all the samples. There is no adapter content also.

What could be the problem?

# Answer
Although generally I recommend people to be wary about the warnings that FastQC gives (it tends to be overly paranoid), the GC content graphs here do look odd. It's good that you've had a look at it with another program to confirm the observation. They should have a distribution that is fairly close to normal, a bit like this:

Without more information about your sample, I can't tell whether or not this is expected from your input. This looks like fairly substantial contamination, but could also be caused by biased selection during sample preparation (e.g. low quality sample, low yield, PCR amplification bias), or just looking at a multi-organism sample rather than a single-organism sample.

As a contamination check, I'd recommend running your RNASeq reads through a kmer-based classification tool (e.g. OneCodex, Centrifuge, BBSketch, or Sourmash). This should tell you the most likely place that contaminant sequences came from.
> 3 votes

---
Tags: rna-seq, sequence-alignment, read-mapping, rna
---


# Question
Title: Get gene exons from UCSC browser by HGNC name
I have a set of gene names in HGNC standard and I want to get all their exons (in BED format) from the UCSC browser but it seems like the browser doesn't index genes by their HGNC name. How can I do that?

# Answer
> 2 votes
You could grab exons via Gencode v27 (hg38) GFF annotations and convert them to BED via BEDOPS `convert2bed`/`gff2bed`:

```
$ wget -qO- ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/gencode.v27.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "exon"' \
    | convert2bed -i gff - \
    > exons.bed

```

Then use `grep` to filter exons against your list of gene names:

```
$ grep -w -f hgnc.txt exons.bed > exons.filtered.bed

```

# Answer
> 2 votes
If you need a 3-column BED, you can just use standard Unix tools:

```
curl ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/gencode.v27.annotation.gtf.gz \
  | gzip -dc \
  | grep -wFf list.txt \
  | awk -vFS="\t" '$3=="exon"{print $1"\t"($4-1)"\t"$5"\t"$9}' > out-3.bed

```

If you need a 12-column BED:

```
# install k8 and paftools.js
curl -L http://bit.ly/k8-024 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` k8
git clone https://github.com/lh3/minimap2
# actual command line:
curl ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/gencode.v27.annotation.gtf.gz \
  | ./k8 minimap2/misc/paftools.js gff2bed - \
  | grep -wFf list.txt > out-12.bed

```

---
Tags: ucsc, genome-browser, exon
---


# Question
Title: How to create artificial or dummy phenotype data
I have whole genome sequence (human) data (SNP data), but don't have phenotype data ready for any traits yet. Hopefully, I will get real phenotype data soon. To run GWAS (genome-wide-association study) or any other analysis, i want to create artificial/dummy data for any trait, let's say human height. Could anyone please recommend what is best option? Does plink generate any dummy phenotype data? or any other programme?

# Answer
> 3 votes
Using wakefield package it is pretty easy, here is an example from GitHub page. This will create 500 dummy variables for id, race, age, etc. :

```
library(wakefield)

myDummyPheno <- r_data_frame(
  n = 500,
  id,
  race,
  age,
  sex,
  hour,
  iq,
  height,
  died
)

# result
myDummyPheno

# # A tibble: 500 x 8
#    ID    Race        Age Sex    Hour           IQ Height Died
#   <chr>  <fct>     <int> <fct>  <S3: times> <dbl>  <dbl> <lgl>
#  1 001   Bi-Racial    27 Male   00:00:00     93.0   71.0 F
#  2 002   White        28 Female 00:00:00    104     64.0 T
#  3 003   Black        28 Male   00:00:00    111     67.0 T
#  4 004   Hispanic     35 Male   00:00:00     90.0   67.0 T
#  5 005   Hispanic     21 Female 00:00:00     94.0   70.0 T
#  6 006   Black        35 Male   00:00:00     74.0   73.0 F
#  7 007   Black        27 Male   00:00:00    102     68.0 T
#  8 008   White        28 Female 00:00:00     95.0   70.0 F
#  9 009   White        30 Female 00:00:00     80.0   71.0 F
# 10 010   White        30 Female 00:00:00    114     66.0 T

```

---
Tags: r, plink
---


title: seqret: Warning: bad /protein_id value

text:
# Question
Title: seqret: Warning: bad /protein_id value
I am using the seqret tool from the Emboss suite to transform produce an EMBL annotation file from a GTF annotation file and the corresponding FASTA sequence:

```
seqret -sequence ecoli.fa -feature -fformat gff -fopenfile ecoli.gtf -osformat embl -auto -outseq ecoli.embl

```

It produces my desired EMBL file but it warns against a 'bad protein\_id':

```
...
Warning: bad /protein_id value '"CCE57732"'
Warning: bad /protein_id value '"CCE57733"'
Warning: bad /protein_id value '"CCE57734"'
Warning: bad /protein_id value '"CCE57735"'
Warning: bad /protein_id value '"CCE57754"'

```

What's wrong with these protein IDs? The organism I am working with is *Escherichia coli*.

Here are a few lines from my GTF file. The first line has a protein\_id.

```
Chromosome      ena     CDS     349     2808    .       +       0       gene_id "HUS2011_0001"; transcript_id "CCQ26880"; exon_number "1"; gene_name "thrA"; gene_source "ena"; gene_biotype "protein_coding"; transcript_name "thrA-1"; transcript_source "ena"; transcript_biotype "protein_coding"; protein_id "CCQ26880";
Chromosome      ena     start_codon     349     351     .       +       0       gene_id "HUS2011_0001"; transcript_id "CCQ26880"; exon_number "1"; gene_name "thrA"; gene_source "ena"; gene_biotype "protein_coding"; transcript_name "thrA-1"; transcript_source "ena"; transcript_biotype "protein_coding";
Chromosome      ena     stop_codon      2809    2811    .       +       0       gene_id "HUS2011_0001"; transcript_id "CCQ26880"; exon_number "1"; gene_name "thrA"; gene_source "ena"; gene_biotype "protein_coding"; transcript_name "thrA-1"; transcript_source "ena"; transcript_biotype "protein_coding";

```

# Answer
> 1 votes
You're supplying a GTF feature file, but telling it that the features are in GFF3 format. I think you need `-fformat gff2`. See https://www.ebi.ac.uk/Tools/sfc/emboss\_seqret/help/

---
Tags: file-formats, format-conversion, gtf, emboss, seqret
---


# Question
Title: When working with Bioconductor and the curatedOvarianData dataset, what are the unit values of the expression data?
I am wondering how to interpret values inside the gene expression levels contained within Bioconductor. For example, if we have the following commands via `R`, we get:

```
source("https://bioconductor.org/biocLite.R")
biocLite("curatedOvarianData")
library(curatedOvarianData)
data(TCGA_eset)
head(exprs(TCGA_eset)["CXCL12",])

TCGA.20.0987 TCGA.23.1031 TCGA.24.0979 TCGA.23.1117 TCGA.23.1021 TCGA.04.1337
    5.437923     6.307250     4.193748     9.071293     5.101532     6.697770

```

The above calls out gene expression levels for the gene `CXCL12`. However, I am wondering how to interpret the values `5.437923`, `6.307250`, etc.

I know they are expression levels, but do they come with units? Do larger values mean greater cancer possibility?

# Answer
> 2 votes
These values represent the expression of a particular gene in each sample. While expression of an immune system gene like `CXLC12` might be associated with the aggressiveness of a cancer, that's not going to be the case for the majority of genes.

Looking at the numbers, these look like they are probably log$\_2$-transformed values. In order to know for sure, plotting a density distribution of the dataset is a good idea:

```
plot(density(exprs(TCGA_eset)));

```

If the values are in the range of around 0..15, then it'll be log$\_2$-transformed values, whereas if they're in the range of 0..100000, it'll be "raw" expression values. The log-transformed values are preferred for data analysis because gene expression tends to be close to a normal distribution in log-space.

These numbers are best treated as relative counts (i.e. the most appropriate *unit* is "expression"), as variation in the preparation of the sample (or the sequencer, or the microarray) can influence the results that come out the other end. Converting them to something more meaningful (e.g. "transcripts per cell") requires normalising against a known reference sample.

# Answer
> 2 votes
The values are `rma` normalized expression values from affy microarrays. They are indeed log2 transformed (as @gringer already expected). You can find more info about the data here, under heading `TCGA_eset`.

To get more info about `rma`, check this out.

Larger values are not necessarily greater cancer possibility, I wish biology was that simple.

---
Tags: r, bioconductor, gene
---


title: Difference between 1D and 1D^2 data

text:
# Question
Title: Difference between 1D and 1D^2 data
I am working on ONT data. Initially, I have worked on data from 1D using Minimap2 aligner. I came across 1D^2 from ONT website. I was wondering how different is the data from two techniques? What is the difference between 1D and 2D? and do we get 2 fastq files (one per strand) for each sample?

TIA

I would greatly appreciate any help regarding my question.

# Answer
1D² data that is called by Albacore generates two overlapping sets of files. Here's what Devon O'Rourke mentioned on the ONT community forums and on Twitter:

> Following Albacore basecalls on a 1d2 library I get two sets of .fq files, summary stats, etc.: one for the 1D basecalling script, and one for the 1D2 script. I'd like to use the reads generated from these scripts in a genome assembly.
>
> After a little bit of `grep` searching, it seems like there are overlapping reads among the 1D and 1D2 directories - namely, those 1d reads that generated 1d2 reads!

Link here, for those who have ONT community access.

I've recommended generating a list of the 1D^2 reads, then excluding those reads from the 1D fastq files. I have a couple of scripts that can help for this:

```
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq

```

More information on this can also be found in this StackExchange answer.
> 2 votes

---
Tags: nanopore
---


# Question
Title: What are the differences between BNGL (BioNetGen Language) and SBML (System Biology Markup Langague) formats?
Per request from meta comment.

I am self-learning about whole cell modeling, specifically An introduction to whole-cell modeling and Fundamentals of Systems Biology: From Synthetic Circuits to Whole-cell Models. While my background is decades of programming, for cellular biology it is only skimming through a few books.

In the process of learning I soon discovered software/apps that use SBML (System Biology Markup language) but in referring back to "An introduction to whole-cell modeling" see that they use BioNetGen (GitHub) which has its own language BioNetGen language (BNGL)

What are the differences between BNGL (BioNetGen Language) and SBML (System Biology Markup Langague) formats?

# Answer
Both formats have formal grammars for describing similar concepts, however they are different and to my knowledge one is not a proper subset of the other. While one could do a comparison of a syntactic nature, a simple Google search of just both acronyms reveals BioNetGen 2.2: Advances in Rule-Based Modeling

In here can be found:

> 2.4 SBML-to-BNGL translation
> SBML is a widely-used model exchange format in systems biology (Hucka et al., 2003). Models encoded in SBML are ﬂat, i.e., their species do not have internal structure, which limits their utility for rule-based modeling. BioNetGen 2.2 includes an SBML-to-BNGL translator, called Atomizer (also available as a web tool at ratomizer.appspot.com), that can extract implicit molecular structure from ﬂat species (Tapia and Faeder, 2013). A full report on Atomizer and its application to the BioModels database (Li et al., 2010) is currently in preparation. However, Tapia and Faeder (2013) reported that an early version of the tool could recover implicit structure for about 60% of species in models within the database that contain ≥20 species. Thus, Atomizer makes available a large set of pre-existing models in a rule-based format, facilitating their visualization (Wenskovitch et al., 2014) and extension (Chylek et al., 2015).

While not directly related to the question this also helped in understanding more about BioNetGen

Wikipedia: Multi-state modeling of biomolecules

> Biological signaling systems often rely on complexes of biological macromolecules that can undergo several functionally significant modifications that are mutually compatible. Thus, they can exist in a very large number of functionally different states. Modeling such multi-state systems poses two problems: The problem of how to describe and specify a multi-state system (the "specification problem") and the problem of how to use a computer to simulate the progress of the system over time (the "computation problem"). To address the specification problem, modelers have in recent years moved away from explicit specification of all possible states, and towards rule-based formalisms that allow for implicit model specification, including the κ-calculus, BioNetGen, the Allosteric Network Compiler and others. To tackle the computation problem, they have turned to particle-based methods that have in many cases proved more computationally efficient than population-based methods based on ordinary differential equations, partial differential equations, or the Gillespie stochastic simulation algorithm. Given current computing technology, particle-based methods are sometimes the only possible option. Particle-based simulators further fall into two categories: Non-spatial simulators such as StochSim, DYNSTOC, RuleMonkey, and NFSim and spatial simulators, including Meredys, SRSim and MCell. Modelers can thus choose from a variety of tools; the best choice depending on the particular problem. Development of faster and more powerful methods is ongoing, promising the ability to simulate ever more complex signaling processes in the future.
> 3 votes

---
Tags: file-formats
---


# Question
Title: Link to UCSC gtexGene track description page for selected gene symbol
I have this link that comes up when I am in the browser window and click on the gene name, which leads me to description page about the gene (in this example it is BOK gene).

How can I generate that link? Below are some of my failed attempts. There is a UCSC FAQ page about the links, didn't help much, maybe I am reading it wrong.

Link that works, generated when clicked on gene on the website:

http://genome-euro.ucsc.edu/cgi-bin/hgc?hgsid=226968178\_C6XmUjsjtt3Mm9fOmePo6VqRN6TV&c=chr2&l=241667282&r=242920780&o=242498135&t=242513546&g=gtexGene&i=BOK

Which gives me this page:

My failed attempts to generated custom link:

http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr2%3A241667283-242920780&gtexGene&i=BOK

http://genome.ucsc.edu/cgi-bin/org=human&db=hg19&gtexGene&i=BOK

http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&i=BOK

http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&hgg\_gene=BOK

http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&i=BOK

http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&g=gtexGene&i=BOK

http://genome.ucsc.edu/cgi-bin/org=human&db=hg19&g=gtexGene&i=BOK

http://genome.ucsc.edu/cgi-bin/hgTracks?org=human&db=hg19&g=gtexGene&i=BOK

# Answer
> 2 votes
You never want to include the hgsid parameter in your url because that is the link to your current browsing session and all its specific contents. So if you send that link to someone else and then they start browsing around that link may not work correctly again in the future if you or the person you sent the link to modify your browsing session.

As to how to create links to hgc pages, the usually follow a somewhat general format: http://genome-euro.ucsc.edu/cgi-bin/hgc?db=hg19&c=chr2&o=242498135&t=242513546&g=gtexGene&i=BOK

You need the db and g parameters to specify the assembly build and track, and then you need to provide the position in the form of the c,o and t parameters. The last thing you want is the i paramter to specify the item name. The reason the i parameter alone is not enough is because we do not require uniqueness on item names (there may be many copies of an item across the genome for instance).

You can obtain the c, o, t and i parameters by selecting the chrom, chromStart, chromEnd and name fields from the gtexGene table:

```
$ mysql --user=genome --host=genome-mysql.soe.ucsc.edu -A -P 3306 -e "select chrom, chromStart, chromEnd, name from gtexGene limit 1" hg38
+-------+------------+----------+---------+
| chrom | chromStart | chromEnd | name    |
+-------+------------+----------+---------+
| chr1  |      11868 |    13052 | DDX11L1 |
+-------+------------+----------+---------+

```

---
Tags: gene, ucsc
---


# Question
Title: efficient counting of dinucleotides/trinucleotides on fastq reads?
What's an efficient way of counting of dinucleotides/trinucleotide pattern on fastq.gz file reads?

I know there are tools like seqtk that will be very efficient at reading through the .fastq.gz files, and tools like cutadapt that are good at finding long(ish) adapter patterns in the fastq.gz PE of files, but, what would be the best tool to count the number of reads where a di-/trinucleotide is present N-number of times in the read?

E.g.

```
countingtool --pattern CGC file_R1.fastq.gz

```

Result: column two is number of times the pattern is seen (exclusive numbers, no overlap)

```
CGC    0     51341234
CGC    1     13948130
CGC    2      1234344
CGC    3        51344
CGC    4         1343

```

# Answer
> 1 votes
Well, I'm sure you can find a more efficient approach, but here's a simplistic solution using standard UNIX tools:

```
zcat file.fastq.gz | grep CGC |
    awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
         sort -rnk2

```

I tested this on a 2.6G fastq file (exome data) and it took:

```
real    1m18.987s
user    1m41.392s
sys     0m10.684s

```

On a 66G WGS sequence file, it took:

```
real    34m6.771s
user    40m16.452s
sys     4m10.176s

```

The extra `grep CGC` isn't needed, you can do the match in awk, but it does speed things up a little bit. You can always get the reads with 0 matches by counting the number of reads and subtracting the ones with 1. Alternatively, remove the `grep`:

```
zcat file.fq.gz | awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
    sort -rnk2

```

---

Following the OP's suggestion, I also tested with `pigz` and, it turns out, pigz can fly!

```
time unpigz -c exome.fastq.gz | grep CGC |
    awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
         sort -rnk2

real    0m44.012s
user    1m19.544s
sys     0m14.900s

```

And for the 66G WGS:

```
real    17m29.935s
user    29m23.816s
sys     5m17.856s

```

---
Tags: fastq, benchmarking, seqtk, cutadapt
---


title: What does it mean to over sample?

text:
# Question
Title: What does it mean to over sample?
One of my studies (A) involve sequencing the microbiome. After selecting the variable region and the primers for the targeted region of the 16S sequence. The samples were sent to a platform.

For another study (B) we used the same strategy (same primers, region, machine...) but another platform. However we obtained from this platform a lot more sequences:
In study A the mean sequences per sample are 10k, while in the study B the mean sequences per sample are 80k

I consulted with others and it seems that the study A we under-sampled and in study B we over-sampled.

I can understand that by under-sampling we didn't obtain enough reads to be sure about the microbiome present in the sample. But what does it mean to over-sampling?

That expert I consulted suggested normalize the data to fewer sequences, discarding information we already have. I can understand that the sequencing has errors and can create chimeras but the resulting consensus sequences of the OTUs have found a match with databases of microorganisms. What or how should I take into account that I over-sampled?

# Answer
> 4 votes
As you sequence a library you are sampling from the pool of all possible oligonucleotides found in it. This is a stochastic process, meaning that you are more likely to sequence a piece of DNA that is more prevalent.

As you increase the number of reads that you sequence for a library, you decrease the chances of seeing a new read that you have not yet seen before. If you plot the number of unique reads per total reads sequenced, it will look something like the cumulative density function of the exponential distribution, where you have diminishing returns on unique sequences as you continue to sequence.

If you sequence a library and quickly stop seeing unique reads, we say that this library has **low complexity**. Therefore, you have "oversampled" it if you sequence too deeply. I am guessing that your metagenomics library actually has a lot fewer species/unique sequences than they anticipated, and that you sequenced deeply anticipating more.

Edit: In response to your comments, the bigger problem in 16S analyses is not sequencing too deeply, it is over-amplifying your samples such that they no longer represent what is actually in your gut microbiome samples. In this case if you used read normalization you would not accurately represent what was in the library that you sequenced since the amplified molecules (16S) are so small that you only have len(16S molecule)-(read len) + 1 possible reads in each direction. You might consider removing PCR duplicates with GATK but this will still have the same effect of not representing what was in the sequencing pool. It is very possible that you have many PCR duplicates after 16S amplification. I generally don't trust amplicon sequencing for quantitive analyses, and don't really trust for low coverage "is-this-species-there-or-not" either, especially when dealing with the 'dark matter' of uncharacterized species. You might be well off just asking "is this species present in high coverage at this time point?"

# Answer
> 1 votes
I don't think there's anything inherently wrong with sequencing to high depth (though naive approaches to dealing with sequencing errors such as singleton removal may fail if the depth is so high that the same errors may be encountered repeatedly by chance). The bigger issue is that they are over sampled *relative* to the undersampled microbiomes. This does complicate comparisons, particularly for things like OTU richness. The obvious solution is to down-sample the oversampled microbiomes so you have the same number of reads for each. This is probably unnecessary if your analyzing the oversampled microbiomes on their own.

---
Tags: ngs, microbiology
---


# Question
Title: Basecalling process carried out by full_1dsq_basecaller.py/.exe
I am unable to understand what exactly the full\_1dsq\_basecaller.py/.exe script does. Is it giving the consensus sequence for linear and complementary strand, or it is just giving the reads that came from the complimentary strand?

If we are getting two fastq files as output which one we should use for alignment?

Thank you in advance. I greatly appreciate any help

# Answer
> 3 votes
My understanding is that the 1D$^2$ process is carried out in three stages:

1. Carry out a 1D calling of reads (i.e. what the standard Albacore does at the moment)
2. Link reads that are similar to each other (in reverse-complement orientation), from the same channel, and have a very short pore transition time
3. Create a 2D consensus base call, using information from the two linked reads

The two outputs from the basecaller are the 1D calls from step 1, and the consensus calls from step 3. Ideally just the 2D consensus calls should be used for assembly (assuming it represents a random distribution of reads), but if that gives too low coverage (or is not random enough), then the non-converted reads from the 1D calls could be included as well.

Because the linking is not a physical link, it's not a good idea to use 1D$^2$ for amplicon reads; the likelihood of two complementary reads from different templates appearing one after another is too high.

---
Tags: nanopore, albacore, 1d2-reads
---


# Question
Title: How many reads has my sequencing run produced on minion?
I am running a 48 h sequencing protocol on a FLO-MIN107, on MinION, with DNA library-prepped with SQK-RAD004. MinKNOW is installed on my Mac and controlling the MinION. I have found the following histogram in MinKNOW. What do the numbers mean? In particular, how many reads have been produced? I understood that the events are consecutive electrical signals and are then transformed into reads? Do I have 4029 reads? The run has started 24 hours ago. 4029 reads surely is really low, isn't it? What does 4,120,722 mean?

The Platform QC of the flow cell was as follows:

* group 1: 480
* group 2: 374
* group 3: 201
* group 4: 59
* total = 1114

# Answer
That's indeed the number of reads and that's quite low. How was your pore occupancy (number of pores sequencing) and your flow cell QC (number of pores good enough for sequencing)? How was your library concentration?

Events is an approximation of nucleotides: it's a guess based on the current signal. As far as I know if you multiply this by 1.8 you get approximately your throughput in nucleotides.

The 4,120,722 is the number of events in the currently selected (mouse hovering?) bin of the histogram.
> 3 votes

---
Tags: nanopore, ngs, minion, long-reads, minknow
---


title: Human Cell Atlas - Data availability

text:
# Question
Title: Human Cell Atlas - Data availability
A news item from 2017-10-18 on the website of the Human Cell Atlas states:

> In addition, the consortium today also announced the impending release of gene expression profiles from the first one million immune cells collected under the HCA, toward an initial milestone of collecting at least 30 million cells representing several tissues and organs for the atlas’ first draft. **These data, to be posted on an online repository by early November, will be freely available for researchers’ use**. (https://www.humancellatlas.org/news/14)

As of now, Feb 2018, I cannot find these data anywhere. Does anybody know where I can download the gene expression data of these 1 million immune cells or when it will be eventually published?

# Answer
> 6 votes
There is now a post on Biostars:

> We are starting to define the release process for the Human Cell Atlas.
>
> The first couple of HCA releases will be primarily single-cell RNA-Seq and we are interested to know what you would expect our releases to contain
>
> Please fill out our survey if you would like to help us define our releases
>
> https://www.surveymonkey.com/r/S6Y753P

And I suppose that answers the original question: the data has not been released yet and probably won't be for a while considering they are still discussing how to make it available.

Update (from the same post):

> The first production release of the DCP is happening in June.
>
> The first data release of HCA is targetting Q4 2018.

---
Tags: public-databases
---


# Question
Title: How to calculate the memory usage of storing kmers in RAM
I want to write a program in C++ that stores kmers in a hash or in a trie. How can I calculate how much RAM I would need for each type of data structure? For this application my kmers are strand-specific, so I cannot reduce the size complexity by considering canonical kmers.

For a hash, for example, I know that there are $4^k$ possible kmers of length $k$. I also am assuming that DNA character with an alphabet $\in \\{G,A,T,C\\}$ can be stored in two bits ($2^2 = 4$). Using this logic the kmer hash should occupy $(4^k) \cdot 2 \cdot k$ bits. Plugging this into Wolfram Alpha (`(4^13) * 2 * 13 bits to gigabytes`), I calculated the amount of RAM that some small kmer sizes would require to store the entire hash in memory in GB.

```
k = 13 requires    0.2181 GB
k = 15 requires    4.027  GB
k = 17 requires   73.01  GB
k = 19 requires 1306   GB

```

This escalated rather quickly, and I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM?

The kmers that I want to store are sometimes sparse compared to the entire space (for k=13, only 29 million or so out of 67 million possibilities).

# Constraints

* I do not need to store the number of occurences of a kmer, just its existence
* I cannot risk false positives. Bloom filters won't work in this case. :^/

# Answer
To follow up on Devon Ryan's answer, I thought it would be a little fun to write a Python script that demonstrates using a bit array to maintain a presence/absence table.

Note: I wrote a C++ port that includes a custom bitset implementation that can be sized at runtime. This and the Python script are available on Github: https://github.com/alexpreynolds/kmer-boolean

```
#!/usr/bin/env python

import sys
import os
import bitarray

# read FASTA
def read_sequences():
    global seqs
    seqs = []
    seq = ""
    for line in sys.stdin:
        if line.startswith('>'):
            if len(seq) > 0:
                seqs.append(seq)
                seq = ""
        else:
            seq += line.strip()
    seqs.append(seq)

# build and initialize bit array
def initialize_bitarray():
    global ba
    ba = bitarray.bitarray(4**k)
    ba.setall(False)
    sys.stderr.write("Memory usage of bitarray.bitarray instance is [%ld] bytes\n" % (ba.buffer_info()[1]))

# process sequences
def process_sequences():
    global observed_kmers
    observed_kmers = {}
    for seq in seqs:
        for i in range(0, len(seq)):
            kmer = seq[i:i+k]
            if len(kmer) == k:
                observed_kmers[kmer] = None
                idx = 0
                for j in range(k-1, -1, -1):
                    idx += 4**(k-j-1) * bm[kmer[j]]
                ba[idx] = True

def test_bitarray():
    test_idx = 0
    for j in range(k-1, -1, -1):
        test_idx += 4**(k-j-1) * bm[test_kmer[j]]
    test_result = ba[test_idx]
    if test_result:
        sys.stdout.write("%s found\n" % (test_kmer))
        sys.exit(os.EX_OK)
    else:
        sys.stdout.write("%s not found\n" % (test_kmer))
        sys.exit(os.EX_DATAERR)

def main():
    global k
    k = int(sys.argv[1])

    global bm
    bm = { 'A' : 0, 'C' : 1, 'T' : 2, 'G' : 3 }

    read_sequences()
    initialize_bitarray()
    process_sequences()

    try:
        global test_kmer
        test_kmer = sys.argv[2]
        if len(test_kmer) == k:
            test_bitarray()
        else:
            raise ValueError("test kmer (%s) should be of length k (%d)" % (test_kmer, k))
    except IndexError as err:
        keys = list(observed_kmers.keys())
        for i in range(0, len(keys)):
            sys.stdout.write("%s found\n" % (keys[i]))

    sys.exit(os.EX_OK)

if __name__== "__main__":
    main()

```

Note that this doesn't look at canonical kmers, *e.g.*, `AG` is considered a distinct 2mer from its reverse complement `CT`.

To use this script, you pipe in your FASTA, specify the *k*, and an optional kmer that you want to test for presence/absence, *e.g.*:

```
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 2 CG
Memory usage of bitarray.bitarray instance is [2] bytes
CG found

```

Or:

```
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 3 AAA
Memory usage of bitarray.bitarray instance is [8] bytes
AAA not found

```

Or if the optional test kmer is left out:

```
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 5
Memory usage of bitarray.bitarray instance is [128] bytes
CATTC found
ATTCT found
TTCTC found
TCTCG found
CTCGG found
TCGGG found
CGGGA found
GGGAC found
TTATA found
TATAT found
TTTAT found
TTATT found
TATTA found
ATTAG found
TTAGA found
TAGAC found
AGACC found
GACCT found
ACCTC found
CCTCT found

```

Or for the ~67M kmers in a 13-mer set, for which a roughly 8.4MB bit array is reserved:

```
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 13
Memory usage of bitarray.bitarray instance is [8388608] bytes
TTTATTAGACCTC found
TTATTAGACCTCT found

```
> 7 votes

# Answer
> I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM?

Hash table based k-mer counters only keep k-mers seen in data. For $16\<k\le32$, you need 64-bit integers to store a k-mer. Given $n$ distinct k-mers, a naive implementation with open addressing hash table would roughly take $2n\cdot 64/8=16n$ bytes. We assume the table is half full – this is where factor 2 comes from. There are various tricks to reduce the memory footprint. You'd better read papers. On a side note, don't use std::unordered\_map. It wastes extra memory.

A bitmap described by other answers takes $4^k/8=2^{2k-3}$ bytes, which is a better choice sometimes, but rarely used by k-mer counters as it does not work with long k, either. Trie or full-text indices are overkilling for fixed k-mers and are usually slower.

By the way, this question reminds me of your question on efficient reverse complement. Note that for k-mer counting, you traverse both strands at the same time. You should not reverse complement the whole k-mer.
> 6 votes

# Answer
If you simply need a set of present/absent k-mers and your kmer length is short, then convert the sequence length to an integer (a 13-mer would require 26 bits or 4 bytes) and use that as an index into a bit array. Initialize the array to 0 and use 1 for "k-mer present". If you have ~67 million k-mers then your array is 67 million bits, which is only about 8.5 megs.

This will scale "OK", but eventually balloon to a huge size if you use larger k-mers. At that point you simply have to use other data structures. The reality is that most applications are fairly tolerant to small amounts of error.
> 5 votes

---
Tags: genome, k-mer, c++
---


# Question
Title: Albacore wheel not supported on this platform (Mac OS X El Capitan)
I download the Albacore wheel from Oxford Nanopore.

I activated my conda environment that has Python 3 as the default python version.

I tried to install Albacore with

```
pip3 install ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl

```

but got the following error:

```
ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl is not a supported wheel on this platform.
You are using pip version 8.1.2, however version 9.0.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.`

```

I then upgraded my `pip` with:

```
pip install --upgrade pip
Requirement already up-to-date: pip in /Users/.../anaconda3/envs/python3env/lib/python3.6/site-packages`

```

So it seems that my pip was already up to date. I do not understand why it complains about my `pip` if I am using `pip3` anyways. I tried to install Albacore again but I got the same message as before.

Can anyone help me install Albacore on my Mac El Capitan laptop, please?

Update based on Wouter's suggestion:

```
python -m pip install ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl
Processing ./ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl
Requirement already satisfied: python-dateutil in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Requirement already satisfied: setuptools in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Requirement already satisfied: h5py in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Collecting progressbar33 (from ont-albacore==2.1.10)
Requirement already satisfied: numpy>=1.12.1 in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Collecting ont-fast5-api>=0.4.1 (from ont-albacore==2.1.10)
  Using cached ont_fast5_api-0.4.1-py2.py3-none-any.whl
Requirement already satisfied: six>=1.5 in /Users/userName/anaconda3/lib/python3.6/site-packages (from python-dateutil->ont-albacore==2.1.10)
Installing collected packages: progressbar33, ont-fast5-api, ont-albacore
Successfully installed ont-albacore-2.1.10 ont-fast5-api-0.4.1 progressbar33-2.4

```

At first sight, it looks like Albacore was successfully installed but if I run: `read_fast5_basecaller.py --help`

I get:

```
Traceback (most recent call last):
  File "/Users/userName/anaconda3/bin/read_fast5_basecaller.py", line 69, in <module>
    from albacore.input_utils import list_input_files
  File "/Users/userName/anaconda3/lib/python3.6/site-packages/albacore/input_utils.py", line 7, in <module>
    from albacore import log_utils
  File "/Users/userName/anaconda3/lib/python3.6/site-packages/albacore/log_utils.py", line 5, in <module>
    from albacore.pipeline_core import get_debug_level
ImportError: dlopen(/Users/userName/anaconda3/lib/python3.6/site-packages/albacore/pipeline_core.cpython-36m-darwin.so, 2): Library not loaded: /usr/local/opt/python3/Frameworks/Python.framework/Versions/3.6/Python
  Referenced from: /Users/userName/anaconda3/lib/python3.6/site-packages/albacore/pipeline_core.cpython-36m-darwin.so
  Reason: image not found

```

Is there a library missing?

# Answer
Not sure about `conda` as I don't use it, but I can highly recommend `pipenv` which is the python.org recommended package management system for python.

I have tested the following on my Mac and it worked like a dream

```
mkdir albacore && cd albacore
mv ../ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl .

# install pipenv if you don't have it
brew install pipenv  # or pip3 install pipenv

# the following will create a virtual environment and install albacore into it
pipenv install ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl

# activate the virtual environment
pipenv shell

read_fast5_basecaller.py -h

```

And that should bring up the albacore help menu and you're ready to go.
> 2 votes

---
Tags: nanopore, minion, software-installation, albacore
---


# Question
Title: How to filter info score post-imputation?
this is a very basic question, but I cannot find it explicitly stated anywhere: When exactly should I filter the info scores after imputation?

I did imputation in Impute2 and am planning on using SNPtest to analyze, but cannot tell if I should filter the scores in Impute2, SNPtest, Gtool, QCtool, or through python.

I'm hoping to get rid of the poorly imputed results (with an info score cutoff \<0.5) in a population of ~4000 (half cases, half controls). After this, I'm running an association analysis -- probably through PLINK, but I'm still figuring out the best way to go from Impute2 to PLINK.

Any help would be greatly appreciated!!

# Answer
In simple GWAS setups, each SNP is analyzed independently. In those cases you can filter out SNPs with poor INFO scores at any point.

For analyses that combine information across SNPs (for example FaSTLMM), I would recommend filtering out SNPs before running the step that combines information across SNPs.

So, if you filter right after running Impute2, you should be safe.
> 2 votes

# Answer
Depending on the size of your samples (1K, 10K, 100K?) or the size of imputation region (1 gene, 1 chromosome, whole genome?) imputation files even after zipping are usually huge.

Instead of filtering SNPs and creating yet more gazillions of subset files. I prefer to have `excludeSNP.txt` file list of SNPs. This list then used as input for further analysis using *SNPTEST* flag `-exclude_snps_g`. This list is part of IMPUTE2 output or could be additional list of SNPs that we wish to exclude for other reasons. In short, filter at the point of analysis not the imputated files.

If *info* file is missing we can run *SNPTEST* with `-summary_stats_only` flag, which gives you the info score.

Here is an old still relevant post at BioStars post-imputation QC.
> 1 votes

---
Tags: gwas, imputation
---


# Question
Title: Why would someone use a CRAM instead of a BAM?
I had this question from a graduate student yesterday, and I was stuck.

What should I say? Why use a CRAM instead of a BAM?

1. When is it a good idea to use a CRAM instead of a BAM?
2. When is it a bad idea?

# Answer
> 17 votes
1. Whenever you want to save space (this can be a substantial savings). Until quite recently (samtools/htslib 1.7), only CRAM supported long CIGAR strings.
2. If you need to guarantee that any random obscure downstream program will be able to handle it.

Uptake of CRAM has been pretty slow. Java programs using htsjdk (e.g., picard, IGV and GATK) have only relatively recently added support for CRAM. If you need to use an old version of those for some very odd reason then CRAM may not be supported.

There are a lot of programs written in python that use pysam to open BAM files and these should, theoretically, support CRAM. The issue is that some of the functions may fail and one can't assume that authors will have always written the code needed to handle this. I'll use deepTools as an example, since I'm one of its developers. One of the things about CRAM files is that they (by default) are made such that you require a reference genome in order to construct the sequence field in each alignment. This works fine if you're using a standard genome (htslib, via pysam, can fetch many standard genomes from the web automatically), but if you're not, then you need to specify a fasta file to use for decompression. Every tool, then, needs to add an option for this. With pysam 0.14 and htslib 1.7 this can be circumvented by not decompressing the sequence, but behavior has to be explicitly requested.

Another issue is that many tools will use features from the file index, such as the `.mapped` accessor, to get the number of mapped reads in a file. CRAM files contain very very little information, so this then fails. Consequently, tool authors need to check for CRAM files and both derive and propagate this information through their functions if it's needed. This can be a time-consuming task (e.g., it took me a couple days to get this implemented in deepTools). Relatedly, `samtools idxstats` is useless on CRAM files, since there are no statistics stored in the index.

That having been said, it's likely that CRAMs slowly gaining acceptance will eventually make it the standard. It's already a convenient archival format, it's just a matter of time before users can assume that most analysis programs are written to handle it.

---
Tags: sam, file-formats, cram
---


# Question
Title: Can I index a compressed FASTA file using STAR?
I am using STAR to align RNA-seq reads to a reference genome. Before the alignment, I need to generate an index of the reference genome. I use the following code to generate the index successfully:

`STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa --sjdbGTFfile ref.gtf --sjdbOverhang 100`

This works fine. However, I would like to keep my reference genome compressed to save disk space. So I am trying the following command:

`STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa.gz --readFilesCommand "gunzip -c" --sjdbGTFfile ref.gtf --sjdbOverhang 100`

but I get the following error:

`EXITING because of INPUT ERROR: the file format of the genomeFastaFile: ref.fa.gz is not fasta: the first character is '' (31), not '>'.
 Solution: check formatting of the fasta file. Make sure the file is uncompressed (unzipped).`

I am using the readFilesCommand successfully with compressed RNA-seq fastq files. Does anybody know if there is a similar option to use compressed references? Is there a workaround using Unix commands (maybe piping?) or do I need to decompress the reference, index it and then compress it again?

# Answer
I can think of two possible workarounds. The simplest would be to try using `zcat` instead of `gunzip -c`, in case that works. Assuming it doesn't, you can use process substitution, if your shell supports it. Something like this:

```
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star \
    --genomeFastaFiles <(zcat ref.fa.gz) --sjdbGTFfile ref.gtf --sjdbOverhang 100

```

The `<()` creates a file descriptor and passes that to the command. It is very useful for passing streams to commands that expect files.

If your shell doesn't support this construct, assuming you are on a system that supports named pipes (Linux does), you can do the same thing a bit more manually:

1. Create a named pipe (FIFO)

```
mkfifo foo

```
2. Set it to read the uncompressed genome

```
zcat ref.fa.gz > foo &

```

The & is important since the command above will hang until something starts reading from `foo`. Alternatively, you can just run this in one terminal and the next command in another.
3. Read the input from the pipe

```
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star \
    --genomeFastaFiles foo --sjdbGTFfile ref.gtf --sjdbOverhang 100

```
> 5 votes

---
Tags: sequence-alignment, genome, star, genome-index
---


# Question
Title: Subset smaller BAM to contain several thousand rows from multiple chromosomes
There are many cases whereby I would like to subset a BAM to create a small file in order to work with (e.g. algorithmic testing, debugging, etc.)

Normally I do the following, which will subset the BAM `file.bam` and keep the header

```
samtools view -H file.bam > header.sam
samtools view file.bam | head -n 5000 | cat header.sam - | samtools view -Sb - > file.unique.bam

```

In this case, I would like 5000 rows in chromosome 1 and 5000 in chromosome 2.

I could first try grepping by individual chromosome, and then combining the two SAMs

e.g. here's complete BAM with grepped chr1 and (incorrect but complete) header

```
samtools view -H file.bam > header.sam
samtools view file.bam | grep "chr1" | cat header.sam - | samtools view -Sb - > file.unique.bam

```

but then I have two problems:

(1) I may not be grepping the alignments to chromsome 2---there may exist BAM rows which contain 'chr2' but are not alignments.

(2) I think one must manually edit the header. There's probably no way around this.

Is there an easy way, Bioinformatics SO?

# Answer
If you're not too hung up on exact numbers like 5000 reads then you can do that with a single samtools command:

```
samtools view -bo subset.bam -s 123.4 alignments.bam chr1 chr2

```

That will select 40% (the `.4` part) of the reads (`123` is a seed, which is convenient for reproducibility). The convenient part of this is that it'll keep mates paired if you have paired-end reads. For 5000 reads per chromosome just change the `.4` part to a sufficiently small number.

In general you don't really need to subset the header. Some tools will perform a bit better if you do, but you'll generally get the same results regardless.
> 3 votes

# Answer
You can use SAMsift:

```
samsift \
   -i file.bam \
   -0 'c={"chr1":5000,"chr2":5000}' \
   -f 'c[RNAME]>0' \
   -c 'c[RNAME]-=1' \
   -m nonstop-remove

```

Explanation:

* `-i file.bam` – input file
* `-0 'c={"chr1":5000,"chr2":5000}'` – initialization (create a count-down dictionary for the chromosomes of interest)
* `-f 'c[RNAME]>0'` – filtering criterion (is the counter for the current chromosome still \>0?)
* `-c 'c[RNAME]-=1'` – code decrementing the counter of the current chromosome (5000 -\> 4999 -\> ... -\> 1 -\> 0 -\> -1 -\> ...)
* `-m nonstop-remove` – remove lines causing Python errors and don't stop (in this case, an error can be caused by accessing a non-existing counter for another chromosome, e.g., for chr3)

See the SAMsift readme for more info.
> 3 votes

# Answer
I'd generally recommend Devon Ryan's answer. However, if you do care about having the same number from each Chr, you could use the following python/pysam code (this will output approx 5000 from every Chr):

```
from pysam import AlignmentFile
from random import random

nreads = 5000

infile = AlignmentFile("mybam.bam")
outfile = AlignmentFile("outbam.bam", "wb", template=infile)

for chr in infile.reference_names:
     reads_in_chr = infile.count(chr)
     frac = float(nreads)/reads_in_chr
     count = 0
     #-- Replace this loop for the exactly 5000 first reads --#
     for read in infile.fetch(chr, multiple_iterators=True):
        if random() < frac:
            count += 1
            outfile.write(read)
     print("ouputted %i reads from %s" %(count, chr))
outfile.close()

```

This uses all chromosomes. If you only want to use, say chr1 and chr2 replace `infile.reference_names` with `["chr1","chr2"]`.

If you want exactly 5000 reads from each chr, but don't care if they are the first ones or not, then you could replace the inner for loop with:

```
for read in infile.fetch(chr, multiple_iterators=True):
    count += 1
    if count >= nreads:
        break
    outfile.write(read)

```

if you wanted the mates of those reads as well, you could add the following after `outfile.write(read)`:

```
mate_read = infile.mate(read)
outfile.write(mate_read)

```

Note that this is slow.
> 2 votes

# Answer
using **samjdk**: http://lindenb.github.io/jvarkit/SamJdk.html

```
$ java -jar dist/samjdk.jar --body -e \
  'Map<String,Integer> c=new HashMap<>(); public Object apply(SAMRecord r) {int n=c.getOrDefault(r.getContig(),0);if(n>=5000) return false; c.put(r.getContig(),n+1); return r;}' \
 input.bam

```
> 1 votes

# Answer
If you just want to create a truncated BAM file with a header then you can significantly simplify your original code:

```
(samtools view -H input.bam; samtools view input.bam | head -5000) | samtools -bo output.bam

```

This command avoids the intermediate file and few gratuitous intermediate command invocations, at the cost of a subshell (invoked by `(…)`).

As above, but with formatting:

```
(
    samtools view -H input.bam
    samtools view input.bam | head -5000 # (*)
) \
| samtools -bo output.bam

```

… this can of course be extended to filter by multiple chromosomes by replacing the line marked with `(*)` above by one or multiple lines that subset by chromosome name (`samtools view input.bam chr*x*`, no need for `grep` if you have indexed the original BAM file!).
> 1 votes

---
Tags: sam, samtools
---


# Question
Title: How to you find read statistics for bam files with terminal mismatch?
I have a bam file aligned to a reference sequence. I am interested in finding the number of reads for read length 18 to 34, where the last base is mis-matched. I got this bam file from bowtie alignment allowing 3 mismatches (-v 0) option, so I am expecting some of the reads to have terminal bases mis-matched. How can we count the number of sequences of a specific length with terminal mismatch in bam file?

\**UPDATE:

I have created a bowtie aligned .sam file with this command: `bowtie -S -p 18 /media/owner/ref.fa test1.fastq test1.sam`

which I then converted to .bam file like this: `samtools view -bS -o test1.bam test1.sam`

From this test1.bam, I was able to get the total number of plus and minus strand of reads aligned to the reference loci using bedtools coverage maping (for positives reads and negative reads) like this: `bedtools coverage -a my_region.bed -b test1.bam -bed -d -s | gzip > test1_region_pos.tsv.gz`

Since I used default bowtie option which allows for three mismatches, I want to know how many reads are aligned with terminal mismatches (i.e,last three, last two and last one base mismatched at 3' region of the aligned reads at their loci position).

# Answer
Could you provide some examples of lines that should be counted? I'm not sure how mismatches are encoded in the version of Bowtie you're using. As far as I remember, Bowtie doesn't distinguish between `=` and `X` in CIGARs (and uses `M` for both), therefore, mismatches must be detected from the MD tag.

You can use SAMsift for filtering the alignments. The following filtering criterion should keep only reads with read length from 18 to 34 that have an MD tag with `0` as the last character, which must be preceded by a non-number (this should signalise a mismatch at the last position).

```
samsift \
  -i test1.bam \
  -f '18 <= len(SEQ) <= 34 and MD[-1]=="0" and not MD[-2].isdigit()' \
  -m nonstop-remove

```

The obtained alignments can be counted by adding `| samtools view | wc -l`.

Explanation:

* `-f ...` – filtering criterion
* `18 <= len(SEQ) <= 34` – the allowed read lengths
* `MD[-1]=="0"` – the last character of the `MD` tag is `0`
* `not MD[-2].isdigit()` – the penultimate character of the `MD` tag is not a digit
* `-m nonstop-remove` – remove alignments without any `MD` tag (probably unaligned, they would cause an error).
> 2 votes

---
Tags: sam
---


# Question
Title: Why is my STAR reference genome indexing aborting on my GNU/Linux server but not on my Mac OS X laptop?
I am running the following command to index my genome:

`STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa --sjdbGTFfile ref.gff3 --sjdbGTFtagExonParentTranscript Parent --sjdbOverhang 100`

and I get the following output:

```
Feb 27 10:27:37 ..... started STAR run
Feb 27 10:27:37 ... starting to generate Genome files
Feb 27 10:27:40 ... starting to sort Suffix Array. This may take a long time...
Feb 27 10:27:42 ... sorting Suffix Array chunks and saving them to disk...
Feb 27 10:28:26 ... loading chunks from disk, packing SA...
Feb 27 10:28:30 ... finished generating suffix array
Feb 27 10:28:30 ... generating Suffix Array index
terminate called after throwing an instance of 'std::bad_alloc'
  what():  std::bad_alloc

```

Searches for the GNU/Linux error lead to some pages about C++ that I don't understand.

Submitting it as a batch job works fine but I would like to know why it does not work on a node. It is quite practical for me to test small code like this without submitting a batch job.

The command above works fine on my laptop running Mac OSX.

The STAR version is STAR\_2.5.3a on both computers: STAR\_2.5.3a

The memory usage of the node is the following:

```
total            used       free     shared    buffers       cached

251              251         0         0       0             232

```

so, is there 0 memory for me to use on this node?

Why is my indexing aborting?

# Answer
A `std::bad_alloc` error tends to mean that you've run out of memory. It's not unusual for head nodes to be fairly limited in capacity, so presumably that's the issue.
> 1 votes

---
Tags: sequence-alignment, rna-alignment, linux, star, genome-index
---


# Question
Title: Create galaxy local toolshed using Planemo
I tried to create galaxy local toolshed using planemo. I used command as follows,

1. create virtual env

```
virtualenv .venv; .venv/bin/activate

```
2. Install planemo

```
pip install planemo

```
3. To configuring a shed account

```
planemo config_init

```
4. Create a repository

```
planemo shed_init --name=seqtk_seq \
                  --owner=planemo \
                  --description=seqtk_seq \
                  --long_description="Tool that converts FASTQ
                                      to FASTA files using seqtk" \
                  --category="Fastq Manipulation"
planemo shed_create --shed_target local

```

This produces the following error:

```
Failed to include files for [{'source': '**'}]
Problem encountered executing action for one or more repositories.
(.venv) varuna@varuna-Inspiron-15-3567:~/git/galaxy_stable/planemo

```

# Answer
`--shed_target local` will only work if you're running a Galaxy toolshed on your local computer, since it's reading the `.shed.yaml` file you made with `planemo shed_init` and trying to create a repository with that name or a local (in this case) server. Planemo itself won't start a Galaxy toolshed, since it's fairly unusual for people to want to do testing involving creating local repositories. What you probably want to do instead is use `planemo test` instead. Once the planemo testing passes then use the various `shed*` commands to upload your wrapper to the toolshed or testtoolshed.
> 1 votes

---
Tags: software-installation, galaxy
---