# Question
Title: How do I estimate branch length of a maximum parsimony phylogeny in R?
I am trying to estimate a phylogenetic tree in R with ape and phangorn. With `pratchet()` and `optim.parsimony()`, I can estimate a phylogenetic tree, but branch lengths are not estimated. Is there a way to estimate a phylogeny with branch lengths with the following made-up data:

```
x <- data.frame("L1" = c(0,0,1,0,0,0,1,1,1,1,0,0,0,0,1), "L2" = c(0,1,1,0,1,1,0,1,1,0,1,1,0,1,1),"L3" = c(1,1,1,0,0,0,1,1,1,0,0,0,1,1,1), "L4"=c(0,0,1,1,0,0,0,0,1,1,0,0,0,0,1))

```

# Answer
> 2 votes
@Nathan S. Watson-Haigh's answer does give branch lengths, but NOT parsimony branch lengths, which may be what the OP wanted. `nnls.phylo` is giving distances calculated from a pairwise distance matrix. To get the parsimony branch lengths (i.e. number of changes mapped on to each branch), score the tree against the dataset using `acctran`:

```
treeRatchet <- pratchet(primates)
treeRatchet.BL <- acctran(treeRatchet,primates)
plot(treeRatchet.BL,main='Phylogeny (Parsimony branch lengths)')

```

Note that the branch lengths in the parsimony phylogram are 'chunkier' than those in the nnls.phylo phylogram -- that's because

1. they are all integers, and
2. they aren't corrected with reference to a model of substitution that considers multiple hits or base composition.

# Answer
> 1 votes
# Terminology

First, a couple of corrections to your use of terminology. A tree *must* have branch lengths associated with it to be called a **phylogenetic tree** as it represents both the relationships between OTU/taxa/leaves *AND* the distances between them (i.e. the branch lengths. A tree where branch lengths have no meaning is called a **cladograms** or a **dendrogram**. You can tell if a tree is a cladogram since all the OUTs/taxa/leaves end at the same position and internal branch lengths are usually (not always) of equal length. Whereas a phylogeny the OTUs/taxa/leaves will typically look "staggered" *and* internal branches will be of different lengths.

Therefore, If the tree you have estimated, truly doesn't have branch lengths estimated and associated with it, it is actually a cladogram.

# Example Derived from `phangorn` Vignette

```
library(phangorn)

# Get example data
fdir <- system.file("extdata/trees", package = "phangorn")
primates <- read.phyDat(file.path(fdir, "primates.dna"), format = "phylip")

# Calculate pairwise distances between OTUs using Maximum Likelihood
dm <- dist.ml(primates)

# Use the distances to cluster the OTUs using UPGMA algorithm
# This creates a dendrogram
treeUPGMA <- upgma(dm)
plot(treeUPGMA, main="UPGMA")
dev.off()

```

```
# Search for different, similar trees to the initial UPGMA tree
# looking for a tree with the lowest parsimony score
# There are 2 search algorithms implemented:
#   optim.parsimony()
#   pratchet()
# These are still cladograms
treePars <- optim.parsimony(treeUPGMA, primates)
treeRatchet <- pratchet(primates, trace = 0)
plot(treePars, main="optim.parsimony()")
plot(treeRatchet, main="pratchet()")
dev.off()

```

```
# Calculate branch lengths using the most parsimonious tree identified by
# the pratchet() searchh algorithm
phylogeny <- nnls.phylo(treeRatchet, dm)
plot(phylogeny, main="Phylogeny")
dev.off()

```

# Choices, Choices, Choices

As you can see from the above example, there are lots of steps to calculating a phylogenetic tree. At each of these steps there may be multiple algorithms to choose from. You choice will depend on the type of data you are analysing, the size of the data set, your available compute resources and the algorithms implemented and available through the tools you use.

Ultimately, you will have to know a fair bit about these algorithms to make an informed choice.

---
Tags: r, phylogenetics, phylogeny
---


# Question
Title: Relationship between sequencing lane and ngs dataset
I am fairly new to NGS data analysis and I am struggling to understand the exact relationship between a sequencing lane and an NGS dataset. I should add I don't work in the lab, I only do bioinformatics.

I understand the basics of how to make NGS libraries (very briefly: shearing, ligation of adapters, optional PCR, ligation to solid surface, following by sequencing-by-synthesis or other method). What I am struggling with is the visualisation of what happens once you have made your libraries. My questions are:

* If you have used adapters containing barcode sequences (that identify each patient), you would presumably end up with one single tube with your libraries in it right?
* If so, you then "put" the contents of that tube along with DNA polymerase and nucleotides into one single lane of a flow cell? Your provider would then give you raw FASTQ files that have been identified using this barcode? So you would effectively received one (or 2, depending if it's single/paired-end sequencing) FASTQ file per patient right?
* If you do not use adapters with barcodes, does that mean you would have to use a sequencing lane per patient? If so, that would vastly increase the price you pay your provider, would it not?

Apologies if these are basic questions.

# Answer
> If you have used adapters containing barcode sequences (that identify each patient), you would presumably end up with one single tube with your libraries in it right?

Correct, you have to prep each sample into a library.

> If so, you then "put" the contents of that tube along with DNA polymerase and nucleotides into one single lane of a flow cell?

You could, or you could pool it with other samples (assuming they have compatible indices). As @swbarnes2 mentioned, you don't see a lot of samples run on a single lane these days. Sequencing centers like to split things across lanes in case one lane fails for some reason or there are bubbles in the flow cell or other technical issues.

> Your provider would then give you raw FASTQ files that have been identified using this barcode?

Typically, yes, and if your sample was run on multiple lanes they'll often be merged for you automatically. This isn't always the case. I've seen a number of people receive samples either still split by lane (annoying, but you can just `cat` things together) or, much worse, not demultiplexed at all. In the latter case, people often receive 3 (or 4) fastq files per lane, one for each read and another for each index read. It's rare that sequencing facilities do this, but it happens.

> If you do not use adapters with barcodes, does that mean you would have to use a sequencing lane per patient? If so, that would vastly increase the price you pay your provider, would it not?

You could sequence a sample per lane. Whether that makes economic sense is up to you to decide.
> 4 votes

# Answer
It is worth noting that the Illumina NovaSeq and NextSeq address all of the lanes of the flow cell from a single loading point (by default). So on those instruments, you end up with all libraries loaded across all lanes. Therefore for each index you will have L001, L002 (+ L003 and L004 if appropriate) FASTQ files (and R1/R2, again if appropriate).

As @Devon mentions, these files can be cat-ted together, but for some processes modelling the lane effect does matter (BQSR, for example).
> 3 votes

# Answer
Indicies matter, lanes do not, unless you have a screw up on the instrument that messes up some lanes but not all. I'm not sure if anyone skips indices just because they put only one sample on a lane nowadays...it's valuable to be able to use the indices to catch sample mix-ups.
> 0 votes

---
Tags: ngs, illumina
---


# Question
Title: Multiway intersection of 8 lists of genes
I have 8 lists of differentially expressed genes for 8 time points (2, 4, 6, 8, 10, 12, 14, 16 hours), the number of genes in each of these lists are not equal.

For example there are 630 common genes between h2 and h4 time points but I have 1143 genes in h2 list and 768 genes in h4 list , so common genes in h2 vs h4 would be 630/1143 and common genes h4 vs h2 would be 630/768.

That would would be great if I have a code to produce a matrix or plot of all possible pairwise combinations of my 8 time points.

# Answer
> 3 votes
EDIT: Sorry, I did not see the R tag and wrote the code in bash. I still left the bash example below and added an R version.

Depends on what you want to do, but there are two common ways in which you could store this. Either as wide format:

```
      h2  h4  h6 ...
Gene1 1   1   0
Gene2 0   1   1
Gene3 1   1   1
...

```

Where 1 denotes presence and 0 absence, or as long format:

```
      timepoint
Gene1 h2
Gene1 h4
Gene2 h4
Gene2 h6
Gene3 h2
Gene3 h4
Gene3 h6
...

```

The best format probably depends on what you want to do. Since you did not specify a method/language, I took the liberty to write a few lines of bash to go from the gene lists to the matrices. Bash examples assume your files are named h2, h4 etc and are tab separated with gene names in the first column. R examples assume `tables` is a list of your data frames and `timepoints` is a vector `c("h2", "h4", ...)`.

**Generate table in wide format:**

1. Concatenate gene names of all tables to a separate file
2. Extract unique values in this file
3. Perform outer join between the separate file and each table

*R:*

```
all_names <- c()
for (tp in timepoints){
  all_names <- append(all_names, tables[[tp]]$genes)
}
all_names <- unique(all_names)

wide_table <- matrix(nrow=length(all_names), ncol=(1+length(timepoints)))
wide_table <- as.data.frame(wide_table)
colnames(wide_table) <- append("genes", timepoints)

for (tp in timepoints){
 tmp <- tables[[tp]]
 tmp$tp <- 1
 tmp <- merge(wide_table, tmp, by="genes", all=T)
 wide_table$tp <- tmp$tp
}

```

*Bash:*

Note: you will might need to add headers to this table

```
for timefile in $yourfiles;do
  # make a file with all gene names
  awk -v tp=$timefile '{print $1}' >> names_file
done

sort names_file | uniq > tmp && mv tmp names_file

for timefile in $yourfiles;do
  join -a1 -a2  -j 1 -o 0,1.2,2.2 -e 0 $names_file <(awk '{print$1,1}' $timefile | sort)
done

```

**Generate table in long format:**

1. For each table, add a column with timepoint
2. Concatenate all tables

*R:*

```
df <- data.frame()
for (tp in timepoints){
  tables[[tp]][[timepoint]] <- tp
  df <- rbind(df, tables[[tp]] )
}

```

*Bash:*

```
   for timefile in $yourfiles;do
      awk -v tp=$timefile '{print $1,timefile}' >> long_table
    done

```

EDIT2: To follow up on this, here is a simple way (using the package dplyr) to count the number of common genes between all possible combination of time points. Assuming you have your data in long format:

```
df <- data.frame(genes=c("G1","G1","G2","G3","G3"), tp=c("h2","h4","h2","h4","h6"))

  genes tp
1    G1 h2
2    G1 h4
3    G2 h2
4    G3 h4
5    G3 h6

```

You can get the co-occurrences using :

```
library(dplyr)
# Make timepoint ordinal factor
df <- mutate(df, tp = factor(tp, ordered=T))
df <- df %>%
  # Compute all pairwise co-occurrences
  full_join(df, by="genes") %>%
  group_by(tp.x, tp.y) %>%
  summarise(N_genes = length(unique(genes))) %>%
  # Exclude self pairs (gene in h2 and h2)
  filter(tp.x != tp.y) %>%
  # Check for duplicate pairs in reverse order (h2-h4 and h4-h2)
  mutate(comb=ifelse(tp.x < tp.y,
                     yes = paste(tp.x, tp.y),
                     no  = paste(tp.y, tp.x)))
# Remove the duplicate pairs
df <- df[which(!duplicated(df$comb)),]

```

And visualise it using for example a heatmap: (ofc it looks stupid because I have a tiny sample where not all combinations are even present)

```
# Visualize number of genes present in each combination
library(ggplot2)
ggplot(df, aes(x=tp.x, y=tp.y, fill=N_genes)) + geom_tile()

```

# Answer
> 3 votes
An UpSet Plot is a useful way of visualising these kind of multi-dimensional overlaps. There's an R package for producing them, with a reasonable set of examples.

The input is a binary matrix indicating set membership:

```
         h2    h4    h6...
GeneA     0     1     0
GeneB     1     1     0
GeneC     0     1     1
...

```

And the output is a pair of bar charts, showing the size of various sets, with a dot matrix showing the set membership. An example from the UpSetR documentation:

```
library(UpSetR)
mutations <- read.csv(system.file("extdata", "mutations.csv",
                      package = "UpSetR"), header=T, sep = ",")
upset(mutations, sets = c("PTEN", "TP53", "EGFR", "PIK3R1", "RB1"),
      sets.bar.color = "#56B4E9", order.by = "freq",
      empty.intersections = "on")

```

# Answer
> 1 votes
I would simply do this...

```
# sepcify your timepoints, each timepoint is a list of genes!!
timepoints <- list(h2,h4,h6,h8, h10,h12,h14,h16)
# initialize an empty array that you will fill with your timepoints
vect = vector()
# go over all combinations and fill the array with the set intersection values you want (I guess...)
for (i in timepoints) {
  for (j in timepoints) {
    vect <- c(vect, length(intersect(i,j))/length(i))}}
# reshape your array to a more ordered matrix
your_matrix <- matrix(vect, nrow=length(timepoints))
# assign columns and rows names
colnames(your_matrix) <- c("h2","h4","h6","h8","h10","h12","h14","h16")
rownames(your_matrix) <- c("h2","h4","h6","h8","h10","h12","h14","h16")

```

---
Tags: r
---


# Question
Title: How to concatenate FASTAs with respect to chromosome?
I have the following two small FASTAs, `dna1.fa` and `dna2.fa`:

```
# dna1.fa
>chr1
ACTGACTAGCTAGCTAACTG
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG

# dna2.fa
>chr1
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

```

If I concatenate these two files using `cat`, the piped result is as follows:

```
$ cat dna1.fa dna2.fa > total1.fa
$ cat total1.fa
>chr1
ACTGACTAGCTAGCTAACTG
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr1
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

```

Here, as expected, these are two files which have concatenated together. It seems this is usually what is done when one concatenates several FASTAs together in bioinformatics.

Question1: Is it possible to concatenate these "based on the chromosome"? Here is the desired result:

```
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

```

Note that these are exceptionally small FASTAs---usually one could be dealing with files of several GB in size.

Question2: The meta question here is, are there any standard bioinformatic analyses whereby the "concatenated by chromosome" FASTA would be more useful than the straightforward concatenated FASTA via `cat` alone?

# Answer
> 2 votes
Using the TblToFasta and FastaToTbl scripts I have posted before, combined with classic \*nix utils, you can do:

```
$ join -a1 -a2 <(FastaToTbl dna1.fa | sort) <(FastaToTbl dna2.fa | sort) |
    sed 's/ /=/2' | TblToFasta | tr = '\n'
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

```

### Explanation

The `FastaToTbl` will put everything on the same line, with the sequence's ID at the beginning of the line, then a tab and then the sequence (no matter how many lines of sequence you have). So, we do this to both files, sort them and pass that as input to `join`, which will join its input on lines where the 1st field is the same (by default). The `-a1 -a2` means "also print non-matching lines from both files" so we don't miss cases where a sequence is in only one of the two files.

Since the `join` will add a space between the two joined lines, we deal with that by replacing the second space on each line with a `=` (`sed 's/ /=/2'`) and then pass the whole thing to `TblToFasta` to revert to fasta format. Finally, we use `tr` to replace the `=` with a newline so you can see where each sequence came from as you show in your example.

For longer sequences where you want to avoid sorting the input, you can do:

```
$ FastaToTbl dna1.fa dna2.fa |
    awk '{
            seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2
         }
         END{
              for(chr in seq){
                printf ">%s\n%s\n", chr, seq[chr]
              }
         }'
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

```

### Explanation

Here, the tbl format sequences (see above) are passed to an awk script which does the work:

* `seq[$1] ? seq[$1]=seq[$1]"\n"$2 : seq[$1]=$2` : if the 1st field (`$1`, the sequence identifier) is already in the array `seq`, then set the value of `seq[$1]` to the current value (so whatever sequence has been already seen for this ID), a newline and the 2nd field (`$2`, the sequence). If it hasn't been seen before, then set the value of `seq[$1]` to the 2nd field, the sequence.
* `END{ for(chr in seq){ printf ">%s\n%s\n", chr, seq[chr]}` : after processing everything, go through each element in the array `seq` and print out the data in fasta format.

However, that will mean keeping both files in memory which may be a problem for huge sequences. You can avoid that by using temp files:

```
$ FastaToTbl dna1.fa dna2.fa |
    awk '{
            if(!seen[$1]++){
                print ">"$1 > $1".tmp"
            }
            print $2 >> $1".tmp"
        }' && cat *.tmp > out.fa && rm *.tmp
$ cat out.fa
>chr1
ACTGACTAGCTAGCTAACTG
GCATCGTAGCTAGCTACGAT
>chr2
GCATCGTAGCTAGCTACGAT
>chr3
CATCGATCGTACGTACGTAG
CATCGATCGTACGTACGTAG
>chr4
ATCGATCGATCGTACGATCG
ATCGATCGATCGTACGATCG
>chr5
CATCGATCGTACGTACGTAG

```

### Explanation

This is the same basic idea as the script above, only instead of storing the information in an array, we print it as it comes into temp files.

* `if(!seen[$1]++){print ">"$1 > $1".tmp"}` : if this is the first time we've seen this ID, print it as a fasta header into the file `$1".tmp"` (e.g. `chr1.tmp`).
* `print $2 >> $1".tmp"` append the current sequence to the corresponding temp file.

Once the temp files have been created, we concatenate them (`cat *.tmp > out.fa`) and delete them (`rm *tmp`).

---

As for why you would want to do this, an obvious example is reconstructing a reference genome from short sequences. You would want all the sequences of the same schromosome together. Of course, you would also need to assemble them in the right order, which you aren't doing here.

Alternatively, you might want to concatenate them this way so as to avoid having the extra, unneeded lines (`>chrN`) which would make your file larger. Or you would want to be able to retrieve all sequences with the same ID at once.

But I don't really know why you would want to do it. You're the one asking for it, after all! :)

---
Tags: fasta
---


# Question
Title: How do I create a for loop to filter through different FDR values?
This is probably a really quick/simple fix but I am a Noob to R programming/programming in general. I am trying to create a for loop that runs the test multiple times with different FDR values (5,10,20). How do I make it run 3 times with different values and then produce 3 different datasets with names like... `atac.glmtop.males.opening.fdr5`, `atac.glmtop.males.opening.fdr10`, `atac.glmtop.males.opening.fdr20`? Also is there a way to make the `write.table` function work so that it exports files with the fdr5, fdr10, fdr20 in the name?

```
fdr.nums = c("5","10","20")
fdr.nums.vector = c(5,10,20)

for (i in 1:length(fdr.list)){
  atac.glmtop.males.opening.fdr[i] <- topTags(
    atac.glmfit.males,
    n = nrow(atac.dge$counts),
    sort.by = "none",
    adjust.method = "fdr"
  )$table %>%
  tibble::rownames_to_column("peakco") %>%
  filter(logFC < 0) %>%
  filter(FDR < (fdr.vector[i]/100)) %>%
  merge(.,annotated.atac.pbmc.narrow, by="peakco" ) %>%
  tibble::column_to_rownames("peakco")

  atac.glmtop.males.opening.fdr[i]

  write.table(
    atac.glmtop.males.fdr[i].opening[,c("chr", "start","end")],
    file = paste0("~/Downloads/sexM.age.groupHO.fdr",i, ".opening.bed"),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
  write.table(
    atac.glmtop.males.fdr[i].opening["GeneName"],
    file = "~/Downloads/sexM.age.groupHO.fdr"+ i + ".opening.GeneNames.txt",
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
}

```

# Answer
> 2 votes
For loops in most languages can be done in a variety of ways. What you have currently coded, loops over values 1 to 3 by using

```
for (i in 1:length(fdr.list)){
...
}

```

However, you probably want to use a different form and loop over each value in that vector using something more like:

```
for (fdr in fdr.nums.vector){
...
}

```

This way you have convenient access to the FDR value in each loop via the `fdr` variable. Which you can then use in constructing your output filenames. Thus your code would become something like:

```
fdr.nums.vector = c(5,10,20)

for (fdr in fdr.nums.vector){
  atac.glmtop.males.opening.fdr[fdr] <- topTags(
    atac.glmfit.males,
    n = nrow(atac.dge$counts),
    sort.by = "none",
    adjust.method = "fdr"
  )$table %>%
  tibble::rownames_to_column("peakco") %>%
  filter(logFC < 0) %>%
  filter(FDR < (fdr/100)) %>%
  merge(.,annotated.atac.pbmc.narrow, by="peakco" ) %>%
  tibble::column_to_rownames("peakco")

  atac.glmtop.males.opening.fdr[fdr]

  write.table(
    atac.glmtop.males.fdr[fdr].opening[,c("chr", "start","end")],
    file = paste0("~/Downloads/sexM.age.groupHO.fdr",fdr, ".opening.bed"),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
  write.table(
    atac.glmtop.males.fdr[fdr].opening["GeneName"],
    file = "~/Downloads/sexM.age.groupHO.fdr"+ fdr + ".opening.GeneNames.txt",
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
}

```

# Advancing From `for` Loops to `apply`

As a beginner, `for` loops seem quite logical and easy to understand. However, R has a more powerful way to achieve the same thing through what is called "vectorisation" using the `apply` function or one of it's relatives `lapply`, `sapply` etc. Once you have understood and mastered `for` loops, you should take a look at `apply`, especially if you ever need to do a `for` loop many times or you find yourself doing nested `for` loops. `R` is notoriously slow with `for` loops and `apply` will likely be many times faster (potentially orders of magnitude). Here are some links for when you are ready to advance from `for` loops to `vectorisation`.

---
Tags: r
---


title: How to convert files to ADAM format?

text:
# Question
Title: How to convert files to ADAM format?
I would like to convert BAM and VCF files to ADAM format.

How do I do that?

# Answer
> 2 votes
These transformations are quite simple with the `adam-submit` script packaged with ADAM.

**Transformation of BAM files**

```
adam-submit -- transformAlignments sample.bam sample.alignments.adam

```

**Transformation of VCF files**

```
adam-submit -- transformVariants sample.vcf.gz sample.variants.adam

```

---
Tags: sam, vcf, adam
---


title: QMean standalone evaluation

text:
# Question
Title: QMean standalone evaluation
I want to use a script to evaluate a protein structure with QMEAN scoring.

I can't understand why the script suggested by SWISS-MODEL (access to QMEAN by json object) runs only if I put the commands on the shell. If I build a module, the result of the second request (request.get), return a QUEUEING status, not COMPLETE status.

Here is the code that I have been using:

```
import json
import requests

def qmean_page(structure):
""" return the Qmean4 score from SWISSMODEL QMEAN tool"""

    url = "https://swissmodel.expasy.org/qmean/submit/"

    # Request the resource
    response = requests.post(url=url, data={"email": "XXX@XXX.XX"} ,files={"structure": open(structure,'r'),},timeout=3,verify=True)

    json.dumps(response.json(), indent=4, sort_keys=True)
    download_json = str(response.json()['download_url'])

    return download_json

def qmean_score(download_url):
# Request the resource
    r = requests.get(download_url, verify=True)

    scores = json.dumps(r.json(), indent=4, sort_keys=True)
    print scores
    jdata = json.loads(scores)
    zscore = jdata['input_data']['models'][0]['global_scores']['qmean4']['zscore']
    return zscore

s = 'P02489.B99990005.pdb'
url_d = qmean_page(s)
print url_d
score  = qmean_score(url_d)
print score

```

**CASE A:** If I comment the last 5 line, I can import qmean.py (the function file) in python terminal.

```
>>> import qmean
>>> s = 'P02489.B99990005.pdb'
>>> url_d = qmean.qmean_page(s)
>>> print url_d
https://swissmodel.expasy.org/qmean/project/TgEjha.json
>>> score  = qmean.qmean_score(url_d)
{
"input_data": {
    "info": {
        "model_001": {
            "chain_mapping": {
                " ": "A"
            }
        }
    },
    "meta": {
        "created": "2018-01-18T10:04:29.978",
        "email": "XXX@XXX.XX",
        "project_name": "Project",
        "seqres_uploaded": false
    },
    "method": "QMEANDisCo",
    "models": [
        {
            "error_status": null,
            "global_scores": {
                "acc_agreement": {
                    "norm": 0.5606936416184971,
                    "zscore": -1.0459143647477258
                },
                "all_atom": {
                    "norm": -0.012626634534455788,
                    "zscore": -2.1825427852787604
                },
                "cbeta": {
                    "norm": -0.004260555683274807,
                    "zscore": -1.9157144117276892
                },
                "qmean4": {
                    "norm": 0.6275190794418422,
                    "zscore": -3.5006904416158258
                },
                "qmean6": {
                    "norm": 0.6339071011931049,
                    "zscore": -3.0507104571131918
                },
                "solvation": {
                    "norm": -0.609099435253941,
                    "zscore": -1.2499924754236265
                },
                "ss_agreement": {
                    "norm": 0.6706229586994028,
                    "zscore": 0.05371252044245716
                },
                "torsion": {
                    "norm": 0.0033300496230023145,
                    "zscore": -2.7188811736675764
                }
            },
            "model_pdb": "https://swissmodel.expasy.org/qmean/project/TgEjha/model_001.pdb",
            "modelid": "model_001",
            "name": "P02489.B99990005.pdb",
            "seqres": [
                {
                    "atomseq": "MDVTIQHPWFKRTLGPFYPSRLFDQFFGEGLFEYDLLPFLSSTISPYYRQSLFRTVLDSGISEVRSDRDKFVIFLDVKHFSPEDLTVKVQDDFVEIHGKHNERQDDHGYISREFHRRYRLPSNVDQSALSCSLSADGMLTFCGPKIQTGLDATHAERAIPVSREEKPTSAPSS",
                    "chain_name": " ",
                    "name": "seq_chain_0",
                    "sequence": "MDVTIQHPWFKRTLGPFYPSRLFDQFFGEGLFEYDLLPFLSSTISPYYRQSLFRTVLDSGISEVRSDRDKFVIFLDVKHFSPEDLTVKVQDDFVEIHGKHNERQDDHGYISREFHRRYRLPSNVDQSALSCSLSADGMLTFCGPKIQTGLDATHAERAIPVSREEKPTSAPSS"
                }
            ]
        }
    ],
    "options": {
        "qmeanbrane": false,
        "qmeandisco": true
    },
    "project_type": "default",
    "sequences": []
},
"status": "COMPLETED"
}
>>>

```

**CASE B:** If I uncomment the last 5 lines and run 'python qmean.py', it doesn't work.

```
  andrea$ python qmean.py
  https://swissmodel.expasy.org/qmean/project/zrcLCK.json
  {
  "status": "QUEUEING"
  }
  Traceback (most recent call last):
  File "qmean.py", line 31, in <module>
  score  = qmean_score(url_d)
  File "qmean.py", line 24, in qmean_score
  zscore = jdata['input_data']['models'][0]['global_scores']['qmean4']['zscore']
  KeyError: 'input_data'

```

If I print 'r' variable of 'qmean\_score' function, It returns:

```
  <Response [200]>
  {
  "status": "QUEUEING"
  }

```

# Answer
> 2 votes
So this is connected to how QMEAN is implemented as a server. It also could be true to many bioinformatics resources so I think this is important.

Posting to QMEAN "submit" triggers the job that is processed on our HPC cluster. That takes time and in the meantime JSON response just tells you that the job is QUEUEING/RUNNING or when it is ready it tells you COMPLETED. The answer suggested by you is the correct assuming you wish to wait for the job completion (which might take quite a lot).

So just check in while loop when the status is COMPLETED and fetch the result eg.

```
import time
import json
import requests
def qmean_score(download_url):
    # Request the resource
    zscore = None
    while True:
        r = requests.get(download_url, verify=True)

        scores = json.dumps(r.json(), indent=4, sort_keys=True)
        jdata = json.loads(scores)
        if jdata["status"] == "COMPLETED":
            zscore = jdata['input_data']['models'][0]['global_scores']['qmean4']['zscore']
            break
        else:
            time.sleep(1000)
    return zscore

```

---
Tags: python
---


title: Error in seq.default in chromPlot

text:
# Question
Title: Error in seq.default in chromPlot
I am using chromPlot to visualise the genome of *C. elegans*.

`library(chromPlot)`

I have created the following data frame with the lengths of *C. elegans* chromosomes.

```
Chrom Start      End    Name
1     1     0 15072434 contigs
2     2     0 15279421 contigs
3     3     0 13783801 contigs
4     4     0 17493829 contigs
5     6     0 20924180 contigs
6     7     0 17718942 contigs
7     5     0    13794 contigs

```

When I use `chromPlot(dataframe)`, I get a plot with 4 chromosomes and the following output/error:

```
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
   1 segment clusters.
Chrom 1 : 15072434 bp
Chrom 2 : 15279421 bp
Chrom 3 : 13783801 bp
Chrom 4 : 17493829 bp
Error in seq.default(minval, maxval, length = howmany) :
  'from' must be of length 1

```

Does anybody know why I get this error?

I don't get this error when I use the data frame provided by the package with human segments. The structure of the data frame is:

```
'data.frame':   457 obs. of  4 variables:
 $ Chrom: chr  "1" "1" "1" "1" ...
 $ Start: int  124535434 121535434 3845268 13219912 17125658 29878082 120697156 120936695 121485434 142731022 ...
 $ End  : int  142535434 124535434 3995268 13319912 17175658 30028082 120747156 121086695 121535434 142781022 ...
 $ Name : chr  "heterochromatin" "centromere" "contig" "contig" ...

```

Compared to my data frame:

```
'data.frame':   7 obs. of  4 variables:
 $ Chrom: chr  "1" "2" "3" "4" ...
 $ Start: int  0 0 0 0 0 0 0
 $ End  : int  15072434 15279421 13783801 17493829 20924180 17718942 13794
 $ Name : chr  "contigs" "contigs" "contigs" "contigs" ...

```

I am on MacOS X.

### Update

The 20th line of the package-provided hg\_gap data frame is:

```
   Chrom Start   End     Name
20     1     0 10000 telomere

```

Note that the Start is 0.

I've also tried changing the start from 0 to 1. The error stays the same.

### Update 2

Sorting the data frame by chromosome name does not improve the situation.

The traceback of the error is:

```
6: stop("'from' must be of length 1")
5: seq.default(minval, maxval, length = howmany)
4: seq(minval, maxval, length = howmany)
3: pretty_ticks(minlab, maxlab, 2, ...)
2: draw.scale(y = 0.8 * chr.length[[chrom]] + 0.2 * xylims[3], minval = margin,
       maxval = maxGeneCount, minlab = annot1_plot_range[1], maxlab = annot1_plot_range[2],
       lwd = 4, col = colAnnot1, cex = cex, title = scale.title)
1: chromPlot::chromPlot(gaps)

```

I have emailed the maintainer of the package Karen Orostica with a link to this post.

# Answer
> 2 votes
Using this data frame:

```
gaps
  Chrom Start      End    Name
1     I     0 15072434 contigs
2    II     0 15279421 contigs
3   III     0 13783801 contigs
4    IV     0 17493829 contigs
5     V     0 20924180 contigs
6     X     0 17718942 contigs
7 MtDNA     0    13794 contigs

```

and this command:

```
chromPlot(bands = gaps)

```

I could obtain this plot:

and then, if you want to add some segments for a subset of bases, you can do:

```
gaps2 <- gaps[1:2, ]
gaps2$End[2] <- 7000000

chromPlot(bands = gaps, segment = gaps2)

```

and obtain the following plot:

A more satisfactory answer could put segments of different colour on the chromosome, so bands but this is the best I've got so far.

---
Tags: r, bioconductor, visualization, chromplot
---


# Question
Title: Reduce number of transcripts in a highly variable de novo transcriptome assembly
I have a *de novo* assembly using both multiple SRA and locally sequenced transcriptomes. I started with 270M PE reads from 9 tissues. Here are the assembly stats generated with `TrinityStats.pl`:

```
################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':  543668
Total trinity transcripts:      1111089
Percent GC: 33.77

########################################
Stats based on ALL transcript contigs:
########################################

        Contig N10: 2117
        Contig N20: 1324
        Contig N30: 945
        Contig N40: 714
        Contig N50: 561

        Median contig length: 344
        Average contig: 496.17
        Total assembled bases: 551293315

#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

        Contig N10: 2172
        Contig N20: 1333
        Contig N30: 925
        Contig N40: 681
        Contig N50: 523

        Median contig length: 323
        Average contig: 473.63
        Total assembled bases: 257499729

```

As you can see, the total Trinity 'genes' number is very high. I suspect that this is due to a high biological variability of the transcriptomes (coming from organisms sampled in Italy, China, Spain...).

* How can I assess whether the problem is really the biological variability?
* How can I reduce this number? I'd rather not downsample the reads. It looks like CD-Hit might do what I'm looking for...

# Answer
> 3 votes
As you've suggested, CD-HIT works for reducing transcript numbers. We used a mixture of expression-based filtering and CD-HIT for reducing transcript counts for our genome-guided transcriptome assembly. This reduced numbers by a lot, without much change in BUSCO scores:

1. Map RNA-seq reads to Trinity-generated transcripts using Salmon
2. Use the expression of BUSCO genes to set a credible signal cutoff
3. Subset transcripts based on this threshold \[for us it was 50 counts\]
4. Identify the longest Met to Stop open reading frame for each transcript (may not be a good idea in all cases)
5. Run cdhit to collapse similar transcripts

Here's our CD-HIT command line:

```
cdhit -T 10 -c 0.98 -i longest_MetStopORF_HC50_TBNOCFED.fasta -o cdhit_0.98_LMOHC50_TBNOCFED.prot.fasta

```

More details here:

https://bmcbiol.biomedcentral.com/articles/10.1186/s12915-017-0473-4#Sec9

# Answer
> 0 votes
Seems to me that your "genes" are very short. What sort organism (prokaryote, lower eukaryote, higher eukaryote) is this? I would look to see if your "genes" contained ORFs, and only keep those with sufficicntly long CDS. You will miss uORFs and ORF fragments this way, but you'd have much more confidence in what you did find.

---
Tags: assembly, transcriptome
---


title: Negative scale-free topology

text:
# Question
Title: Negative scale-free topology
I am trying to build gene co-expression networks with WGCNA by combining in-house and publicly available RNA-seq datasets. I am interested in identifying gene networks associated with different environmental conditions. I used SVA to adjust for batch effects and "protected" different environmental conditions in the model as I don't want to lose information about them. However, when I used the "batch-adjusted" data to create networks, all the scale-free topology indices were negative (max = -0.05).

I tried generating the networks with softpower = 16 as it is suggested for sample sizes 20-30. However, only five networks were generated which is far less than I expected.

What might be the reason for getting negative scale-free toplogy indices?

Can combining studies from different environmental conditions cause the networks to be not biologically meaningful?

# Answer
Whenever you have batch or other nuisance variables driving variation you're going to get odd networks. You tried using SVA, which is a good method to use, but it seems in this case that it's able to optimally handle this dataset. I also strongly suggest that you use `pickSoftThreshold()` to choose a `softpower`, since that may produce better results (though since you already have negative values I'm not holding by breath). In general, you're going to need to play with either SVA or filtering to remove the batch-related genes that are causing problems.
> 1 votes

---
Tags: rna-seq, networks
---


title: How to download dbSNP database?

text:
# Question
Title: How to download dbSNP database?
I am bench-marking the performance of ADAM/Avocado on NGS whole-genome sequencing datasets. I am trying to repeat the tests reported in Frank Nofthat's (author of ADAM and Avocado) thesis, however, I am running into an issue acquiring the dbSNP database to use for Base Quality Score Recalibration.

The Base Quality Score Recalibration command (page 128-129):

```
adam-submit -- transformAlignments aln.adam bqsr.adam -recalibrate_base_qualities -known_snps dbsnp.adam

```

requires `dbsnp.adam`, which is described as "the dnSNP database" (converted to ADAM format) with reference to the NCBI database of genetic variation. This "database" can be found at this FTP site however I do not know which parts of these files to download and convert to ADAM in order to perform this step of the transformation pipeline.

Also, in order to use this genetic variation databaes do I need to make sure that the reference genome of "known variants" is the same genome that the sequencing data was aligned to?

Thank you.

# Answer
In the dbSNP FTP site you linked, you need to go into organisms and select your organism of interest (human obviously). You can then select which release of the SNP database you want (eg. 150 vs 151) \[the *modified* dates hint at when they were released\] and which genome build you want (eg. hg37 vs hg38).

Below I am looking at the *.vcf* files. I assume this what you want because they can be converted to ADAM format via the `adam vcf2adam` command.

For example human (newest genome and newest dbSNP release): ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/human\_9606\_b151\_GRCh38p7/VCF/

Your next decision is whether or not you want to use ALL the annotated SNPs in the database or just the common ones. For the former you choose the *All\_20180418* file and for the later the *common\_all\_20180418* file. I'm pretty sure the convention here is that the date in the filename is when it was uploaded/updated.

In your case I'm not sure what makes the most sense. Firstly, are you trying to replicate the initial conditions of when that software was published or do you want simply to benchmark the method with the best available data. Besides that, I imagine using all them is computationally much more taxing and if you are just coming up with general trend then using the common ones seems viable.
> 2 votes

---
Tags: adam
---


# Question
Title: How can I extract only insertions from a VCF file?
I'm looking to subset a standard VCF file to generate one which only includes insertions (i.e. not indels).

I can get part of the way there with:

```
bcftools view -v indels <vcf> | awk '{if(length($4) == 1) print}'

```

However this wouldn't catch an insertion that was part of a multi-allelic record with an insertion and deletion, where the reference length would be greater than 1 bp. Potentially then one way to go is a chain of decomposition, normalisation and then the same reference length filtering — surely there's a better way in one of the many VCF manipulation utils?

# Answer
> 4 votes
A one-liner:

```
zcat my.vcf.gz |
  perl -ane '$x=0;for $y (split(",",$F[4])){$x=1 if length($y)>length($F[3])}print if /^#/||$x'

```

or equivalently

```
zcat my.vcf.gz |
  perl -ane '$x=0;map{$x=1 if length>length($F[3])}split(",",$F[4]);print if /^#/||$x'

```

For simple VCF operations, I generally recommend to write a script. This may be faster than those using heavy libraries. With a script, you only parse fields you care about; most libraries unnecessarily parse every field.

---

On a related note, I recommend **not** to decompose multi-allelic sites unless necessary. Decomposing is tricky, makes VCF harder to parse and to understand and may be a potential source of errors. Here is an example:

```
#CHROM POS ID REF  ALT               QUAL FILTER  INFO  FORMAT  S1  S2  S3  S4
11     101 .  GCGT G,GCGA,GTGA,CCGT  199  PASS    .     GT      0/1 1/2 2/3 2/4

```

vt decompose+normalize produces the following VCF:

```
#CHROM POS ID REF  ALT QUAL FILTER INFO FORMAT  S1   S2   S3   S4
11     101 .  GCGT G   199  PASS   .    GT      0/1  1/.  ./.  ./.
11     101 .  G    C   199  PASS   .    GT      0/.  ./.  ./.  ./1
11     102 .  CGT  TGA 199  PASS   .    GT      0/.  ./.  ./1  ./.
11     104 .  T    A   199  PASS   .    GT      0/.  ./1  1/.  1/.

```

In theory, you can reconstruct the original VCF from this output. However, it is very challenging for a program to do that. When you compute allele frequency line-by-line, this VCF will give you wrong results. `bcftools norm -m-` replaces "." with "0". You can get a correct ALT allele frequency from the bcftools output, but a wrong REF allele frequency. Furthermore, vt is also imperfect in that "CGT=\>TGA" is not decomposed.

My preferred output is:

```
#CHROM POS    ID     REF    ALT    QUAL   FILTER INFO   FORMAT S1     S2     S3     S4
11     101    .      GCGT   G,<M>  0      .      .      GT     0/1    1/2    2/2    2/2
11     101    .      G      C,<M>  0      .      .      GT     0/2    2/0    0/0    0/1
11     102    .      C      T,<M>  0      .      .      GT     0/2    2/0    0/1    0/0
11     104    .      T      A,<M>  0      .      .      GT     0/2    2/1    1/1    1/0

```

Here we use a symbolic allele `<M>` to represent "another ALT allele". You can calculate the allele frequency by looking at one line, and won't confuse other ALT alleles with REF. bgt can produce such a VCF indirectly. However, it discards all INFO, so is not a practical solution, either.

In summary, it is very difficult to decompose multi-allelic sites. When you get decomposing wrong, your downstream analyses may be inaccurate. Decomposition should be used with caution.

# Answer
> 5 votes
Just to highlight that all the steps can be done within bcftools capabilities, and since I can't just comment on @blmoore 's answer:

```
bcftools view --types indels <vcf> |
  bcftools norm -m - |
  bcftools filter --include 'strlen(REF)<strlen(ALT)' |
  bcftools view -H

```

more bult-in functions for 'bcftools filter' here

# Answer
> 3 votes
One method is to decompose multi-allelic records so that they're represented as one-allele, one-record using vt or similar:

```
bcftools view -v indels <vcf> |
  vt decompose - |
  bcftools view -H |
  awk '{if(length($5)>length($4)) print}'

```

Some decomposed alleles will come with excess reference padding. To left-shift and trim these, add a `normalize` step (NB matching these back to your input VCF becomes non-trivial):

```
bcftools view -v indels <vcf> |
  vt decompose - |
  vt normalize -r <reference.fasta - |
  awk '{if(length($4)==1) print}'

```

**edit:** As gringer suggests, this can also be done without vt:

```
bcftools view -Ou -v indels <vcf> |
  bcftools norm -Ou -Nm - |
  bcftools view -H |
  awk '{if(length($5)>length($4)) print}'

```

To also include complex alleles, use `view -V snps` (!snvs) instead of `-v indels`

# Answer
> 2 votes
Using vcffilterjs

* get the length of the REF;
* loop over the ALT, ignore the symbolic
* accept the variant if it's an insertion , eq: len(ALT)\>len(REF)

.

```
java -jar dist/vcffilterjs.jar -e 'function accept(vc){var a=vc.getAlleles();var lenRef=a.get(0).length();for(i=1;i<a.size();++i) {var alt=a.get(i);if(alt.isSymbolic()) continue;var lenAlt=alt.length(); if(lenRef<lenAlt) return true; } return false; }accept(variant);' input.vcf

```

.

# Answer
> 0 votes
If you're doing set operations, you could use `vcf2bed`:

```
$ vcf2bed --insertions < in.vcf > out.bed

```

Based on the VCF v4.2 specification, `--snvs`, `--insertions`, and `--deletions` are options available to filter input. In each case, the length of the reference and alternate alleles is used to determine which type of variant is being handled.

---
Tags: vcf
---


# Question
Title: Why is FPKM still used for gene expression studies?
Perhaps my understanding is wildly misguided, but I'm seriously confused about why people are still interested using FPKM values for cross sample gene expression analysis. My understanding is, that for some gene ${ i }$, its $FPKM$ is calculated as follows:

$$ FPKM\_i \propto \frac{ \mathrm{fragments\ mapping \ to\ gene \ }i } { \sum\_j \mathrm{ fragments\ mapping\ to \ gene\ } j } \bigg/ \mathrm{effective\\\_length}(i) $$

For comparing gene expression within a single sample, this doesn't seem too harmful to me. In fact, this seems like a good idea.

However, the use of $FPKM$ or even the famous $log\_2(FPKM+1)$ in a similarity analysis seems to me plain wrong. True differences in RNA abundances of a couple of genes would significantly impact the FPKM values of the other genes leaving perhaps confused about the true differences between the samples.

Am I missing something here? Is it ever meaningful to use FPKM values for analyzing across samples?

# Answer
> 1 votes
> Is it ever meaningful to use FPKM values for analyzing across samples?

One should never use FPKMs for anything important. They can occasionally be useful for plotting, but even in that case one needs to construct the FPKMs from properly normalized data and not use the original definition.

It's becoming ever rarer to see people using FPKMs these days, thankfully.

---
Tags: rna-seq, fpkm
---


# Question
Title: Calculating L2 norm (Euclidean norm) to set of genes in time series
I have gene expression RNA-seq data of 9 time points (0 hour, 2 hour, ..., 16 hour). I need a plot on which x axis is time point and each data point at a certain time point represents the L2 norm (Euclidean norm) of the difference between the gene expression vector at that time point and the previous time point.

This is the link of my expression data

https://drive.google.com/file/d/1crJSqVdGujy1tFyeoV6E7e-xO9SXteJZ/view?usp=sharing

# Answer
If you read your matrix into R and saved it as `m`:

```
diffMatrix = m[,2:9] - m[,1:8]  # The change per-timepoint
l2n = apply(diffMatrix, 2, function(x) sqrt(sum(x^2)))  # The L2-norm

```

R also has a `norm()` function, which can directly produce the L2-norm, but it turns out that this is slower than the simple method shown here. If you click through to the other question you can find somewhat faster implementations of this, but anything will suffice for small datasets like yours.
> 1 votes

---
Tags: r
---


# Question
Title: Is there a way to tell which chromosome a gene is on, by looking at the "Chromosome/scaffold name"
I recently got a data set, from which I need to figure out which chromosome a gene is from, but the head of the data reads like:

```
Gene ID Description     Gene type       Gene End (bp)   Gene Start (bp) Strand  Associated Gene Name    Chromosome/scaffold name

```

while one entry of the data is like:

```
ENSG00000252303 RNA, U6 small nuclear 280, pseudogene [Source:HGNC Symbol;Acc:HGNC:47243]       snRNA   67546754        67546651        1       RNU6-280P       CHR_HG2128_PATCH

```

It looks like the last column is the information about which chromosome the gene comes from, but how to interpret the "CHR\_HG2128\_PATCH"?

# Answer
What you are looking at is a patch, not an actual chromosome location. HG-2128 is an issue ID, something that the Genome Reference Consortium uses to track issue with references. Information regarding that specific patch is here. The issue is described as a "Possible misassembly or indel variation in GRCh38 within AL049860.8."

During the alignment step before you received your data, someone used a top level assembly version of the reference genome. This includes extra, accessioned scaffolds tacked on to the end of the file that contain the fix. These scaffolds represent what will be found in the next major release of the genome. The fixes aren't applied directly to the reference as soon as they are found due to the inevitable disruption of coordinates. So reads mapped to this FIX patch rather than the primary original scaffold in the primary assembly, on which the RNU6-280P gene resides.

The coordinates are slightly different. In the current GRCh38.p12 assembly the RNU6-280P gene is located at chr6:67546651-67546754 (source). You can see that the coordinates are very close to those indicated in the table you provided (Gene Start, Gene End) but are indeed different.

Here is some more information on patches.
> 2 votes

---
Tags: gene, ensembl
---


# Question
Title: Increase number of threads for GATK 4.0 HaplotypeCaller
I am using GATK version 4.0, I tried to use multiple threads for calling variants using HaplotypeCaller using following command

```
gatk --java-options -Xmx90G  -nt 28 HaplotypeCaller -I output.bam -R wheat_ref.fa -O final.vcf

```

and the error is

```
'-nt' is not a valid command

```

I have also used **--num\_threads** but again the error was the same, as written here this parameter is available in previous versions of GATK but I don't know what's its status in current version of GATK i.e GATK 4.0. Whether it is deprecated or replaced by some other parameter?

Is there any way of calling multiple threads for variant calling with GATK 4.0?

# Answer
I am assuming you are using the non-SPARK version of the method. To specify the number of threads you wish to use with HaplotypeCaller, include `--native-pair-hmm-threads` (documentation). This will only parallelize the pair hidden Markov models (pair HMM) process. This is an argument for the method HaplotypeCaller itself, so it should be located after the method is called as opposed to one of the `--java-options`.

I do recommend trying out the SPARK implementations of the GATK tools. The whole method will take advantage of the threads you specify. I've been running SPARK locally (i.e. not on Firecloud) on computational servers. I've been impressed with the enhanced performance.
> 3 votes

---
Tags: sequence-alignment, variant-calling, gatk
---


# Question
Title: Genotyping with Illumina HumanOmni1-Quad before whole genome sequencing
I read a paper about (whole genome sequencing analysis of 100 Southeast Asian Malays. It says "Prior to library preparation, each sample was genotyped on the Illumina HumanOmni1-Quad as an initial DNA QC and for in-house assessment of gender and familial relations". I can hardly understand how Illumina HumanOmni1-Quad worked here for initial DNA QC and for in-house assessment of gender and familial relations. Could you please kindly explain to me a bit about this sentence?

Many thank in advance

# Answer
> 2 votes
The Illumina HumanOmni1-Quad beadchip is a microarray device consisting of 1,140,419 markers which have derived from the 1,000 genomes project. The markers chosen are is high-value regions of the genome: ABO blood typing SNPs, cSNPs, disease-associated SNPs, eSNPs, SNPs in mRNA splice sites, ADME genes, AIMs, HLA complexes, indels, introns, MHC regions, miRNA binding sites, mitochondrial DNA, PAR, promoter regions, and Y-chromosome.

The markers come together to form a sort of individual fingerprint. Germline SNPs are heritable and so can be traced up the line. Some SNPs are more prominent in subpopulations, something that Ancestry.com and other similar services leverage to predict ethnicity.

This paper utilized genomic DNA extracted from the Singapore BioBank. BioBanks tend to record age, lifestyle, family history, and health information alongside the genomic DNA. I (conservatively) assume that the Singapore BioBank at least has familial relations and gender information tied to the DNA. My interpretation is that the ivestigators SNP genotyped using the Illumina Beadchip the genomic DNA from the bank and compared their results with the "ground truth" provided by the genetic bank.

Illumina brochure for the HumanOmni1-Quad beadchip.

---
Tags: illumina, genome-sequencing
---


# Question
Title: Output of Seurat FindAllMarkers parameters
I compared two manually defined clusters using `Seurat` package function `FindAllMarkers` and got the output:

Now, I am confused about three things:

1. What are `pct.1` and `pct.2`?
2. How come `p-adjusted values` equal to `1`? What does it mean?
3. If we take first row, what does `avg_logFC` value of `-1.35264` mean when we have cluster `0` in the `cluster` column? Is it that in cluster `0` the `Cttnbp2` gene downregulated by a factor of `2^1.35264`?

# Answer
> 3 votes
1. `pct.1`– The percentage of cells where the gene is detected in the first group
2. `p_val_adj`– Adjusted p-value, based on bonferroni correction using all genes in the dataset.

> This is **not** also known as a false discovery rate (FDR) adjusted p-value.

An adjusted p-value of 1.00 means that after correcting for multiple testing, there is a 100% chance that the result (the logFC here) is due to chance.

3. You are correct. The expression of `Cttnbp2` donwn-regulated by that factor. If you were to look in the table at the cluster 1 DE genes you would see `Cttnbp2` with a logFC value of 1.35264. When you are looking at genes differentailly expressed in only two clusters, the sign is the only thing that changes.

Vignette posted above by @heathobrien.

---
Tags: r, seurat
---


title: How to filter for "CAF=0.5747,0.4253"

text:
# Question
Title: How to filter for "CAF=0.5747,0.4253"
In my gzipped VCF file I have many variants, and they report allele frequency such as:

```
##INFO=<ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">

```

and

```
CAF=0.5747,0.4253

```

I would like to filter out the variants \>= 1% allelle frequency. How do I do that?

# Answer
using **vcffilterjdk**:

you want to keep the variant if there is no INFO/CAF.

Otherwise, get All values of CAF as a String list (because it is declared as a string in the VCF header...), stream those data, skip the first value for REF , convert each string to double, you want to keep the variant if any allele is rare (\< 0.01)

```
java -jar dist/vcffilterjdk.jar -e 'return !variant.hasAttribute("CAF") || variant.variant.getAttributeAsStringList("CAF","1.0").stream().skip(1L).map(A->Double.parseDouble(A)).anyMatch(A->A<0.01);' input.vcf

```
> 4 votes

# Answer
This could be done with `bcftools view`. To include variants where the first element in CAF is \> *x*:

```
bcftools view -i 'INFO/CAF[0] > 0.1' ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/common_all_20180418.vcf.gz

```

However, it appears that CAF is encoded as a list of strings in this input so the `>` operator cannot be used. In addition, you may want more control on the filtering. In this case, I would suggest using a dedicated library in your favourite language for dealing with VCF files. In python, you could use pysam along these lines:

```
python -c "import pysam
invcf= pysam.VariantFile('ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/common_all_20180418.vcf.gz')

# Write to stdout
outvcf= pysam.VariantFile('-', mode= 'w', header= invcf.header)

for line in invcf:
    if float(line.info['CAF'][0]) > 0.1:
        outvcf.write(line)
outvcf.close()
invcf.close()
"

```

I would stay away from parsing VCF files with regex, awk, or custom text processing scripts as even if the format is quite well defined they can be quite tricky.
> 0 votes

---
Tags: vcf, file-formats
---


# Question
Title: Difference between samtools mark duplicates and samtools remove duplicates?
What is difference between samtools mark duplicates and remove duplicates ? Is it necessary to mark duplicates before removing duplicates with samtools?

# Answer
> 3 votes
`samtools rmdup` and `samtools markdup -r` do the same thing. Without the `-r` flag `samtools markdup` only flags the duplicates.

You'll have to run `samtools fixmate -m` and sort the output to add ms and MC tags to prepare the file for `markdup`..

According to the current documentation, `rmdup` is obsolete. Please do not use `rmdup`.

# Answer
> 2 votes
**TL;DR:** just use `markdup`.

`rmdup` removes duplicates from BAM, while `markdup`, like Picard's `MarkDuplicates`, marks duplicates by default without hard removal – the latter is usually the desired behavior. In addition, `markdup` implements a better algorithm that takes care of more corner cases and gives more consistent results.

PS: to emphasize – I see no point to use `samtools rmdup` nowadays. It has been declared obsolete as @kohlkopf said.

---
Tags: sequence-alignment, ngs, samtools
---


# Question
Title: Batch detection of CRISP proteins in fasta file
Probably a naive question. I am inexperienced.

I am interested in identifying potential CRISP (Cysteine-rich secretory proteins) in a certain tissue transcriptome (ca. 20k sequences in fasta). I have detected signalP and estimated % of cysteine in sequences. However CRISP rely on a certain pattern, and I am not aware of batch-detected algorithms.

Please would anyone know the best/simplest way of detecting CRISP candidates?

\[EDIT\]

I have been asked for the pattern. I summarise what the literature says:

* Several of them have a CAP domain;
* There should be a region called 'cysteine-rich domain' (CRD) upon the carboxyl terminal half of the protein containing 10 of roughly 16 conserved cysteines.

I am particularly interested in venom-derived CRISPs, but I am not sure whether they cluster in a certain CRISP family.

References: Roberts et al. (2007) Structure and function of epididymal protein cysteine‐rich secretory protein‐1. Yamazaki & Morita (2004) Structure and function of snake venom cysteine-rich secretory proteins.

# Answer
> 1 votes
Is your goal to detect the actual proteins from the sequences? If you'll settle for detection of transcripts from the genes encoding these proteins as opposed to detection of the proteins themselves, here is an approach.

Aligning your transcriptome to a reference genome (probably with STAR) will be the first step. You have lots of short reads in your FASTA, but in their current state they are no more informative than a box of puzzle pieces. Aligning your reads to the reference will allow you to extrapolate which genes that you short reads correspond with. From here you can quantify the features (genes) in your transcriptome.

Lucky for you, there are only 3 genes in the cysteine rich secretory protein family, `CRISP1`, `CRISP2`, and `CRISP3`. Alignment with a reference (and annotations) and quantification of genes that your reads align to will give you all you need. You'll just have to find the three genes in your final counts table.

I recommend STAR as an aligner with the `--quantMode GeneCounts` argument. It's very quick if you've got the resources and is geared towards RNA-seq reads.

---
Tags: sequence-annotation, blast, transcriptome
---


# Question
Title: GRCm37-designed exon target enrichment, which reference to use?
I have exomes from 24 individual mice. The exomes are the product of Roche's SeqCap EZ HyperCap target enrichment kit. I see that the mouse exome design comes from mm9/NCBI37, the previous major reference genome released in July 2007.

Which reference should I use to align these reads? What artifacts may arise from mapping these reads designed from mm9 to the mm10 reference? I have already aligned the reads to mm10 and noticed an erroneous exon right in the middle of an intron, perhaps a target that was an exon in the mm9 release but revised out of the mm10 release. My goal is to investigate structural variation.

# Answer
The only possible issue with using GRCm38 as the reference is that since GRCm37 was used in the bait design there are likely a few probes that will show artificial structural variations, due to errors in GRCm37. However, these will be present in all of your mice and if you use liftOver to look at the regions around the probes then such regions should become obvious. In fact, that'd be a convenient thing to do if you want to keep your (likely better) GRCm38 alignments, flank each probe site by a few hundred bases and then use liftOver to get the GRCm37 coordinates. If you go through the result you'll be able to find the false-positive structural variants due to the baits being designed with an older reference.

Alternatively, just use GRCm37 and hope that any issues with it don't result in false-negative SV calls.
> 2 votes

---
Tags: sequence-alignment, variant-calling, reference-genome, exome
---


# Question
Title: pysam pileup: what reads appear in the pileup?
Cross-posted on Biostars (with no answers currently). I hope that's OK.

---

I thought that ***iff*** a read matches the reference at a position, it would appear in the pileup column for that position. But tracking a single read in pysam shows that that's not the case.

This read appears in some pileup columns for positions to which it does not map, and does not appear in some pileup columns to which it does map.

What's wrong in my understanding?

---

Updated code (and updated output) moved to Gist.

---

**Edit**: Here's the raw SAM read. I wasn't separating the read from it's mate when I calculated the figures above.

```
E00489:44:HNNVYCCXX:1:2120:26524:19135  163     5       345248  0       48M102S =       345446  350     TGTGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTGTGGGCGGGAGTGCGTGGGGGAGGGTGAGAGTGTGGATGTCGGTGGTCGCGGCTGCGTGAGAGAGGGGGTGTGTGGGGGGGGGGGGGGGGGGGGGTGTGGGTGGGCG  AAAFFJFJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJ-J-7-777-7-7-7--7-77-77F---7-7AFJ7----7-7FAA7-7---7----77-777------7-7)7-AAA-A)7FA)<)<--)<)--<-FAJFJFFFF-<  XA:Z:5,+345312,48M102S,0;5,+345554,48M102S,1;5,+345452,4S44M102S,0;5,+345192,48M102S,1; MC:Z:6M2D144M   MD:Z:48 RG:Z:HNNVYCCXX.1       NM:i:0   MQ:i:60 AS:i:48 XS:i:48
E00489:44:HNNVYCCXX:1:2120:26524:19135  83      5       345446  60      6M2D144M        =       345248  -350    TGGGGATGTGTGTGTGTGTGTCGGATGATGTCCCTGGTTGTGTGTGGGGATGTGTGTGTGTGTGGGTTGATGGTCCTGGCTGTGTGTGTGTGTGGGTGTGCGTGTGTGGGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTG  JJAF7-FJJJFJFF-AFJF7--<A--<F<7-7--FA--F-J--A-7--7-FJFJ7JFJJF<JA77-7----7----<F7-<JJF-<<FAJJJJA-7-<JAJJ<7FFJJ<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA  MC:Z:48M102S    MD:Z:6^TG31C25C2A5T0C76 RG:Z:HNNVYCCXX.1        NM:i:7  MQ:i:0  AS:i:119        XS:i:69

```

I guess this explains why the read appears in more than `150` pileup columns—in some of the columns, it's really its mate. But given the CIGAR strings for these paired reads, I would expect one of these mates to appear in `198` = (`48` \+ `6` \+ `144`) pileup columns, at each position where one of them maps?

---

**Update:** I re-ran this script, using @Bioathlete's suggestion of using the flags to distinguish the read and its mate. In case `pileup()` was filtering reads with low base quality, I also calculated how many positions the read has `baseQ >= 13`. It's `150` positions, but this read still appears in pileup columns for only `121` positions.

# Answer
The two entries in the sam file represent mate pairs for the given ID. You can tell the difference based on the sam flags. Picard has a nice tool to determine the meaning of the flags.

You are missing coverage within the first read. I think that it is related to the low quality bases. samtools mpileup filters bases with a quality below 13 and pysam just wraps the samtools C functions.
> 2 votes

---
Tags: python, sam, sequence-alignment, pysam
---


# Question
Title: plotting two heatmaps with the same order of genes
I expected the same pattern but here I am not able to compare the patterns as the order of genes does not seem the same. I mean how I can have two heat maps on which the order of genes are the same so I woulds be able to compare the block of yellow and blue colors in two heat maps.

I have used this code

```
require("RColorBrewer")
myCol <- colorRampPalette(c("dodgerblue", "black", "yellow"))(100)
myBreaks <- seq(-2, 2, length.out=101)
heat <- t(scale(t(sc_DEGG)))
hr <- hclust(as.dist(1-cor(t(y), method="pearson")), method="complete")
heatmap.2(
    heat,
    Rowv=as.dendrogram(hr),
    Colv=as.dendrogram(hc),
    col=myCol,
    breaks=myBreaks,
    main="Title",
    key=T,
    keysize=1.0,
    scale="none",
    density.info="none",
    reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean),
    trace="none",
    cexRow=0.2,
    cexCol=0.8,
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

```

I have made a pair of heat maps based on your code in complexheatmap package. The right heat map is lighter than the left one. I don't know which one has been clustered based on the another one. How I can make the right one darker so that both seem more similar??

# Answer
> 11 votes
ComplexHeatmap is built for plotting side-by-side heat maps with the same clustering - you use the `+` notation, similar to ggplot2. To use the same example data as @b.nota:

```
library(ComplexHeatmap)
# First matrix
set.seed(2)
m <- replicate(8, rnorm(26))
rownames(m) <- letters[1:26]

# Second matrix
set.seed(3)
m2 <- replicate(8, rnorm(26))
rownames(m2) <- letters[1:26]

Heatmap(m) + Heatmap(m2)

```

And the results:

There are lots of examples in the extensive documentation.

EDIT:

Since you have a few genes with high values in the 'T' matrix, they are washing out the colour scheme for the rest of the heatmap. You could set the max to be lower, and accept that those genes will saturate - the resulting heatmap will be less 'honest' about the data, but might capture what you want to see. To change colour scheme requires a call to the `colorRamp2` function of the `circlize` package.

```
library(circlize)
mycol <- colorRamp2(c(-2,0,2), c("dodgerblue", "black", "yellow"))
Heatmap(m, col=mycol, cluster_columns = FALSE) +
    Heatmap(m2, col=mycol, cluster_columns = FALSE)

```

This gives:

# Answer
> 6 votes
You expect that the order of genes is the same for both heatmaps, but the chance is small that the 4000 genes will exactly be clustered in the same order.

What you can do, is use the order of genes (rows) from one heatmap and create a new heatmap with the other data set. You'll lose the tree of course because it will not be clustering anymore.

Here a small example:

```
library(gplots)

# First matrix
set.seed(2)
m <- replicate(8, rnorm(26))

rownames(m) <- letters[1:26]

hm <- heatmap.2(
    m,
    main="m",
    key=T,
    keysize=1.0,
    scale="none",
    trace="none",
    density.info="none",
    Colv = FALSE,
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

# Second matrix
set.seed(3)
m2 <- replicate(8, rnorm(26))

rownames(m2) <- letters[1:26]

hm2 <- heatmap.2(
    m2,
    main="m2",
    key=T,
    keysize=1.0,
    scale="none",
    trace="none",
    density.info="none",
    Colv = FALSE,
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

# Second matrix with row order of first
reorder_m2 <- m2[rev(hm$rowInd),]

hm2_m <- heatmap.2(
    reorder_m2,
    main="m2 in m's order",
    key=T,
    keysize=1.0,
    scale="none",
    trace="none",
    density.info="none",
    Rowv = FALSE,
    Colv = FALSE,
    dendrogram="none",
    distfun=function(x) dist(x, method="euclidean"),
    hclustfun=function(x) hclust(x, method="ward.D2")
)

```

# Answer
> 2 votes
You have 2 options provided that both datasets have the same genes (rows):

1. Keep the original order of rows in the heatmap. If they’re in the same order in the dataset then they will be in the heatmap. You can suppress reordering with any of the following arguments to `heatmap.2`

`Rowv = NULL`

`dendrogram = “col”` or `dendrogram = “none”`

`reorderfun = function(x) return(x)`

2. You can compute the order of the genes (rows) in advance and pass them to the heatmap. The code shown in the question does this already with `Rowv=as.dendrogram(hr)`. To have the same order in a second heatmap, you have to pass the same same dendrogram `hr` to both heatmaps. This can be computed from either dataset or another set of samples altogether.

Another option is to compute separate dendrograms (from `hclust` of each dataset) for the samples and combine them using `merge`. You can combine the datasets with `cbind` and plot the datasets together as one heatmap (with 2 sets of samples clustered separately) by passing the merged dendrogram to `Colv`:

```
#compute same dendrograms separately
hc1 <- hclust(as.dist(t(data1), method="euclidean")), method="complete")
hr2 <- hclust(as.dist(t(data2), method="euclidean")), method="complete")
mergedcols <- merge(as.dendrogram(hr1), as.dendrogram(hr2))
mergeddata <- cbind(data1, data2)
#compute gene dendrogram across all samples
hr <- hclust(as.dist(1-cor(t(mergeddata), method="pearson")), method="complete")
heatmap.2(
mergeddata,
Rowv=as.dendrogram(hr),
Colv=as.dendrogram(mergedcols),
col=myCol,
breaks=myBreaks,
main="Title",
key=T,
keysize=1.0,
scale="none",
density.info="none",
reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean),
trace="none",
cexRow=0.2,
cexCol=0.8,
distfun=function(x) dist(x, method="euclidean"),
hclustfun=function(x) hclust(x, method="ward.D2")
)

```

You should be able to show the division between these samples with the `colsep` argument.

---
Tags: r, heatmap
---


title: Find paralogs in a draft genome

text:
# Question
Title: Find paralogs in a draft genome
We generated a (diploid, chordata, highly heterozgous) genome using PacBio and we wanted to see whether it contains lineage-specific duplications (paralogs, basically). The genome is not in Ensembl yet.

The only data we have at the moment are:

* genome
* transcript annotation
* RNAseq

We found some methods from papers:

* use Blast
* detect segmental duplications in complete genomes with SDDdetector
* detecting putative recent segmental duplications or diploid homozygous organisms based on NGS data DuplicationDetector
* (I just came up with this, but should work) map back the reads to the assembly & analyze the read depth to detect duplicated portions

I'll gladly take advice.

# Answer
> 1 votes
This is really tough to do with highly heterozygous animals. What are your genome assembly stats? Specifically, what is your number of contigs, scaffolds, assembly size, and the N50?

If you have proximity ligation data it will be easier to determine if potential paralogs are truly from different regions of the genome or are just from homologous locations on sister chromatids that ended up in uncollapsed in your final assembly.

If I were you and had a good genome, I would start by annotating gene models using RNAseq data, then blasting each gene in the model against all genes in the model. This would give you potential paralogs to start looking at more closely.

---
Tags: genome, assembly, duplications, paralogs
---


title: Associating SNP and GENE

text:
# Question
Title: Associating SNP and GENE
Assuming I have SNPs data using hg19, how can I know which SNP belongs into which Gene?

The data looks like:

```
chr10_103577643
chr10_124712463

```

and so on.

I want to add a column of Gene, which would tell to which Gene the SNP belongs. The file is not a VCF file but a Txt file.

# Answer
> 3 votes
Use the Variant Effect Predictor(VEP)

# Answer
> 2 votes
You can use snpEff or snpSIFT java annotator also for hg19.

---
Tags: gene, snp, hgp
---


# Question
Title: How many reads do I need to cover the entire genome?
Suppose my genome is 3 million bases and that my reads are 100 nucleotide long. I need to know how many reads I need to cover the entire genome.

I start from using the equation $C = \frac{N \cdot L}{G}$ where C is the coverage, N the number of reads, L the length of a read and G the length of the haploid genome. I also know that the average coverage follows a Poisson distribution, so I need this for C.

On my slides, in fact, I have the following:

$$ C = \frac{N \cdot L}{G} \approx ln \left( \frac{G}{L \cdot \epsilon}\right), $$

if we assume that this is the average coverage needed to cover the entire genome with probability $1-\epsilon$.

The problem is that I do not understand how to arrive at the expression on the right. I thought that if we want to be sure to have a coverage of at least 1 on average, we would have had to compute $1 - P(X=0) = 1 - e^{-\frac{N \cdot L}{G}} = 0.99$, if we take $\epsilon=0.01$. But I do not get the same result as with the equation above.

# Answer
Let's see what your slides are claiming $\epsilon$ to be, using:

$$ C = \frac{NL}{G} \approx \ln \left( \frac{G}{L\epsilon} \right) $$

We can rearrange to get

$$ \epsilon \approx \frac{G}{L} \times e^{-C} $$

Now, what does it mean for an entire genome to be covered? Having a clear understanding of this definition is crucial, and your misunderstanding seems to partially stem from this.

I claim it's natural to consider a genome as consisting of $G/L$ bins. In this setting, we have $N$ balls that each randomly goes into a bin. Then, the entire genome is covered when all bins contain at least $1$ ball/read.

From the Poisson distribution, we know that P(a bin has 0 balls) = $e^{-C}$.

Then, we define $\epsilon$ to be the probability that at least one of the $G/L$ bins has 0 balls. Given the way we've modeled the problem, this is exactly equal to

$$ \epsilon = 1 - (1 - e^{-C})^{G/L}. $$

From above, your slides are claiming $ \epsilon \approx \frac{G}{L} \times e^{-C}$. This might be strange, as it appears to be equal to the expected number of bins with 0 balls, instead of the probability that there is at least 1 bin with 0 balls. Furthermore, it's strange because it contains only one exponential, compared to our exact derivation which has a double exponential.

The mismatch in the number of exponential is a hint that another mathematical approximation is at play. We can use a first-order approximation of the binomial series to obtain

$$ (1-e^{-C})^{G/L} \approx 1 - (G/L)e^{-C} $$ which holds particularly when $|e^{-C}| \times \frac{G}{L} \<\< 1$.

In your case, $\frac{G}{L} = 30,000$, so the approximation only really starts getting accurate when $C \> 15$ or so.

Following through with our approximation, we get

$$ \epsilon \approx 1 - (1 - (G/L) e^{-C}) $$

$$ \epsilon \approx \frac{G}{L} \times e^{-C}. $$
> 3 votes

---
Tags: ngs, coverage
---


# Question
Title: RNA seq fasta file annotation from alignment to reference matches
I've got a fasta file with some RNA seq data and another csv file with the output from plast where I've aligned it to a reference using plastn. I'm struggling with figuring out a command to append my original RNA seq data fasta headers with the annotations.

I've looked at awk and sed but can't wrap my head around how to use them in this instance.

CSV output example:

```
Contig11    ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA   78.95   95  20  0   371 465 4627    4721    1e-14   82.4    75  625 1   0   15.2    0   5692    1   0   15.2    0

```

And example fasta contig:

```
>Contig11
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA

```

How would I be able to parse through the csv to find a matching contig with the fasta file, and append the second csv column to the fasta header. Ideally to produce:

```
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA....

```

Cheers for any help!

# Answer
> 1 votes
An `awk` approach (assuming your input "csv" file is space-separated):

```
$ awk '(NR==FNR){a[">"$1]=$2; next;} (/>/){$1=$1" "a[$1]}1;' foo.csv foo.fa
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA

```

### Explanation

* `NR==FNR` : `FNR` is the line number of the current input file and `NR` is the line number of all input. When processing more than one file, the two will only be equal while reading the first file.
* `(NR==FNR){a[">"$1]=$2; next;}` : if we're reading the first file, the csv, save the second field as the value in the associative array `a`, whose indices are a `>` followed by the 1st field. I use the `>` only to simplify matching it later in the fasta file. The `next` moves on to the next line so the rest of the command isn't run for the 1st file.
* `(/>/){$1=$1" "a[$1]}`: Now we're in the 2nd file. If the line starts with a `>`, then set the 1st field to be itself, a space and then whatever was stored in the array `a` for this 1st field.
* `1;` : this is awk shorthand for "print this line".

You could expand the above into:

```
awk '{
        if(NR==FNR){
            a[">"$1]=$2;
        }
        else{
            if(/>/){
                $1=$1" "a[$1]
            }
        }
        print;
      }' foo.csv foo.fa

```

---

If your files are small enough to be sorted easily, you can use the `FastaToTbl` and `TblToFasta` scripts I have posted before and then `join`:

```
$ join -o 1.1,2.2,1.2 <(FastaToTbl foo.fa | sort ) <(sort  foo.csv ) | TblToFasta
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
>Contig12 mamamama|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
AKHSALKJHDLKAJBDLKAJDBADADKJALDKJBAD

```

---
Tags: rna-seq, sequence-annotation, awk, functional-annotation
---


# Question
Title: Bioinformatics approach to dentifing potential PCR primer sequences for transcribed gene
I have an annotated transcriptome and would like to develop PCR primers for particular transcribed genes. My species is a non-model plant. Can I use BLAST or another tool to identify potential PCR primer sequences? Or more generally, is there a workable bioinformatics approach to identifying potential PCR primers for transcripts arising from a non-model organism? Any literature refs would be great too.

# Answer
It depends if you want to amplify your target from genomic DNA or cDNA.

Depending on intron size, primers that anneal to different exons may not be able to amplify across the intron(s), and if you're unlucky enough to design a primer that spans a splice junction, you'll definitely be stuffed if you're working with genomic DNA. I would try blastn or translate to an amino acid sequence and use tblastn to see if you can identify the intron/exon structure in a related species. If you don't have a genome sequence available for anything close, I think your only option is to target as small an amplicon as possible and hope for the best.

If you are working with cDNA, any standard primer design tool should do
> 2 votes

---
Tags: blast
---


# Question
Title: Selection of differential expressed genes
I'm working with RNA-seq data. I have 40 tumor samples and 5 Normal samples. Differential analysis with `Deseq2` based on Fold change \> 1.2 and alpha \< 0.05 gave very low number of differentially expressed genes. Only 2 upregulated genes.

```
res <- results(dds, lfcThreshold = log2(1.2), alpha = 0.05)

```

To get more number of differential expressed genes I have few questions now

1. Instead of FDR \< 0.05 can I use FDR \< 0.1 (or) FDR \< 0.5. Will there be any low confidence with this?
2. Can I select differential expressed genes only based on FDR \< 0.05 without any fold change cutoff?
3. As I get very low number of DEGs with FDR \< 0.05 & Foldchange \> 1.2, Can I select DEGs based on Foldchange \> 1.2 and p.value \< 0.01 or 0.05 ?

PCA of the samples looks like this:

# Answer
1. You can usually get away with FDR \< 0.1, but that's as high as you can go. This all presumes you're doing follow-up experiments of some sort, of course. I guess you could increase the FDR more, but you're then really increasing the odds that your follow-up experiments will fail. Obviously increasing the FDR will decrease the confidence in the results, it's called "false discovery rate" for a reason.
2. Yes, fold-change cut-offs are mostly useful when you need to prioritize hits further. Small changes are almost impossible to validate with less-sensitive methods (qPCR) and less likely to be biologically relevant.
3. You're just going to waste your time if you ignore FDR.

Remember that all thresholds are intended to increase the odds of success in follow-up experiments. If you have infinite time and resources then you can set these to whatever you want, but you have neither so you better prioritize reasonably.

It's likely that you have a large batch effect in your data, effect sizes in cancer tend to be much larger than this.
> 3 votes

# Answer
I think since your PCA shows that your normal samples cluster right alongside the tumor ones, you have to conclude that there is much more variation between tumors than between tumors and controls. DESeq's differentially expressed gene list was generated properly.

Your results are what they are, sorry they are so useless.
> 1 votes

# Answer
In general, messing with FDR or log fold change cutoffs beyond the ones you used gets pretty sketchy. I usually only adjust these values in the more conservative direction, although you may be able to get away with FDR \< 0.1 if you have a large effect size (which you don't).

You might consider investigating why you don't seem to get any differentially expressed genes before adjusting your cutoffs: PCA on the samples based on their top 1,000 most variably expressed genes would be a good start. Do the tumor and normal transcriptomes cluster apart from each other?

With so few normal samples, it is possible that you are not fully characterizing the variability of that population and it is affecting your ability to identify meaningful differences between normal and tumor.
> 0 votes

---
Tags: r, rna-seq, differential-expression, deseq2, rna
---


# Question
Title: GMAJ error: Exception in thread "main" java.lang.ClassCastException: java.lang.StringBuffer cannot be cast to java.lang.String
I have installed GMAJ using:

```
wget http://globin.bx.psu.edu/dist/gmaj/gmaj.zip  -O output/software/gmaj.zip
unzip -u output/software/gmaj.zip -d output/software/gmaj/
rm output/software/gmaj.zip

```

I ran lastz on two chromosomes from two different assemblies of *C. elegans*:

`lastz c_elegans.PRJNA13758.WS263.genomicI.fa I.fa --format=maf > MAF.maf`

and got a MAF file.

Then, I tried to visualise the file using:

`java -jar output/software/gmaj/gmaj.jar MAF.maf`

but I got the following error:

```
Exception in thread "main" java.lang.ClassCastException: java.lang.StringBuffer cannot be cast to java.lang.String
    at edu.psu.bx.gmaj.MajGui.setDefaults(MajGui.java:272)
    at edu.psu.bx.gmaj.Maj.<init>(Maj.java:53)
    at edu.psu.bx.gmaj.MajMain.main(MajMain.java:86)

```

Can someone help me resolve this error or suggest an alternative way to visualise my alignment?

I work on Mac OS X.

Cheers.

# Answer
It's a programming error that should look like this:

```
StringBuffer sb = (...)
(...)
String s = (String)sb;

```

while it should be

```
StringBuffer sb = (...)
(...)
String s = sb.toString();

```

( see also https://www.biostars.org/p/177992/ ) . You cannot do anything about it but asking the authors of the software.
> 1 votes

---
Tags: sequence-alignment, visualization, software-installation, maf
---


# Question
Title: Where to download baseline/average gene expression level of all human coding genes?
I am looking for the most appropriate dataset for downloading baseline gene expression level across all human coding genes during development. I am aware that EMBL Expression Atlas is one of the resources that provide such information, but I am ideally looking for species averaged (and probably scaled) expression data for all human coding genes during early stages of development! Any insight is much appreciated!

# Answer
> 3 votes
There is not going to be a public database with this data. Apart from anything else, generating data for early human development is difficult and ethically tricky. Also, most people who care about development would probably think that a whole body average would be meaningless.

The best you could do would be to select tissues you were interested in from something like Expression Atlas and calculate the averages yourself.

---
Tags: public-databases, data-retrieval, gene-expression
---


title: Psipred installation

text:
# Question
Title: Psipred installation
I am trying to install psipred in a unix server.

Psipred is installed, but fails because it can not find the correct database files. I have downloaded uniprod90.fasta. Psipred also loads the old legacy ncbi BLAST-2.2.26 and BLAST+.
`runpsipred` and `runpsipredplus` have been patched to call the correct version of BLAST.

I have run a few commands to index the uniprod90 database:

```
pfilt uniref90.fasta > uniref90filt
formatdb -t uniref90filt -i uniref90filt
makeblastdb -in uniref90.fasta -dbtype prot -input_type fasta -out uniref90.fasta

```

I then ran `psipred` and got this error:

```
$ ./psipred/example/example.fasta
Running PSI-BLAST with sequence  ./psipred/example/example.fasta ...
[blastpgp] WARNING: Unable to open uniref90.pin
[blastpgp] WARNING: Unable to open uniref90.pin
FATAL: Error whilst running blastpgp - script terminated!

```

`makeblastdb` created 26 pin files, but not uniref90.pin.

# Answer
Try with:

```
makeblastdb -in uniref90.fasta -dbtype prot -input_type fasta -out uniref90

```

The `-out` is the name of the DB. So if you use `uniref90.fasta`, it creates uniref90.fasta.pin`while the script seems to be expecting`uniref90.pin\`.
> 4 votes

---
Tags: software-installation
---


title: Odds ratio calculations in GWAS paper

text:
# Question
Title: Odds ratio calculations in GWAS paper
Posting here because the reddit thread I made didn't gain any traction.

I'm having a little bit of trouble understanding how the odds ratio was calculated in this table, specifically the OR marked in red:

To my understanding, an odds ratio = $$\[(\text{# of cases with exposure}) / (\text{# of controls with exposure})\] / \[(\text{# of cases without exposure}) / (\text{# of controls without exposure})\]$$

In the red-highlighted case, my odds ratio is

$$(0.661 / 0.6580) / \[(1-0.661) / (1-0.6580)\] = 1.0134$$

which is very off from $0.98$. Also, I know I'm using frequencies rather than absolute counts, but it shouldn't matter because if I were to multiply each frequency by the corresponding population values, all the population values would cancel out.

So, what have I done wrong?

Link to the paper

# Answer
> 2 votes
Those odds ratio aren't calculated from the raw allele frequencies. They're calculated from coefficients estimated from complex linear models with a bunch of covariates that are run on subsets of the data then meta-analysed together.

It is a bit surprising that the odds ratio is less than one when the frequency is higher in cases, but the p-value indicates that the OR is not significantly different from 1 so it may be down to statistical noise. It's possible that within populations one allele is more common in cases than in controls but that cases are primarily sampled from populations where that allele is uncommon, causing it to be more common in controls overall. See here. Alternatively, 1/0.98 is pretty close to the value you calculated, so it's possible that the reference allele got switched around at some point.

---
Tags: gwas
---


# Question
Title: size of the pathways for analysis and filtration
I have recently started working on a substance's effect on a cell line in different dosages. for this, there is a tool called bmdexpress2 that I am using. Its input is the normalized counts from RNASeq for each dosage as a big matrix. When it comes to the pathway analysis step, unlike hallmarks of GSEA, this tool uses some databases for defining pathways which involve pathways of even 2-3 genes.

So what I want to discuss here is the strength of small pathways and filtering thresholds. How can we bioinformatically decide on the importance of a pathway made of 3 genes? Should we just filter them out as there are much more comprehensive pathways? Or do their GO Levels matter? Also there are many cases that out of those 3 genes, 1 of them is differentially expressed; is having 1 gene diff. expressed out of X genes in a pathway enough to say that the pathway is enriched?

Here is an example GO term for a small pathway.

And here is an example of the output I get from the tool. all genes platform column is the number of genes involved in that pathway, input genes column is how many of the genes in that pathway are significantly differentiated. I believe the rest is self explanatory.

Looking at the table, can we really say, ESCRT III Complex, for example is a good pathway to focus while having 1 gene out of 13 in a pathway but that change is significant?

Thanks

# Answer
The problem here is that any one GO term with a small number of members is unlikely to include any differentially expressed genes by chance, but if you test a large number of GO terms, some of them definitley will. The solution is to correct the p-values for multiple comparisons. If you do this, then none of those GO terms with a single DE member will be significant. The larger GO terms with multiple DE members in the chart you show probably won't survive correction either. Hopefully you have some that are more significant.
> 4 votes

---
Tags: rna-seq, pathway
---


title: Convert Cotton Probe ID to Gene Symbol

text:
# Question
Title: Convert Cotton Probe ID to Gene Symbol
I am new to bioinformatics, my background is in Electrical Engineering. I am trying to convert Affymetrix Cotton Probe IDs to gene symbols. I have a gene expression dataset and I need the expressions for only certain genes. So in the dataset, I have a probe id such as Ghi.10301.1.S1\_s\_at, how can I get the gene symbol associated with it? Thanks

EDIT: This is the specific dataset I am looking at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50770

# Answer
> 1 votes
Your dataset was produced on a microarray. The probe IDs in the FASTA file are similar to read IDs produced from other NGS experiments like RNA-Seq. According to the workflow provided by Ensembl the first step is to align your reads to a reference. You'll only keep the reads that match the reference completely. Once the probes are aligned you'll be able to quantify the number of reads mapping to the gene that the probe was designed to attract.

I did find a file providing the database entries associated with each probe ID (one of which is the ID you provided) for the Affymetrix Cotton Genome Array.

# Answer
> 1 votes
The first stop for decoding probe IDs for a microarray platform is usually the GEO Platform Description: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL8672

This description contains a table of annotation, and is provided (in this instance) by the array manufacturer. You can view, or download, the full table: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc=GPL8672&id=1538&db=GeoDb\_blob35

GEO Platform annotations do date, and so for up-to-date annotations, download from NetAffx: http://www.affymetrix.com/support/technical/byproduct.affx?product=cotton

The most recent annotations (dated 2016-04-13) are available here: http://www.affymetrix.com/Auth/analysis/downloads/na36/ivt/Cotton.na36.annot.csv.zip.

For the example in your question, looking at this probe in the CSV:

```
"Ghi.10301.1.S1_s_at","Cotton Genome Array","Gossypium hirsutum","Mar 30, 2016","Consensus sequence","GenBank","Ghi.10301.1","gb:DT050231 /DB_XREF=gi:72275634 /DB_XREF=COT_DW_H07 /TID=Ghi.10301.1 /CNT=21 /FEA=EST /TIER=ConsEnd /STK=1 /UG=Ghi.10301 /UG_TITLE=Transcribed locus, moderately similar to NP_568148.1 ubiquitin conjugating enzyme ubiquitin-like activating enzyme (Arabidopsis thaliana) /REP_ORG=G. hirsutum","DT050231","Ghi.10301","Ghi.10301","---","---","---","---","---","est","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","---","This probe set was annotated using the Accession mapped clusters based pipeline to a UniGene identifier using 21 transcripts. // false // Accession mapped clusters // E","AI730871,DN803875,DR452973,DR457278,DT050231,DT545642,DT546719,DT546982,DT546983,DT547780,DT547781,DT555433,DT559909,DT564450,DT566924,DT569293,DT569974,DT571307,DT571318,DT571472,DT573548","AI730871 // /DB_XREF=BNLGHi8130 // gb // --- // --- /// DN803875 // /DB_XREF=G.hir-8-10 DAA bolls drought stressed751 // gb // --- // --- /// DR452973 // /DB_XREF=CM008C10 // gb // --- // --- /// DR457278 // /DB_XREF=CM053F01 // gb // --- // --- /// DT050231 // /DB_XREF=COT_DW_H07 // gb // --- // --- /// DT545642 // /DB_XREF=EST1056282 // gb // --- // --- /// DT546719 // /DB_XREF=EST1057359 // gb // --- // --- /// DT546982 // /DB_XREF=EST1057622 // gb // --- // --- /// DT546983 // /DB_XREF=EST1057623 // gb // --- // --- /// DT547780 // /DB_XREF=EST1058420 // gb // --- // --- /// DT547781 // /DB_XREF=EST1058421 // gb // --- // --- /// DT555433 // /DB_XREF=EST1066073 // gb // --- // --- /// DT559909 // /DB_XREF=EST1070549 // gb // --- // --- /// DT564450 // /DB_XREF=EST1075090 // gb // --- // --- /// DT566924 // /DB_XREF=EST1077564 // gb // --- // --- /// DT569293 // /DB_XREF=EST1079933 // gb // --- // --- /// DT569974 // /DB_XREF=EST1080614 // gb // --- // --- /// DT571307 // /DB_XREF=EST1081947 // gb // --- // --- /// DT571318 // /DB_XREF=EST1081958 // gb // --- // --- /// DT571472 // /DB_XREF=EST1082112 // gb // --- // --- /// DT573548 // /DB_XREF=EST1084188 // gb // --- // ---","---"

```

---
Tags: gene-expression
---