# Question
Title: Extracting expression data from GSE dataset downloaded from GEO
I have downloaded GSE16146 dataset from GEO using GEOquery R package.

I would like to extract "Data table" from downloaded GSE16146.

```
>library("GEOquery")
>GSE16146 <- getGEO("GSE16146")
>Table(GSE16146)[1:5,]

```

This returns the following error:

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

```

Can anyone help?

# Answer
> 8 votes
According to the manual, all you need to do is:

```
library('GEOquery')
gseGSE16146 <- getGEO('GSE16146', GSEMatrix=FALSE)

```

As explanation, `getGEO()` outputs by default to `GSEMatrix=TRUE` and returns a list of `ExpressionSet` objects. You should get what you were looking for with:

```
Table(GSMList(gseGSE16146)[[1]])[1:5,]

```

The manual has also a paragraph about this:

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

I tested the code I posted on `R 3.3.2` and I get access to the Data Table just fine.

update:

Answering how to extract the actual expression data (not what asked in the original question):

```
library(GEOquery)
data = getGEO("GSE16146")
datExpr = exprs (data[[1]])

```

---
Tags: r, bioconductor, microarray, geoquery
---


title: JSmol - hide "JSmol" logo

text:
# Question
Title: JSmol - hide "JSmol" logo
I'm using JSmol for my SimRNAweb server. I would like to hide "JSmol" logo from all small visualizations that I have there. I would like to just put information about JSmol under the panel in text.

# Answer
> 5 votes
OK, I got an answer here. Use:

`frank off`

to get:

---
Tags: rna-structure, 3d-structure
---


# Question
Title: Problems to Extract uniquely mapping reads from BWA MEM alignment
I did a mapping of genomic paired-end reads to a reference assembly using bwa mem. I need to extract the reads that mapped only once to my reference. For that, I have tried to follow the method proposed in another post Obtaining uniquely mapped reads from BWA mem alignment and unfortunately in my case, after using:

```
$ samtools view -h output.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d

```

I still obtain some outputs. Here is an example:

```
NB501461:7:HH3GGAFXX:4:21612:26866:20400    99  ctg7180000256968    441 0   147M4S  =   560 230 AAGCTTTAGTCGCATGGCCATTCACAAGACCAGTGATTTTTTATAGAAGTCGTTCAAGAGATACTGATGAGGAGTACAGCCATTCTTTAAGCTATACAAAAAGCACCACCCACACCACACAATCCTCAAAAAGATCAAGGCAGTACAGATT AAAAAEEAEEAEAAEEEAAEEEAEEEEE/EEEEEEAA/EEEEEEEEEEEEEEAEE/EEEE///E6EAEEEEAEA<EE/EA/</EEAE/EEEEE//EEEEEE6E/<EAE/EEE/EE</E/E//AE<6EA//A<A///<A/6A/6//A///// NM:i:7  MD:Z:34A56T11A4A3A2A22T8    AS:i:112    XS:i:117
NB501461:7:HH3GGAFXX:4:21612:26866:20400    147 ctg7180000256968    560 0   32S111M8S   =   441 -230    CCAGCTCCACCACCCGACCCCCCCACCCACCCCAATCCTCAAACAGATCCATGCAGTACATTTTCATTTGTTCTGGCCCTGCGGCTCTTAGTTTAGGAAATATTCAATGACTTGCTACCTAGTAAGAACATTATTGTGCACAACGAATAGA ///E/////A////A/6A/EE/<//////A//AA/EEE/A/<E///A/A<//</////A/////<///EE////E/E/E/E</E//EA//E//A/</E//E//AE<//A<///EE///////E//E//EA/A//EA/AE/A/E/E//6/A/ NM:i:7  MD:Z:11A5A22A4A4T27A7A24    AS:i:76 XS:i:79
NB501461:7:HH3GGAFXX:4:21612:6399:20400 113 ctg7180000282074    357 8   73S50M28S   ctg7180000432624    28  0   CTTATCCAATTTTTAAACATAGGATTAGGGTACCTTCCAAATGTTATTGAAACTTAATCTGAAGATGACCTCAAATGGAAGATGACCTCCAATGGAAGACAACCTTCCATGGAAGACGACCTCAAAAATGAAGAAGAAGAAAACATCTAAG /A//<//<<//E/66A//<A////A/<</E/EA/EE//<E<E//<AE//A/EEAEEEE/E//A/A/</EEE//E<EAEEEEE/EEEEEE<EEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEAAAAA NM:i:2  MD:Z:32C1A15    AS:i:40 XS:i:35
NB501461:7:HH3GGAFXX:4:21612:6399:20400 177 ctg7180000432624    28  26  76S75M  ctg7180000282074    357 0   CTCCTTACTGTTGTTTCTTCTTTTTTCTTCTCTTTCATTTTTGAGGTTGTCTTCCATGTTTTGTTGTCTTCCATTTGAGGTCGTCTTCATATTAAGTCTCTATAACATTTGGCTTTGTCCCCTAAGCCTAAGTTTCAAATTTGGAAAAGAT /////</////6/////A////EAEE<//////<A////E/E<//A//<///E//<E/////E/A/E//E/EAEEAE/EA///E//E/AE/<EAEE/6E6EEA/EE/AEEAE6EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA NM:i:4  MD:Z:5T15T17G29G5   AS:i:55 XS:i:44

```

As I do not see any XA or SA flags there anymore, I wonder what those can be?

I have tried to continue the steps I want to pursue until the "ValidateSamFile" from Picard tools, and there I do not know if it is linked or not, but I get the message: "Error: mate not found" for a few thousands of reads. Any idea what could be wrong, and how to fix it?

I hope I was clear enough, otherwise please feel free to ask me for specific details. Any help is welcome, Thanks!

# Answer
> 1 votes
The issue is that you are filtering reads without regard to the pairing so you might filter out one of the mates if one aligns uniquely and the other does not. You will need to filter out the un-mated reads from the file to have a valid BAM.

> As I do not see any XA or SA flags there anymore, I wonder what those can be?

I am not sure what this means. The grep should remove those flags and only leave records without them.

As a side note if you don't want split alignments I would recommend using a different aligner such as BWA aln or BowTie2

# Answer
> 1 votes
> extract the reads that mapped only once to my reference

Every read can be mapped to multiple places if you allow enough mismatches. You heard "unique mapping" from time to time. Strictly speaking, however, this concept is either ill defined or practically useless given 150bp reads. For a proper aligner, the best way to measure the reliability of an alignment is simply mapping quality. All the other approaches, including the answer in the other thread, have problems.

By the way, many tools allow you to set a threshold on mapping quality. You can apply the threshold on the fly. It is a common and good practice to avoid writing filtered alignments to another BAM.

---
Tags: read-mapping, bwa
---


# Question
Title: How can I output an identity matrix in progressiveMauve?
I'm just getting started with the Mauve aligner and I'm finding the documentation a bit lacking. I'm using the `progressiveMauve` tool from the command line and would like to output an identity matrix file. The set of output files it generates by default for `progressiveMauve --output=outfile` are `outfile`, `outfile.backbone` and `output.bbcols`.

The documentation for the tool lists an `--input-id-matrix` option but this doesn't seem to have any effect. The documentation page for output files does describe an identity matrix file but doesn't say anything about how to generate it.

# Answer
> 1 votes
mauve documentation have a similar command but the usage is `--id-matrix=<file>`. So very similar to the progressivemauve command (just doesn't include --input). Try it and see if it works.

I think it's a bug for sure as several other people on biostar have complained.

# Answer
> 0 votes
It indeed seems that the `--input-id-matrix` option doesn't work as expected, as Amar has pointed out.

I tried the following command:

```
./progressiveMauve genome1.fasta genome2.fasta genome3.fasta --output=test > test.log

```

And saw the following bits in `test.log` (truncated):

```
[...]
Genome conservation distance matrix:
0       0.0767545       0.187454
0.0767545       0       0.193385
0.187454        0.193385        0

[...]
genome content distance matrix:
0       0.0767545       0.187454
0.0767545       0       0.193385
0.187454        0.193385        0

bp distance matrix:
0.9     0.0125  0.075
0.0125  0.9     0.0375
0.075   0.0375  0.9
[...]

```

Could any of those three tables contain the information you are looking for?

p.s. as an aside, if you are looking into a way to obtain a whole genome similarity matrix you might want to look into mash or sourmash.

---
Tags: sequence-alignment
---


# Question
Title: Counting a specific consecutive character with its occurrence position and length
I have a sequence file and want to count consecutive character "N", with the tandem's position of occurrence and its length. Say I have a file named `mySequence.fasta` like this:

```
>sequence-1
ATCGCTAGCATNNNNNNNNNNNNNNCTAGCATCATGCNNNNNNATACGCATCACANNNNNNNNNCgcatATCAC

```

and anticipated output should be like this:

```
Position 12 N 15
Position 38 N 6
Position 56 N 9

```

Kindly help me to solve this by `awk` or `sed` providing my file name `mySequence.fasta`

# Answer
> 2 votes
I couldn't come up with a one liner but have a two liner instead...

To get the start positions of the N repeat:

```
grep N mySeq.fasta | grep -b -o N | awk -F: '{print $1}' |awk '$1!=p+1{print $1}{p=$1}'
12
38
56

```

To get the length of the repeats:

```
grep N mySeq.fasta | awk -F '[^N]+' '{for (i=1; i<=NF; i++) if ($i != "") print length($i)}'
14
6
9

```

With thanks to stack overflow.

---
Tags: fasta, text-processing
---


# Question
Title: Read and write FASTA files with more information than id and sequence
I am trying to read a fasta file, manipulate is in Python (using BioPython) and then write it back.

The format of my sequences is like:

```
>k119_5 flag=0 multi=141.0706 len=473
AGGTTAGTCAGCACCGTTTCCGTGGTGCTGCCTTTCGCTTCAAAACCGACGGCGTCTATTACTGCATCCACGCCCCGGTGACCTGCCGTTTGTTCAATAATTGACTGTGCCGGATCGCTGTCTTCATCAAAATTAATCGGGATCGCGCCGTAGCGGTCGGCGGCGAAATGCAAGCGGTAGGGATGATGATCAACAACAAAAATCTGTTCCGCACCGAGCAACCGTGCACAGGCGATTGTCAACAATCCCACAGGACCAGCACCATAGACTGC
AACGCTTGAACCTTGTTGGATCTGCGCATTTTTTGCTGCCTGCCATGCCGTTGGCAGAATATCAGAAAGGAAAAGCGCTTTATCATCTGAAAGCAAAGGCGGTACTTTAAACGGCCCCACATTCCCTTTAGGGACGCGGACATATTCCGCCTGCCCACCAGGAACGCCGCCATACAGGTGACTATAACCAAACAATGCCGC

```

which is the result of the assmebly using Megahit. While trying to read, it only recognizes `k119_5` as the id, while I want the whole first line as the id of the sequence. Is vital information I want to keep!

Any ideas how to work with this?

# Answer
> 8 votes
The `description` field in the SeqRecord object has the information you are looking for:

```
>> from Bio import SeqIO
>> s = SeqIO.read('genome.fasta', 'fasta') # single sequence fasta file
>> s.id
k119_5
>> s.description
k119_5 flag=0 multi=141.0706 len=473

```

Edit: as an aside, if you write a `SeqRecord` object using `SeqIO`'s `write` method, the FASTA header should contain the whole original header.

---
Tags: fasta, python
---


# Question
Title: What will be an appropriate mathematical distribution for SNP data?
I found that several papers describe SNPs as a binomial distribution with the probability of "success" equals to minor allele frequency.

However, in my experiments, when I generate SNP array following this distribution, the simulation results behave very differently from the SNP array generated following the random mating procedures. I feel like there is some missing information about the description of a binomial distribution.

I wonder what will be a better mathematical description of the SNP array data? Besides a direct answer, any thoughts, or any suggestion of papers are also highly welcome.

**EDIT:**

Thanks for these suggestions, after reading relevant materials, I think what I observe is about the Hardy-Weinberg equilibrium affected by evolutionary influences like mutation and migration. I will continue to look and list findings here if I find anything interesting. Please feel free to suggest some mathematical descriptions.

Also thanks for the suggestion of tools. They look very helpful, but I really hope to understand some mathematics behind it.

# Answer
> 6 votes
Simulating genotypes with realistic correlation structures is indeed not so simple, and there's quite a few papers dedicated entirely to that (e.g. https://bmcgenet.biomedcentral.com/articles/10.1186/s12863-015-0173-4). Also, DEPICT (https://data.broadinstitute.org/mpg/depict/index.html) comes with a number of simulated GWASs to generate the nulls, so that's a simple way to obtain some ready-made fake data.

On the other hand, indepedent SNPs *should* behave as draws from $Binomial(2, MAF)$. Note that the standard QC procedure of filtering out genotypes that deviate from Hardy-Weinberg equilibrium is just a goodness-of-fit test against $Binomial(2, \hat{MAF})$. If you employ this filter, some LD-pruning, and still get different results, consider posting more details about the discrepancies you see - that would be indeed unexpected.

---
Tags: statistics, snp
---


# Question
Title: Database of copy number alterations in different cancer types
Is there a database of somatic copy number alterations (CNAs) in different cancer types? That is, annotations of which CNAs are more commonly seen in different cancer types.

# Answer
I don't know if there's a database that does *exactly* what you want, but there are some places that might help you figure this out, especially if you already have a list of CNAs/genes/regions in mind.

ICGC is the first that comes to mind, as it has samples from many types of cancer and copy number data for most of them. Depending on what you want though, you might have to apply for access to download the data you want. It's rather annoying, and their data download tools are atrocious to use, but they have a **lot** of data.

DGV might be another option, though their data browser is rather dated feeling. They publish a lot of gold-standard germline CNVs and such.

COSMIC is a large database for cancer-associated mutations, including CNVs, so it might be worth a look for you.
> 7 votes

---
Tags: public-databases, cancer, cnv
---


title: Filtering step for read counts data

text:
# Question
Title: Filtering step for read counts data
I have around 1200 samples as columns and 60,000 genes with Htseq-Counts data. Before normalization with voom function I want to do filtering step.

I want to remove genes whose expression is == 0 in at least 10 samples.

Can I do this with read counts itself or should I transform counts to cpm?

If filtering with counts data is fine May I know how to do this.

# Answer
If you have your counts in a data.frame called `counts`, something like this might work:

```
filtered.counts <- counts[rowSums(counts==0)<10, ]

```

For example lets assume the following data frame.

```
> A <- c(0,0,0,0,0)
> B <- c(0,1,0,0,0)
> C <- c(0,2,0,2,0)
> D <- c(0,5,1,1,2)
>
> counts <- data.frame(A=A, B=B, C=C, D=D)
> counts
  A B C D
1 0 0 0 0
2 0 1 2 5
3 0 0 0 1
4 0 0 2 1
5 0 0 0 2

```

And now you would want to remove the genes with at least 3 zeroes or more:

```
> filtered.counts <- counts[rowSums(counts==0)<3, ]
>
> filtered.counts
  A B C D
2 0 1 2 5
4 0 0 2 1

```

Every gene with 3 or more zeroes is now removed.
> 8 votes

# Answer
While this answers explains how to do it I want to address when and why and which thresholds to do it.

Filtering the genes with low counts is usually done because the counts are not reliable it would be noise, specially when there are low number of samples these genes disturb power of the analysis. However with 1200 samples having 10 or more samples without any count is valuable information because you have enough samples to estimate if this is due to a latent or known variable or if it is due to noise.

The CPM are used because the size of the library per sample affects the power of the experiment. The higher the library size the easier is to find if a gene is expressed or not. CPM are calculated to take into account that effect. Filtering by counts or CPM will result in different genes filtered if you use other thresholds than 0 (for instance I have filtered those genes that the average CPM where below 1, which would be different to filtering by those that the average counts is below 1.)
> 3 votes

---
Tags: r, rna-seq, reads, filtering
---


title: checks for spike-in sequence controls

text:
# Question
Title: checks for spike-in sequence controls
I would like to know what do people verify when designing/using spike-in controls, to be used in sequencing experiments (mainly Illumina). So far I came up with this list:

* Does it align only to a given genome synthetic reference?
* Does it contain G-quadruplexes? E.g. QGRS Mapper
* Does it contain known Illumina error motifs? (Do this matter anymore for recent platforms/chemistries?) See this article
* Does it contain other motifs known to be polymerase inhibitors? Ref

Anything else?

# Answer
> 4 votes
**Disclaimer**: I'm a developer for http://www.sequin.xyz.

Sequin is a new set of spiked-in controls for next-generation sequencing, and that includes Illumina. We design controls for RNA-Seq, genome sequencing, metagenomics etc.

Please study our papers if you want more details.

* **Reference Standards For Next-Generation Sequencing** by Hardwick should be a good start.

Here's my list:

1. Measured counts against input concentration. Pearson's correlation, spearman correlation and regression slope.
2. PCA analysis
3. Whether the synthetic aligned to the synthetic reference
4. Spiked-in amount (e.g. dilution)
5. Diagnostic statistics
6. ROC Curve
7. Detection limit

(1) is important because the controls should give you correlation = 1, and slope = 1 in a perfect flawless experiment. This is a sample output (you can read more about it in our papers):

(2) This is our sample PCA (our paper has the details)

(3) is also important but you can only do that in simulations.

(5)

(6)

(7)

---
Tags: ngs, illumina, errors, spike-in
---


# Question
Title: How to build a BAM header file with htslib in C++?
I'd like to use C++ to generate a new BAM file programmatically. This is an example how to use htslib to generate a new BCF file on the fly.

> https://github.com/samtools/htslib/blob/develop/test/test-bcf-translate.c

The two key functions are:

* `bcf_hdr_append`
* `bcf_hdr_add_sample`

Now when I navigate to the header file for SAM/BAM functions:

> https://github.com/samtools/htslib/blob/develop/htslib/sam.h

The only "append" function is `bam_aux_append`. But that's for optional fields.

My questions:

* Given a `bam_hdr_t` header (potentially empty), how to add a new chromosome reference? (For example, add a new chrAB to the BAM file)
* Given a non-empty `bam_hdr_t *` header, how to remove a chromosome? (For example, remove chr4 from the BAM header)

My sample code:

```
// Empty BAM header, nothing inside
bam_hdr_t * hdr = bam_hdr_init();

// I'd like to add a new synthetic `chrAB` into the header, how?

```

# Answer
You're just looking in the wrong header file. The BCF stuff will be in `vcf.h` (you'll see the two functions you mentioned there).

With BAM files one doesn't typically need to add/remove headers, which is why there are not convenience functions for them. To add a header:

1. Increment `hdr->ntargets`
2. `realloc(hdr->target_len, ...)` and `hdr->target_name`.
3. You likely want to append the name/length to the end of `hdr->target_name` and `hdr->target_len`, respectively.

I don't think you need to update the dictionary or plain text for most uses (or at least I've never needed to).

Removing a chromosome would follow the reverse instructions, (you'd need to shift everything in the name and length arrays if you remove entries in the middle.

# WARNING, doing this is prone to error!

I would strongly cautious you to be very careful when doing this. Please remember that the chromosome that entries are mapped to isn't actually stored in alignment entries themselves. Rather, they simply have an index into `hdr->target_name`, which saves space. Consequently, if you start adding/removing headers to/from a BAM file with alignments then you're going to have to modify **all** of the alignments to keep the correct alignment:chromosome associations. Otherwise you'll end up with completely useless results (this is the same caveat as with `samtools reheader`).
> 4 votes

---
Tags: samtools, c++
---


# Question
Title: Using the t-SNE algorithm on microarray data + an error bonus
I'm trying to use the t-SNE algorithm on some microarrays data. More specifically my data frame has 18600 columns with genes (features) and 72 rows with conditions with replicates ( 10xWt , 10xTg , etc ). The expression values are in log2 scale.

Here is the code that I'm trying to run.

```
# t-SNE implementation
library(Rtsne)

set.seed(1)

for(i in 1:15){

  tsne = Rtsne(data.T[,-18601], dims = 2, perplexity=i, verbose=TRUE, max_iter = 1000, pca=T)

  colors = rainbow(length(unique(data.T$classes)))
  names(colors) = unique(data.T$classes)
  plot(tsne$Y, t='n', main="tsne")
  text(tsne$Y, labels=data.T$classes, col=colors[data.T$classes])

  readline(prompt="Press [enter] to continue")
}

```

Please note that I'm not counting the column 18601 because this colums contains the labels/classes for each condition.

The think here is that when I execute this script, R returns me this error:

> Error: protect(): protection stack overflow

Should I change the `--max-pp-size` or it's a bug in Rtsne package?

Also I was wondering if it is more meaningful to run the tSNE algorithm using not the log2 values of the expression level but the log fold change values in respect to the Wt (wild type) condition. I'm asking because I couldn't find a such other implementation of the tSNE on microarray data.

For the configuration of the Rtsne function I read this article

Any other suggestion on the implementation is welcomed.

# Answer
Converting your `data.frame` to a `matrix` (and then removing the `data.frame`) will often free up enough memory that you won't run into this. Note that a `matrix` is more memory efficient than a `data.frame` and you're requiring `Rtsne()` to hold both in memory at the same time (many math-centric functions will end up converting things to a matrix at some point for efficiency).

For what it's worth, it's never been entirely clear to me what the interaction is between a `data.frame` and the pointer protection stack, but it's often the case that this solves this sort of error.
> 13 votes

---
Tags: r, visualization, microarray, clustering
---


# Question
Title: True bam file encoding and viewing as text
We all have been using samtools to convert between sam and bam files. SAM is tab-delimited file. My question resides on the ability to visualize bam files as text. Is there any way to visualize any bam file as text in linux? Why not work only with sam files and we need to convert them to bam?

# Answer
> 8 votes
> Why not work only with sam files and we need to convert them to bam?

SAM is for human only, a binary format like BAM is smaller and parsed much faster by a program (and we want speed and optimize spaces for those huge files, don't we ?).

Furthermore, bam files are compressed using a method that allow them to be accessed using **random-access**, that is how you can query a sub-part of bam file for a given genomic location.

you can always go back to SAM using samtools.

# Answer
> 3 votes
bam file is a compressed sam file, comparable with zip or tar. They are binary, which means we (humans) cannot read them. So you need to unpack them (convert back to sam) to read it as a text file.

bam files are smaller than sam files, which is the reason most people like to use bam files instead of sam files.

# Answer
> 3 votes
The way to visualise BAM files as text in Linux is with tools like samtools.

If you want to see the raw binary information as slightly more human-readable, you can run `hexdump -C <file.bam>`, but that type of text output is completely useless for further processing.

If you want more specifics about the BAM file format, you can have a look at the SAM file format specification, which also includes details on the related BAM and CRAM formats.

# Answer
> 3 votes
Assuming you want to visualize the pileup of reads across the reference, you may want to look into ASCIIgenome, a command line visualization tool for sam/bam/cram files, similar in spirit as IGV, but working in the command line.

```
ASCIIGenome my_sorted.bam

```

---
Tags: sam
---


# Question
Title: How to deal with correlated SNPs in GWAS?
I have population SNP data in VCF format, and I found that some SNPs have a great similarity across samples(\> 99%). For example:

```
CHROM_POS   s0  s1  s2  s3  s4  s5  s6  s7  s8  s9
chr1_1  A   G   G   A   G   G   A   A   G   A
chr1_2  A   G   G   A   G   G   A   A   G   T
chr1_3  C   C   C   A   C   A   C   C   C   C

```

The similarity between chr1\_1 and chr1\_2 is 0.9, because they have a different SNP in sample **s9**.

Is there a good way to remove these similar SNPs before I put the data into another pipeline?

Added: These SNPs will be used for GWAS analysis.

# Answer
What you are attempting to do is known as **LD-pruning**.
As @Emily\_Ensembl said, it is not customary to do this for standard association tests: it is possible that one of the SNPs you remove is causal, or a better proxy for the causal locus, and would give (slightly) better association signal than the other. Even for SNPs in perfect LD, pruning is unwise because it would complicate interpretation: if SNP `chr1_2` had some obvious function, but you removed it and are left only with its proxy `chr1_1` (which could be far away), you will have a hard time connecting the dots. Typical GWAS tests are handled pretty quickly on modern computers, and are linear in SNPs, so there is no computational reason to prune either.

On the other hand, many fancier procedures would require LD-pruned data, because they either
a) are very computationally expensive, or
b) assume independent SNPs.

I suggest reading PLINK's LD pruning documentation, as it pretty much covers the standard ways to do this.
> 6 votes

---
Tags: snp, gwas
---


# Question
Title: Adjusting phenotypes by regressing out covariates
I'm trying to use the bfGWAS tool, which analyses GWAS data and integrates functional annotations to identify casual SNPs (paper and github). In the user manual, it states:

> We recommend first regressing covariants out from the original phenotypes, and then provide bfGWAS the corrected phenotypes (i.e., residuals from the regression model with covariates).

It is unclear to me what this means and how to go about determining the corrected phenotypes. Typically, in an association test with covariates, a regression would be performed using a model like

```
pheno ~ snp + covariate

```

I could then get the residuals from this model. However, I'm not certain how to go about it for all SNPs, and then adjust the phenotypes. I have case-control data, so the phenotypes are binary.

# Answer
Instead of the "observed phenotype", it's asking for a "corrected phenotype", which is the residuals from `pheno ~ covariates`. In R, one would get that as follows:

```
fit = lm(phenotype ~ covariates, data)
corrected_phenotype = resid(fit)

```

Of course `covariates` is just a place hold for something like `age + gender` or something along those lines.
> 2 votes

---
Tags: statistics, gwas
---


# Question
Title: Interpreting Intergrative Genomic Viewer (IGV)
I was following a tutorial on "Tuxedo Genome Guided Transcriptome Assembly Workshop" and was wondering how to interpret the following:

From what I understand from 'Color Legends', the color blue represents something that is below normal and red represents something that is above normal.

Also, my understanding is that merged.gtf is the reconstructed transcript and genes.bed is the reference. They look to be similar for the most part, but seems to have breaks between. What's the appropriate way of interpreting this result?

# Answer
You've apparently colored your alignments by read strand. In this case, red indicates "+ (watson) strand" and blue indicates "- (strand) strand". This strand association is determined by the orientation of read #1 in a pair (or just "the read" if you have single-end data). This isn't actually documented in the user manual, but you can find one of the Broad institute software engineers replying here.

Your next question is likely why the color of the reads and the orientation of the genes are opposite of each other. This is because in modern directional protocols, read #2 in a pair matches the orientation of the gene from which it arose.
> 6 votes

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


title: Meaning of BWA-MEM MAPQ scores

text:
# Question
Title: Meaning of BWA-MEM MAPQ scores
Does anyone know what the `MAPQ` values produced by `BWA-MEM` mean?

I'm looking for something similar to what Keith Bradnam discovered for Tophat v 1.4.1, where he realized that:

> 0 = maps to 5 or more locations
>
> 1 = maps to 3-4 locations
>
> 3 = maps to 2 locations
>
> 255 = unique mapping

I'm familiar with the notion that MAPQ should be theoretically be related to the probability of an "incorrect" alignment (10^(-`MAPQ`/10)), although this is vaguely-specified enough that aligners in practice tend to actually just use something like the above.

Note that this is different than the `BWA` `MAPQ` scoring interpretation, because `BWA-MEM` gives `MAPQ` scores in the range $\[0,60\]$, rather than $\[0,37\]$ as has been established for `BWA`.

# Answer
> 16 votes
First of all, if you want to understand mapping quality (mapQ), ignore RNA-seq mappers. They often produce misleading mapQ because mapQ is not important to RNA-seq anyway.

Strictly speaking, you have two questions, one in the title: the meaning of mapQ; and the other in a comment: how mapQ is computed. On the meaning, mapQ is nearly the same as baseQ – the phred scaled probability of the alignment/base being wrong. It often amuses me that we question mapQ but take baseQ for granted. BaseQ is also scaled and discretized differently; even fewer people know how Illumina/pacbio/nanopore/historical sequencers estimates baseQ.

On the second question, Section 2 of the MAQ supplementary explains the theoretical aspects of mapQ, which is still correct today. Briefly, mapping quality consists of three components: 1) the probability of contamination; 2) the effect of mapping heuristics and 3) the error due to the repetitiveness of the reference. Only 3) can be modeled theoretically.

In case of bwa-mem, if we assume the matching score is 1, type-3 error is estimated with:

$$10/\log10\cdot\[\log4\cdot(S\_1-S\_2)-\log n\_{\rm sub}\]$$

where $S\_1$ is the best alignment score, $S\_2$ is the second best and $n\_{\rm sub}$ is the number of suboptimal alignments. Factor $\log 4$ comes from the scoring matrix. Factor $10/\log10$ is the Phred scale. This equation assumes gap-free alignment and is very close to Section 2.5.2 in the MAQ supplementary. It is ok-ish for short reads, but often overestimates for long reads. I am not aware of a practical approach in general cases. In addition to this method, you can estimate mapQ by read simulation: just try to find a function such that it fits empirical mapQ. Some have tried machine learning, too.

# Answer
> 7 votes
Those numbers are not arbitrarily picked (well... maybe 255/60/40 is arbitrarily picked).

To convert from log10 Q values like these (also used for error rates in FASTQ files) to probabilities, divide the number by 10, negate it, then raise 10 to the power of the result.

Another way of looking at it is to consider the decade to mean the number of 9s in the \[mapping\] accuracy, e.g.:

Q10: 90%
Q20: 99%
Q30: 99.9%
Q40: 99.99%

With the formula adjusted to make all the numbers in between a smooth transition.

With this transformation, a MAPQ score of 3 corresponds to a probability of an incorrect mapping of almost exactly 50%, which is as confident as you can be about a specific mapping if you know that it maps equally well to two locations:

```
> 10^-0.3
[1] 0.5011872

```

A MAPQ score of 2 is an incorrect mapping probability of 63% (i.e. close to 3 different locations), and a MAPQ of 1 is an incorrect mapping probability of 79% (close to 5 different locations). These fit with the numbers that you have described.

However, as you point out, BWA doesn't do this; it seems to ignore the number of mappable locations when considering the MAPQ score. I explore the reason why this is a problem in my question about multi-mapped reads.

I notice that there's a bit of an attempt to answer this in a similar previous question:

> $C * (s\_1 - s\_2) / s\_1,$
>
> where $s\_1$ and $s\_2$ denoted the alignment scores of two best alignments and `C` was some constant.

See the answer by user172818 for more insight into the details of BWA-MEM MAPQ scores.

---
Tags: sequence-alignment, read-mapping, bwa
---


# Question
Title: Missing data mappings in mygene.info while trying to convert Genes Ensembl Ids to Entrez Ids
I need to convert a lot of Ensembl Ids to the relative counterpart in Entrez (e.g., `ENSG00000157764` \> `673`).

I found mygene.info and it seems what I needed. Let's see the query about `ENSG00000157764` in action. We can easily find the key-value `"entrezgene": 673`

Unfortunately, the mapping seems not complete in some cases or, maybe, I'm missing something...

**Example:**

* Gene: `NDUFV2P1` *(NADH:ubiquinone oxidoreductase core subunit V2 pseudogene 1)*
* Expected Mapping `ENSG00000267809` \> `4730`

If I query *mygene* for `ENSG00000267809`, I get this result, but this time no data about the entrez id.

What it sounds pretty odd is that if I query another *mygene*'s endpoint for `NDUFV2P1`, 2 identical records are returned, respectively with Ensembl and Entrez Ids.

Is there a way to leverage such a service to get clean data? Otherwise, is there another robust conversion mapping strategy out there?

# Answer
> 8 votes
As others have mentioned in their answer, bioMart is usually the best place to go for this infomation because it draws its data directly from the Ensembl database.

However, you will find you do not get a full translation of Ensembl IDs to Entrez IDs even there. The reason for this is there simply does not exist a simple 1-to-1 mapping of Entrez IDs to Ensembl IDs: they are not just different names for the same thing, but are actually different things. At some point someone (or more likely some algo) has decided that the structure Ensembl ID A is referring to is similar enough to the structure that Enztrez ID B is referring to to link the two together, but such a linking is not always possible.

In the example you gave, the Ensembl database does not list ENSG00000267809 as being equivalent to 4730, neither does the Entrez database list the opposite, even though both genes have the same name. Because of this, it is very unlikely that any service you use to convert the IDs will pick up this match.

In this particular case, the reason for this is probably that the gene is a processed pseduo-gene and is not annotated in Refseq/NCBI as producing any transcripts.

---
Tags: gene, public-databases, conversion
---


# Question
Title: genes with highest RPKM/FPKM for a human RNA-seq experiment?
Which genes/transcripts show up as having the highest RPKM/FPKM on a human RNA-seq experiment? Is it usually the same genes/transcripts?

EDIT: mRNAs (e.g. poly-A RNA-seq preps) and not including small RNAs.

# Answer
As asked as it is, the answer is probably no.

Indeed, the highest expressed RPKM/FPKMs will be different from one condition (or tissue) to another one for example.

Then, you may also have technical artefacts due to the wet-lab part or to the normalisation. For example, mitochondrial genes are often reported in the top expressed genes.

Now, if you want to compare your results to an outside source, I recommend checking the EBI gene expression Atlas* and more specifically the Baseline expression experiments.

Pro-tips:

* You can filter the genes based on their nature (eg, mRNA) and/or their expression.
* You can also visualise the overall expression of a specific gene across several experiments.

\*Disclaimer: I am related to this project myself.
> 6 votes

---
Tags: rna-seq, ngs, illumina, human
---


# Question
Title: How to represent a deletion at position 1 in a VCF file?
I am writing a small script to write a VCF file from short multiple sequence alignments, and stumbled into an interesting corner case.

This is the way a deletion should be represented in VCF format (I assume version 4+):

```
NC_000016.9  2138199 rs1799758   GTGAG   G  .   .  RS=1799758;...

```

But what happens if the deletion is at the first position? Immagine the following two aligned sequences:

```
>reference
AGTC
>sample
-GTC

```

How should it be encoded in VCF format? According to this thread on SEQAnswers, it might be necessary to include the first position after the end of the deletion, like this:

```
chromosome  1  AG G  .  .  [...]

```

Is this a VCF 4+ correct way to encode such variants?

# Answer
> 7 votes
From the latest VCF specs (page 8):

4. REF - reference base(s): Each base must be one of A,C,G,T,N (case insensitive). Multiple bases are permitted. The value in the POS field refers to the position of the first base in the String. For simple insertions and deletions in which either the REF or one of the ALT alleles would otherwise be null/empty, the REF and ALT Strings must include the base before the event (which must be reflected in the POS field), unless the event occurs at position 1 on the contig in which case it must include the base after the event;

So I think you're right

---
Tags: vcf, variation
---


# Question
Title: Transform traditional blast output to `--outfmt 6`
I have run a blastx of metagenomic databases (raw illumina reads) using the nr database. Unfortunately, I forgot to add the `--outfmt 6` argument to the code and got the traditional output.

Could I parse this output into `--outfmt 6`?

# Answer
> 3 votes
Biopython has a parser for Blast's default text format (which I assume is the one you have). Just use the `blast-text` string to specify this format in `SearchIO`'s `parse` method.

```
from Bio import SearchIO
results = SearchIO.parse('myfile.txt', 'blast-text')
for entry in results:
    # extract the information you want
    query = entry.id
    for hit in entry.hits:
        subject = hit.id
        for hsp in hit.hsps:
            identity = hsp.ident_num *100 / entry.seq_len
            aln_len = hsp.aln.get_alignment_length()
            mismatches = aln_len - hsp.ident_num # not sure about this
            gaps = hsp.gap_num
            q_start = hsp.query_start + 1
            q_end = hsp.query_end
            s_start = hsp.hit_start + 1
            s_end = hsp.hit_end
            e_value = hsp.evalue
            bit_score = int(hsp.bitscore)
            # print in outfmt 6 format
            print('\t'.join([str(x) for x in [query,
                                      subject,
                                      identity,
                                      aln_len,
                                      mismatches,
                                      gaps,
                                      q_start,
                                      q_end,
                                      s_start,
                                      s_end,
                                      e_value,
                                      bit_score]]))

```

Have a look at Biopython's docs for more information on how to parse Blast results.

# Answer
> 3 votes
The `blast_formatter` command allows you to produce output in any of the formats supported by BLAST. However, it requires that you run `blastx` with `--outfmt 11` (ASN.1 format). I'm not sure there are any tools that will do the conversion from the full default BLAST report. This could potentially be scripted, but would be complex enough that you're probably better off re-running BLAST.

---
Tags: python, blast, metagenome, parsing
---


# Question
Title: Can we describe synthetic biological sequence data as belonging to a "species" or "taxon"?
For the purpose of testing phylogenetics software, we're using synthetic sequence data (i.e., it's computer generated, and there's no actual biological organism with that e.g. DNA sequence). For real-world biological sequence data, we would say the sequence belongs to a "species" or "taxon" (or "taxonomic group"), but this feels inaccurate since there is no actual species or taxon.

**Question**: Can we describe synthetic biological sequence data as belonging to a "species" or "taxon"?

I'm wondering if this is standard practice, despite the synthetic nature of the data.

# Answer
Phylogenies can be made with/of non-species groups. I think "taxon" is perfectly fine. You might also consider the use of the "leaf / leaves" or "terminal" if you want to emphasise the tip nature of the organisms. Even "sequence" would be fine.

Ninja edit to my own suggestion: since you're making a phylogeny from sequences of unknown context or relationship (i.e. you can't say that they're sequences from a "species" or "OTU" or from anywhere), agnostic terms like "taxon" or "sequence" are more correct.
> 6 votes

# Answer
When looking at genomic data for bacteria, taxa above the species level are operational taxonomic units (OTUs). OTUs are designated based on percent identity of the sequence (specifically the 16S rRNA gene sequences in most cases). The classification is usually based on rule of thumb, and some rough guidelines are:

* 97% for Genus
* 94% for Family
* 88% for Order

Using OTUs is common practice. I don't see much of an ideological difference using OTUs with synthetic data.

*Konstantinidis, K. T., & Tiedje, J. M. (2005). Towards a Genome-Based Taxonomy for Prokaryotes. Journal of Bacteriology, 187(18), 6258–6264. http://doi.org/10.1128/JB.187.18.6258-6264.2005*
> 2 votes

---
Tags: terminology, phylogenetics
---


# Question
Title: Where to download JASPAR TFBS motif bed file?
I am interested in determining if any transcription factor binding site motifs are enriched in some BED files from a DNA methylation experiment.

I am looking for a database that has BED Files containing regions enriched for specific transcription factor binding motifs in homo sapiens.

I know that JASPAR contains TFBS motifs, but I can't find a download for regions of hg19 that are enriched with those motifs.

EDIT: I suppose I could determine TF motif enriched regions using HOMER if there are no curated databases with this information, but I would rather not. Your help is appreciated!

# Answer
You could make your own such file with FIMO, the JASPAR MEME files available via the MEME download site, and your genome build of interest (per-chromosome FASTA files from UCSC goldenpath, say), e.g., for assembly `hg19` and UCSC chromosome naming scheme:

```
$ for chr in `seq 1 22` X Y; do echo ${chr}; wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg19/chromosomes/chr${chr}.fa.gz | gunzip -c - > hg19.chr${chr}.fa; done

```

For each chromosome (omitting the `for` loop code for brevity), you can use the `fasta-get-markov` in the MEME suite to generate a 5th-order Markov model of kmer frequencies (why 5th order?):

```
$ fasta-get-markov -m 5 hg19.chr${chr}.fa > hg19.chr${chr}.5th_order_background_model.txt

```

Here is one way to download a non-redundant, vertebrate JASPAR database (MEME-formatted):

```
$ wget http://jaspar.genereg.net/download/CORE/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.txt

```

To run a FIMO query with the 5th-order background model and JASPAR database, using a p-value threshold of 1e-4:

```
$ fimo --verbosity 1 --bgfile hg19.chr${chr}.5th_order_background_model.txt --thresh 1e-4 --text JASPAR2018_CORE_vertebrates_non-redundant_pfms_meme.txt hg19.chr${chr}.fa \
  | tail -n+2 \
  | awk -vOFS="\t" -vFS="\t" '{ print $3, $4, $5, $1":"$2">"$10, $8, $6 }' \
  | sort-bed - \
  | starch --omit-signature - \
  > hg19.chr${chr}.bed.starch

```

The output of a FIMO search is a text file which can be modified into a sorted BED file via sort-bed, an example of which is shown above.

To get a faster answer, each per-chromosome query can be put onto its own node on a computational cluster.

Using a p-value threshold of 1e-4 will result in very large files (when uncompressed). Using tools to compress BED files is useful here, as is using a more stringent threshold (1e-5, etc.).

Use of BEDOPS `starch` to compress BED files (as shown above, for instance) will give better compression performance than `bgzip` with per-chromosome random access and integration with BEDOPS `bedops` and `bedmap` CLI tools. The `bgzip` option offers integration with `tabix` and offers finer-grained random access, but it will use more disk space and doesn't integrate with BEDOPS tools.
> 3 votes

---
Tags: database, sequence-annotation, chip-seq, motifs, ucsc
---


# Question
Title: Using BLAT command line tool to blat split sequences
I have the nucleotide sequence:

*AATTGAGGCACATTTTTTTTTAGACAGTCTTGCTCTGTTGCCCAGGCTGGAGTGCAGTGGTGTGATCATAGCTCACTGCAGCCTCGACCTCCTGGGCTCAACAAAGCACACAGTGGGCGGATCCCCACCAG*

When I blat this on UCSC Genome Browser, the first hit is a match which spans 6572 basepairs on chromosome 19, and essentially matches the first portion of the sequence to one spot on the genome, and the second portion to another part of the genome approximately 6.5kb downstream.

However when I use the BLAT command line tool, it doesn't find the second half of this match, it only matches the first portion.

Here are the flags I'm using:

**-t=dna**

**-q=dna**

**-minIdentity=90**

**-out=blast8**

How do I get the command line tool to mimic the behavior of BLAT on UCSC Genome Browser?

# Answer
> 5 votes
The FAQ offers an answer:

> **I'm setting up my own Blat server and would like to use the same parameter values that the UCSC web-based Blat server uses.**
>
> We almost always expect there to be some small differences between the hgBlat/gfServer and the stand-alone command-line blat. The best matches can be found using pslReps and pslCDnaFilter utilities. The web-based blat is tuned permissively with a minimum cut-off score of 20, which will display most of the alignments. Other than to confirm that your command-line blat is working, there is little use in perfectly replicating the web-based blat results. We advise deciding which filtering parameters make the most sense for the experiment or analysis. Often these settings will be different and more stringent than those of the web-based blat. With that in mind, use the following settings to replicate the search results of the web-based blat:
>
> \[ . . . \]
>
> standalone blat:
>
> blat search:
>
> ```
> blat -stepSize=5 -repMatch=2253 -minScore=0 -minIdentity=0 database.2bit query.fa output.psl
>
> ```

---
Tags: human-genome, genomics, blat
---


# Question
Title: Error rate setting in Canu error correction
I want to use Canu to correct my nanopore long read (version: MinION R9.5), but I am not quite sure how to set the correctErrorRate. Should I follow the Canu manual (Nanopore R7 2D and Nanopore R9 1D Increase the maximum allowed difference in overlaps from the default of 14.4% to 22.5% with correctedErrorRate=0.225), or you guys have a better option?

# Answer
It's generally a good idea to trust the "official" suggestions. You can also adjust the error rate based on coverage according to parameter reference:

> For low coverage datasets (less than 30X), we recommend increasing correctedErrorRate slightly, by 1% or so.
>
> For high-coverage datasets (more than 60X), we recommend decreasing correctedErrorRate slighly, by 1% or so.
>
> Raising the correctedErrorRate will increase run time. Likewise, decreasing correctedErrorRate will decrease run time, at the risk of missing overlaps and fracturing the assembly.

That being said, the sequencing quality can vary from between libraries or flow cells. The error rate is not a constant.

According to the developer (full thread):

> if you want to assemble 1D data or have a bad run, it is possible you would need to increase the error rate. You could run with a high error rate (0.1) and look at the distribution of overlap error rates in the unitig step to look for a peak in the distribution. If you have a near neighbor you could also map the corrected reads to it and estimate the residual error that way.
> 3 votes

# Answer
There seems to be some confusion here.

For the first stage of the Canu pipeline (read correction), i think the parameter to use is **rawErrorRate**.

From the docs:

**rawErrorRate** The allowed difference in an overlap between two uncorrected reads, expressed as fraction error. Sets corOvlErrorRate and corErrorRate. The rawErrorRate typically does not need to be modified. It might need to be increased if very early reads are being assembled. The default is 0.300 For PacBio reads, and 0.500 for Nanopore reads.

I think that **correctedErrorRate** is what you use at the next stages of the pipeline (trimming and assemble).
> 0 votes

---
Tags: nanopore, canu, quality-control
---


# Question
Title: How to determine the primary Uniprot accession number from a set of accession numbers?
Given a list of Uniprot IDs that are linked to an Ensembl gene ID, is there a way to systematically determine the UniProt accession corresponding to the protein product of the gene's primary transcript **with no other information**?

According to ExPasy

> Researchers who wish to cite entries in their publications should always cite the first accession number. This is commonly referred to as the 'primary accession number'. 'Secondary accession numbers' are sorted alphanumerically.

But what if the order has been jumbled or compiled from a different source and resorted?

For example:

Ensembl:

ENSMUSG00000035642

Uniprot:

Q8R0P4, Q8CF11, D6RJK8, D6RJJ4, D3Z442, D3Z1Q3, D3YZD8, D3YY39, D3YX09, D3YWY5

This question is cross listed on the Biology stack-exchange.

# Answer
> 4 votes
Reposting my answer from the related ticket in the Biology section:

I think there is an issue with the terminology. The "primary" accession number, is the first accession number in cases where an entry has more than one accession number, as described in http://www.uniprot.org/help/accession\_numbers:

> Entries can have more than one accession number. This can be due to two distinct mechanisms:
>
> ```
> a) When two or more entries are merged, the accession numbers from all entries are kept. The first accession number is referred to as the
>
> ```
>
> ‘Primary (citable) accession number’, while the others are referred to as ‘Secondary accession numbers’. These are listed in alphanumerical order.
>
> ```
> b) If an existing entry is split into two or more entries (‘demerged’), new ‘primary’ accession numbers are attributed to all
>
> ```
>
> the split entries while all original accession numbers are retained as ‘secondary’ accession numbers.
>
> Example: P29358 which has been ‘demerged’ into P68250 and P68251.

Both reviewed and unreviewed entries can have primary accession numbers.

What you probably mean, as previous posters understood, are accession numbers of reviewed entries as opposed to unreviewed ones.

In that case, you can indeed add "reviewed:yes" to your query, e.g. when you are using the UniProt ID mapping, http://www.uniprot.org/help/uploadlists

# Answer
> 3 votes
Consider the circumstances in which entries will have more than one accession (taken from the same page you linked in the original post).

> Entries will have more than one accession number if they have been merged or split. For example, when two entries are merged into one, the accession numbers from both entries are stored in the AC line(s).
>
> If an existing entry is split into two or more entries (a rare occurrence), the original accession numbers are retained in all the derived entries and a new primary accession number is added to all the entries.

As mentioned in the comments, this *does not* imply that there is or should be a single "primary" UniProt accession associated with (for example) an Ensembl gene ID. In some cases, a single protein record may have multiple accessions due to annotation changes, and the suggestion from your original post is to only report the primary accession. However, there may be multiple proteins associated with an Ensembl gene ID, and each protein will have its own accession (or set of accessions).

Choosing a single representative isoform for a gene is an entirely different question.

# Answer
> 0 votes
Not with no other information, no. There is no way of looking at a list of accessions and guessing which one is primary. Primary accessions are just accessions like any other; they don't have any specific characteristics.

This is further complicated by the fact that each gene can (and usually does, for higher eukaryotes) have multiple transcripts, each of which can produce a different protein product. So in order to go from gene identifier to a protein, you need to first find the relevant transcript.

So, as a first step, get the primary transcript IDs for each of your genes and then retrieve the protein accessions associated to them.

---
Tags: proteins, public-databases, identifiers
---


# Question
Title: Server for finding kmers in set of sequences
Is there a server/website somewhere where I can submit a list of DNA/RNA sequences and find the list of kmers hits and organisms where it's found? I checked the kraken website, but they don't have a webserver for it.

I am specifically looking for kmers, as a complement of just doing BLAST on the NCBI nt database or similar.

EDIT: So far I have seen:

* the Kraken app on Illumina BaseSpace, but it's not trivial to upload FASTQ data that hasn't been produced from an Illumina instrument.
* the CLARK trio of apps in InsideDNA, tried once with default parameters but produced empty results.

# Answer
One Codex is powerful and user friendly, and works in a similar way to Kraken underneath https://www.onecodex.com/

The Kaiju web server might be of interest and is more sensitive than exclusively k-mer based methods http://kaiju.binf.ku.dk/server

The title of this question is misleading – it sounds like you're actually looking for a webserver for taxonomic classification using *k*-mers.
> 6 votes

---
Tags: k-mer
---


# Question
Title: SNP-phenotype association analysis only using the SNPs on a specific gene
I know that the GWAS association p-value threshold is 1e-8. This makes sense because in GWAS, you make ~1 million hypothesis (i.e. use that many SNPs in the association test). However, let's say I do not do genome-wide analysis but only use the 100 SNPs from a specific gene to check for association with my phenotype. Then, do you think a raw p-value of 1e-4 (which becomes 1e-2 after multiple hypothesis correction) is a reliable p-value to come to a conclusion that there is an association between that SNP and the checked phenotype? Obviously, if I checked all 1 million SNPs instead, this p-value would be insignificant after correction, but I am not checking 1 million SNPs, as I mentioned.

# Answer
> 5 votes
> I know that the GWAS association p-value threshold is 1e-8

This may be a common threshold of statistical significance that is used, but it's definitely not an absolute value. It's a hack to try to work around many issues with GWAS associated with testing millions of SNPs. Unfortunately, the most relevant issue in GWAS (for spurious significance) is unexpected shared genetic structure in the cases or controls, and this cannot be excluded by a p-value threshold, regardless of how low it is set.

In your 100-SNP example, you are correct that the bonferroni correction would be a threshold of 1e-4 (for a desired significance threshold of 0.01), but it's important to be cautious about GWAS results, even when the p-value indicates otherwise. Trust no one, especially yourself. P-values should never be used as an indication of importance (which suggests that most GWAS results should probably be reconsidered), and are better used in conjunction with other evidence for identifying relevant / significant tests.

My recommendation for GWAS (or any similar bulk genotyping study, as in your 100-SNP subset) is that at least a bootstrap sub-sampling process is carried out, to make sure that observations found within the comparison of interest are at least well replicated when comparing smaller sub-groups of the cases and controls. More information on that can be found in the poster I presented at Queenstown Research Week last year, and detail about the implementation in my draft preprint (which will probably never be properly published, given how old the research is).

---
Tags: gene, snp
---


title: Clarification on Gene Enrichment

text:
# Question
Title: Clarification on Gene Enrichment
When I run a GSEA analysis on two conditions from the same RNaseq (negative control PBS injection VS positive control CpG injection) from the same dataset/same gene list, I get results that look something like this:

Notice in my example that many gene sets are significantly enriched in PBS VS CpG, but none are very significantly enriched in the inverse, CpG VS PBS.

I have a basic question:

If we're comparing two items, say A and B, shouldn't a up-regulated gene set that's enriched in A have a corresponding down-regulated gene set that's enriched in B?

**I'm consistently confused about how to interpret the fact that one of my two conditions has many more enriched gene sets than the other.** What am I missing/misunderstanding?

Thanks!

# Answer
If you compare A vs B the genes's fold change will have the opposite sign to B vs A. So will be the gene set up or down-regulated depending on the comparison to take.

The gene set test analyze if a given group of elements is sorted in a certain way in a list (I am talking about the GSE like the one you performed or the one in the Broad Institute). Usually it is used to say that gene set X is up-regulated in the comparison A vs B, that is, gene set X (in general) is more expressed in A than in B. The enrichment is measured in a enrichment score (or a normalized enrichment score, the higher, the clearer trend/distribution group X has. To asses X is really differentially expressed a p-value is also calculated.

Despite that the number of the up-regulated gene sets sum up all the gene sets tested as pointed by juod, the number of gene sets with a a p-value below 5% and 1% is not the same in both comparisons which indicates that you are not really using the same samples in both gene sets.

Check how you performed your comparison. Are the genes in each comparison with the same value and opposite sign?
Check also the software used to run the GSEA analysis (if I recongize the type of report is from a package in Bioconductor \[I don't recognize which one now\]) as you (or me) might have misunderstood something or there is an error if the comparisons are done properly. For instance it could be that the report of p-values below the threshold would be from only the up-regulated gene sets.
> 8 votes

# Answer
While I'm not familiar with GSEA software in particular, I believe your problem is that it only tests for *upregulated* gene sets. Notice:

> * 866/4408 gene sets are upregulated
> * 3542/4408 gene sets are upregulated

3542+866 = 4408. I.e., 866 sets have higher mean expression in positive condition, the rest have higher mean expression in negative condition. To avoid this confusion, I personally would reserve using the terms "up/down-regulated" only for sets where this change is at least remotely significant.
> 5 votes

---
Tags: rna-seq, gsea, go-enrichment
---


title: Ultimate reproducibility in R?

text:
# Question
Title: Ultimate reproducibility in R?
I'm looking for a convenient and reliable way to make an R analysis reproducible, either at different times or across collaborators.

Listing the package versions or a `sessionInfo()` output is not very helpful, as it leaves the work of re-creating the environment entirely to you.

As I see, in order to make analysis reproducible, the following has to be the same:

* the **R environment** (R, repos (Bioconductor), packages, and maybe RStudio)
* **data**
* **code**
* **seed** for random processes: `set.seed()`

I'll post an answer, but would like to hear about other options.

# Answer
You're right about trying to get the R environments to be as similar as possible, as easily as possible. It's not always necessary to have the exact environment to get the same results, but it definitely helps if you find yourself struggling to replicate results.

# General tips for ensuring reproducibility, not specific to R

See Sandve *et al.*, PLoS Comp Bio, 2013. They mention a number of things that you do, like recording seed values and session info. I'll highlight a few points:

1. **Environment**: Include a small script that lists all the necessary dependencies
    * This should include all version numbers of languages and packages used
    * *Better implementation*: make it a script that automatically installs dependencies
2. **Data**: Share raw and/or processed data where possible
    * At the very least, details on how to obtain the data you used
    * This can be an accession number in a database, or a link to the source material
    * *Better implementation*: again, make it a script that automatically downloads the raw data you used
3. **Code**: Comment and document your code
    * People (yourself included) need to understand what you've done and why
    * Additionally, your work may be reproducible, but it may not be correct. If you (or anyone else) attempts to fix a problem in your code, they will have a difficult time doing so if they can't understand what your code is doing
4. **Code**: Use version control
    * Services like GitHub and BitBucket make it simple to share your exact work with collaborators
    * Version control also ensures that you have a history of all the changes done to your work
5. No manual manipulation of data
    * If something is done to your data, write it down

# R-specific tips

1. The `packrat` package, as you mentioned in your answer
2. Using a number of custom functions for analyzing your data? Compile those functions into an R package
    * Allows others to easily load/unload the exact functions you used through tools like CRAN, BioConductor, or `devtools::install_github()`
    * Makes it easy to see whether these functions are working as intended, since they are separated from other parts of the analysis
3. Make unit tests for your work with the `testthat` package
    * Helps ensure that your code is working as intended, and gives others confidence in your work as well
    * Write the tests using `testthat::test_that()`, and check that they pass using `R CMD check` or `devtools::check()`
    * This works especially well with R packages
    * This also works especially well with TravisCI, since you can automatically have your code get tested after you change it and make a commit
4. Consider Docker images for your environment
    * This bundles everything about the R environment together into one easily executable image (and it's not just limited to R)
    * These images can be easily shared through Docker Hub
    * **Note** this might not always be the best solution, since few people know Docker and it's not always simple to use. Docker also requires admin privileges, which may not be possible for work on computer clusters
    * There is lots of current discussion about Docker in science. See Boettiger, 2014, for an example, using R
> 16 votes

# Answer
Another option I'm trying at the moment is the combination of snakemake, conda.io and knitr.

* **Conda** is a "Package, dependency and environment management for any language" and could be used to set up the **R-environment** as well as other software. Note, most bioinformatics software can be found in the related **bioconda** project.
* **Snakemake** is a workflow management system, which can be used to automate the preprocessing of your data in a reproducible manner (also it works nicely with conda since a conda environments can be specified in the Snake-Makefile)
* Knitr for dynamic report generation with R. Well suited to document and reproduce your R workflow into reports, ticking your **code** and seed **boxes**.
> 5 votes

# Answer
One way to achieve this is via `packrat`, a dependency management system developed by RStudio.

Setup:

1) Create project

2) Create script in project directory to start RStudio with specific R version:

```
export RSTUDIO_WHICH_R=/path/to/R
rstudio

```

3) Bring project under `packrat` control. Install bioconductor, install packages...

---

Every time you work on the project, run the above script in project directory. https://github.com/rstudio/packrat/issues/342#issuecomment-315059897
> 1 votes

# Answer
*checkpoint* is a system for ensuring you get packages of the same time / version. However, I confess that I haven't used it and thus can't vouch for it.

My feeling is that the rest will have to be handled outside this: code and data can be versioned in git (and dat). A whole solution might have use something like docker as others have suggested.
> 1 votes

---
Tags: r, reproducibility, workflow-management
---


# Question
Title: Melt p-values for CpG sites mapping to the same gene
I have some data I am working with, and I am curious if I am able to combine p-values from a paired t-test for CpG sites in the genome using Fisher's Method to get one p-value for each unique gene. Linked here is the Wikipedia page for Fisher's Method. I understand that an assumption of the method used is that each individual p-value being combined must be independent. I'm relatively new to biostatistics, so I'm curious if using CpGs from the same gene would violate this assumption.

# Answer
> 6 votes
Methylation levels have high local correlation, so Fisher's method would be problematic. Having said that, you have no reason to use Fisher's method after a paired t-test. A paired t-test will give you a single p-value per gene, which is what you want. Do be sure to only include CpG with some minimal coverage in both group.

# Answer
> 1 votes
What is typically done in methylation analysis is to assess the islands of methylations.

Check this workflow, in the section linked it uses some predefined islands for instance. I am no expert on this area but you could asses if the islands or certain regions are more methylated than expected.

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


title: Dealing with absence of coverage

text:
# Question
Title: Dealing with absence of coverage
In my snakemake workflow, I use deeptools `bamCoverage` (in a `make_bigwig` rule), then `computeMatrix` and `plotProfile` (in different `plot_*_profile` rules) to create coverage profiles for various kinds of small RNAs. It sometimes happens that there are very few reads in a particular category of small RNAs. This results in `bamCoverage` exiting with an error (when I use the `--skipNAs` option, in order to make coverage profiles not take into account empty regions):

```
Error: The generated bedGraphFile was empty. Please adjust
your deepTools settings and check your input files.

```

The fact that `bamCoverage` fails make things a little annoying for me, because the later the failure, the easier to handle things in my workflow: I will need to catch the error and generate a dummy output file in the `make_bigwig` rule, and then handle the special dummy file in the `plot_*_profile` rules.

I would be fine having an empty coverage file: the absence of reads (either in the bam file, or matching selection criteria in `bamCoverage`) is not due to an error.

**Is this just not possible, by the format definition?**

**Did I miss an option to have `bamCoverage` accept making empty files?**

**Would this make `computeMatrix` fail?**

# Answer
I'll fix the `bamCoverage` issue in the 2.6.0 release (it may already be fixed since I rewrite how bigWig files are made for that release). You can follow my status on the issue page. In principle you should be able to do this, it's just a deepTools limitation.

This should work fine with `computeMatrix`, though keep in mind that you're just going to get a big block of NA values.

**Update**: This error ends up coming from pyBigWig. I've just found that `bedGraphToBigWig` also can't make a bigWig file with now entries, so I suspect this is sort of undefined behavior. Regardless, I'll alter pyBigWig to produce a file that it can open and use in a meaningful manner.
> 3 votes

---
Tags: snakemake, deeptools, workflow-management, errors
---


# Question
Title: How do I efficiently perform a metagenome screen of “all” species?
I’ve got an RNA-seq dataset with a large proportion of environmental RNA “contamination”. BLASTing random reads reveals that much of the data comes from bacterial, plant and viral RNA. My target organism only accounts for ~5% of the RNA-seq read data.

I would like to get a more or less (more, if possible) comprehensive overview of what species are found in my sample — bacteria, plants, animals (?), but also viruses. How can I perform this?

I have been unable to find a “standard” way of performing a metagenome screen. BLAST online services all seem to be severely rate limited (certainly unable to handle uploading ~80 M reads). Installing BLAST (or a similar tool) locally is of course not an issue, but I can’t find a comprehensive database that spans all phyla — the best I’ve been able to find are databases, such as NCBI-NR, that are restricted to single phyla or classes, or to bacteria.

Ideally I’d like a workflow that I feed with my RNA-seq data, and get out a list of species with % of reads mapped.

# Answer
> 7 votes
As you point out, BLAST is not going to work well for reads from a typical RNASeq experiment.

From an RNASeq dataset, you might be able to get away with just a ribosomal screen using Silva, because the ribosome is the most abundantly transcribed gene. This will usually be the case even if polyA selection or ribosomal depletion was carried out on a sample.

However, bear in mind that there will be a lot of false results due to database incompleteness. I used Kraken when I did something like this a few years ago, and ended up with one sequence annotated as a fruit bat (where my subject was a water-living planarian flatworm). Centrifuge would be a better option now; the authors even claim that they can store the entirety of nr in a 70 GB index with Centrifuge (which would be better / more comprehensive than just a ribosomal screen).

# Answer
> 4 votes
I think you're just looking for `nr`. It is absolutely not limited to prokaryotes, far from it. It is, according to the blurb on the NCBI's blast page:

> The nucleotide collection consists of GenBank+EMBL+DDBJ+PDB+RefSeq sequences, but excludes EST, STS, GSS, WGS, TSA, patent sequences as well as phase 0, 1, and 2 HTGS sequences. The database is non-redundant. Identical sequences have been merged into one entry, while preserving the accession, GI, title and taxonomy information for each entry.

The link you gave in your question doesn't describe the `nr` database, but the collection of RefSeq non-redundant sequences. The two are not the same thing at all. The nr databases contains sequences from all available species so it should be perfect for what you're after.

Now, I admit I haven't tried this with millions of redas, but I did used to run many blast searches with several thousand sequences using the command line blast tools and the remote nr database. I just checked the blast+ docs now and see no mention of a rate limit. Which isn't to say there isn't one, for sure, but since each query is run sequentially when you give a multifasta file, there seems to reason to limit since the servers should be able to queue the jobs as needed.

So, don't use the API or roll your own, just download NCBI's blast client and try that:

```
blastn -query query.fa -db nr -task blastn

```

If that complains, you'll have to download the `nr` DB (well, since you want nucleotide blast, you probably want its nucleotide cousin, `nt`) which you can find here: ftp://ftp.ncbi.nlm.nih.gov/blast/db/. Just get all the `foo.nt*` files.

# Answer
> 4 votes
You might want to look into minihashes as a more efficient way of clustering and comparing your sequences to known species and samples. In the mash paper the authors show how it would be feasible to create minihash sketches that you could search against from large metagenomics datasets (check the "Clustering massive metagenomic datasets" section in the results):

> For a large-scale test, samples from the Human Microbiome Project \[36\] (HMP) and Metagenomics of the Human Intestinal Tract \[37\] (MetaHIT) were combined to create a ~10 TB 888-sample dataset. Importantly, the size of a Mash sketch is independent of the input size, requiring only 70 MB to store the combined sketches (s = 10,000, k = 21) for these datasets.

Moreover, the creators of sourmash (a tool using the same approach) provide precomputed sketches for the whole of NCBI's RefSeq database (includes bacterial, viral and fungal genomes). Creating a similar sketch for the whole of ENA/SRA (or a subset of it) could require some time but would allow for very fast and accurate searches.

# Answer
> 4 votes
One thing that hasn't been suggested so far would be phylogenetic placement. First option there would be placing against one larger tree roughly touching the domains you mentioned to get an overview. Second option would be having three trees (bacteria, plants, viruses) made up of relevant taxa for the type of sample you have.

This would give you a detailed and accurate view of your sample, but it does involve some work. Sadly the computational throughput for placement is still a bit behind BLAST, and not many complete pipelines exist, but in my opinion it is worth the extra mile (I might be biased. I'm working to resolve as many of these issues as possible).

If you want to look into it, I can reccomend this excellent paper by Erick Matsen. Since your samples sound vaguely soil-ish, I will also link this recent Nature paper using placement to analyze rainforest soil samples.

# Answer
> 1 votes
I have had good success using the MEGAN software to detect and visualize contamination. My input data was WGS not RNA-seq But I would check it out. It can create a very nice weight phylogenetic tree for all the detected organisms.

https://ab.inf.uni-tuebingen.de/software/megan6

---
Tags: rna-seq, metagenome
---


# Question
Title: Bacterial genome annotation of a clinical isolate strain?
So I'm basically so new to this that I'm just trying to find out what tools, methods, and keywords I should go look up by myself.

I have a unique strain of a bacteria.

I was given RNAseq data for this unique strain.

I want to analyze the RNAseq data, but I need an annotated genome to do that.

As a resource, I do have the annotated genome of the common lab strain.

The common strain should be pretty similar to the unique strain, but I don't know how similar.

Can anyone tell me what tools/methods might be applicable so I can go look them up?

# Answer
For a quick (but reliable) analysis, I'd recommend using Kallisto or Salmon to quantify isoform read counts using the transcriptome of the lab strain.

If you have a concern about transcripts that are in your sample but not in the lab strain, you can do two Trinity assemblies: one fully *de-novo*, and one genome-guided. These transcripts can then be used as a reference assembly for Kallisto/Salmon, and counts compared to the lab-strain results. Note that the Trinity assemblies will only apply to the particular environments that are in your samples; if a transcript isn't expressed, it won't be in the assembly.

Any transcripts in the *de-novo* assembly that are not in the genome-guided assembly are potentially novel to your strain. However, care must be taken in interpreting this: if there is any other contamination in the sample, those contaminant transcripts will also be assembled.

I'm not too familiar with annotation pipelines, but based on a couple of Twitter posts it looks like Prokka might be a reasonable start. Torsten Seeman (lead developer of Prokka) has done a post about alternative annotation pipelines:

http://thegenomefactory.blogspot.co.nz/2013/03/bacterial-genome-annotation-systems.html

Update: Torsten Seeman now recommends Bakta as a "worthy successor to Prokka":

https://github.com/oschwengers/bakta

> Bakta is a tool for the rapid & standardized annotation of bacterial genomes and plasmids from both isolates and MAGs. It provides dbxref-rich, sORF-including and taxon-independent annotations in machine-readable JSON & bioinformatics standard file formats for automated downstream analysis.
> 5 votes

---
Tags: rna-seq, sequence-annotation, reference-genome
---


# Question
Title: Recommendations for missing value imputation - DNA methylation data
I'm looking for some options for imputation for a high-dimensional dataset of DNA methylation (bisulfite sequencing) data. Dimensions on the order of 50-100 samples x ~500,000 CpG loci/features.

I've used K-nearest neighbors, but It seems that this method is not very accurate. As far as I can tell, it limits the minimum number of "genes/features" to impute at one time to something like 1500. which is kind of a lot.

I've also used missForest to impute smaller datasets with greater accuracy, but it seems computationally infeasible to do it on the full 50x500,000 dataset.

I need to impute values because I'm using the data for some statistical modeling which require complete cases.

Anyone know of better alternatives than K-nearest neighbor for large scale imputation?

# Answer
> 2 votes
I think this is still an active field of research. I have heard of Phenix, which might be appropriate.

---
Tags: methylation, sequence-analysis, imputation
---


# Question
Title: What is the best distribution to model the FPKM values from normalized RNA-Seq data?
I know that the discrete raw counts from the RNA-Seq data are usually modeled by a negative binomial or a Poisson distribution, but what I am working on are the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values which represent a transformation of the raw values. They are continuous except many actual 0s. What is the probability distribution FPKM values are likely to follow, or are assumed to follow? Even after taking logarithm twice (after replacing 0s with 1, of course), the data does not follow a normal distribution since there are a lot of 0s.

# Answer
FPKM/TPM values are generally log-normally distributed. Reference : Gene expression profiling in single cells from the pancreatic islets of Langerhans reveals lognormal distribution of mRNA levels
> 5 votes

# Answer
Cuffdiff is one the tools that deals FPKM data for analysis of differential expression.By default they use the negative binomial with some condition specific parameter to avoid biases. There are some other models also available in Cuffdiff. Check the following documents for more details,

http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/index.html#library-types

https://www.reddit.com/r/bioinformatics/comments/3bx3em/fpkm\_vs\_raw\_read\_count\_for\_differential/
> 0 votes

---
Tags: rna-seq, statistics
---


# Question
Title: Retrieving a list of human genes having GO associations
I need a command in R to retrieve all human genes associated with a Gene Ontology entry. I tried to look for it online but did not find it.

# Answer
Here's an example for the mouse genome:

```
library(org.Mm.eg.db)
select(org.Mm.eg.db, c("GO:0048406"), c("GENENAME","SYMBOL"), c("GO"))

```

You get output like:

```
                                                    GENENAME SYMBOL
1                                     pregnancy zone protein    Pzp
2 nerve growth factor receptor (TNFR superfamily, member 16)   Ngfr
3 nerve growth factor receptor (TNFR superfamily, member 16)   Ngfr
4                                             neurotrophin 3   Ntf3
5             neurotrophic tyrosine kinase, receptor, type 1  Ntrk1
6            furin (paired basic amino acid cleaving enzyme)  Furin
7              proprotein convertase subtilisin/kexin type 6  Pcsk6
8                                                 sortilin 1  Sort1
9                                      alpha-2-macroglobulin    A2m

```

**Edit**: Given your comment that you want instead a list of all genes with a GO term associated with them (for human):

```
library(org.Hs.eg.db) # Install it from bioconductor yourself
unlist(as.list(org.Hs.egSYMBOL)[mappedkeys(org.Hs.egGO)])

```

The `mappedkey()` part gets the index into the bimap of all entries with a valid GO mapping. That's then used to subset the symbol bimap (after converting to a list).
> 6 votes

# Answer
I guess the following code will help,

```
source("https://bioconductor.org/biocLite.R")
biocLite("biomaRt")
library("biomaRt")
ensembl = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
#listAttributes(ensembl)
mapping <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol", "go_id"), mart = ensembl)
head(mapping)

```
> 6 votes

---
Tags: r, human-genome, gene-ontology
---


# Question
Title: Is it possible to do alignment within Python? Check variants against reference?
I'm currently adding a few SNPs randomly into a FASTA within python using BioPython. In the following example from BioPython, I add an SNV at location "5"

http://biopython.org/DIST/docs/api/Bio.Seq.MutableSeq-class.html

```
from Bio.Seq import MutableSeq
from Bio.Alphabet import generic_dna
my_seq = MutableSeq("ACTCGTCGTCG", generic_dna)
my_seq[5] = "A"
print(my_seq)
## outputs ACTCGACGTCG

```

In my case, I'm using the reference hg38 as the input FASTA to manipulate.

After randomly inserting these SNPs, I would like to find a region whereby there are no SNPs within (let's say) +/- 50Kb, or the next nearest thing.

Is it possible to align the FASTA with variants against the reference within python? Otherwise, is there a quick way to check the variant FASTA against the reference FASTA, finding a region with no variants? What would be an efficient way to execute this task?

# Answer
> 3 votes
As others have said, doing whole-genome (or even whole-chromosome) alignments is the wrong solution. Simply track where you're creating SNVs. If you wrote your SNV locations to a sorted BED file you could use something like the following:

```
#!/usr/bin/env python
minWidth = 100000  # Only report regions of >=100kb
last = [None, None]
for line in open("some bed file"):
    cols = line.strip().split()
    if cols[0] != last[0]:
        if int(cols[2]) >= minWidth:
            print("{}:1-{}".format(cols[0], cols[2]))
    else:
        if int(cols[2]) - last[1] >= minWidth
            print("{}:{}-{}".format(cols[0], last[1], cols[2]))
    last = [cols[0], int(cols[2])]

```

This will miss the last part of each chromosome/contig, since it doesn't know how long they are, but you should be able to figure out how to handle that. The memory requirements of this are trivially small.

I suspect there's also a bedtools or bedops command that will directly give you the distance to the next downstream region. One could alternatively parse that (it's all the above code is doing, but using `bedtools` or `bedops` would produce more validated results than a quickly thrown together script).

# Answer
> 2 votes
If you store the positions of your random SNPs in a `sort-bed` sorted BED file, you can filter with BEDOPS `bedmap`:

```
$ bedmap --count --echo --range 50000 --delim '\t' SNPs.bed | awk '$1==1' | cut -f2- > SNPs.filtered.bed

```

The file `SNPs.filtered.bed` will only contain those elements from `SNPs.bed` which do not overlap within 50kb.

Or you can use BEDOPS `bedops` and `closest-features` to do the same thing, but the logic isn't nearly as pretty:

```
$ bedops --everything --range 0:50000 test.bed | closest-features --dist test.bed - | awk -v FS='|' '{ if ($NR > 50000) { print $1; } }' > SNPs.filtered.bed

```

The result is the same as with the `bedmap` pipeline. IMO, the `bedmap` approach is cleaner.

Once you have `SNPs.filtered.bed`, you could convert it to FASTA via `samtools faidx` and `bed2faidx.pl` or similar. Then you have FASTA upon which you can apply your random perturbations.

---
Tags: sequence-alignment, python, reference-genome
---


# Question
Title: What is the R/Bioconductor preferred way of accessing TCGA data?
It seems there are two packages providing access to TCGA data these days:

Will *GenomicDataCommons* supersede *TCGAbiolinks*? Is there anything that I might miss if I decide to use only *GenomicDataCommons*?

# Answer
There isn't **a** R/Bioconductor person or view, so there isn't a preferred way to access TCGA. Use what suits you.

Bioconductor is build by contributors which are encourage to not overlap in their contributions. Usually the packages have different options/features even if they allow to access to the same data.

In this case, the TCGAbiolinks package is specifically for TCGA dataset while GenomicDataCommons is for "access the NIH / NCI Genomic Data Commons RESTful service". So the functionality they bring to R and Bioconductor is different.

The fact that one project is finishing doesn't mean that the data or the packages build around that project are obsolete. Some functionality different between them is the functions to harmonize and prepare (`GDCprepare`) data from the TCGAbiolink package, or see a boxplot correlation and AAIC plot to define outliers `TCGAanalyze_Preprocessing`.
> 1 votes

---
Tags: r, bioconductor
---


# Question
Title: Can FASTA files have nucleotide and protein sequences within them; or must they only have 1 type?
Can FASTA files have nucleotide and protein sequences within them; or must they only have 1 type? For example, a FASTA file has 2 sequences. Can the first one encode amino acids while the second one encodes bases?

Thank you

# Answer
While there's nothing stopping anyone from doing that with the FASTA format (after all, it's just a text file with '\>' defining header lines), I don't know of any software that would support such a file structure. At best, it would interpret the nucleotide sequences as protein sequences (A/C/G/T are all valid 1-letter protein codes).

A better question to ask would be "Does ultra-cool-bioinformatics-tool X support combined nucleotide and protein sequences in the same FASTA file?" In which case the answer would most likely be, "No."
> 9 votes

---
Tags: fasta
---


# Question
Title: Secretome and membrane receptor profiles
I am looking the secretome profile and the membrane receptor profile for a given cell type.

In my specific case, this should be the secretome and outer membrane receptor profiles of dorsal root ganglion.

What I've done in the past is taken proteomics datasets or RNAseq datasets and used this as a reference cell surface proteome or secretome.

I wasn't sure what the best way to filter secreted proteins or membrane receptors were, so I used the Human Protein Atlas and this turned out to be a big headache. A possible alternative method might be to use DAVID, but there must be a more efficient method via the command line.

My questions are these:

* Is there a "consensus" database containing protein expression profiles for different tissue types in human?
* If not, what is the best way to filter out the **secreted** and **transmembrane** proteins by gene symbol?

# Answer
> 1 votes
I don't think there is a database with expression data for your particular cell type. I am afraid your "headache" approach doesn't sound so bad...

So get expression data e.g., from GEO for your cell type, and define which genes/proteins you consider as "expressed" (e.g., minimal RPKM value).

Then, you can use GO annotation to find your genes involved with membrane or secretion (or what ever). You can use biomaRt in R for this.

For example, you want to know all "membrane" genes.

```
> library(biomaRt)
>
> ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
>
> membrane.genes <- getBM(attributes=c('hgnc_symbol',
> 'ensembl_transcript_id', 'go_id'),
> filters = 'go', values = 'GO:0016020', mart = ensembl)

```

And now you want to know it only from genes that are expressed in your ganglion cells (with dummy example of two genes).

```
> gene.list <- c("KIR2DL3", "GPR1")
>
> membrane.genes.list <- membrane.genes[membrane.genes$hgnc_symbol ==
> gene.list,]
>
> membrane.only <- membrane.genes.list[membrane.genes.list$go_id ==
> 'GO:0016020',]

```

# Answer
> 1 votes
> Is there a "consensus" database containing protein expression profiles for different tissue types in human?

https://www.ebi.ac.uk/gxa/home

> If not, what is the best way to filter out the secreted and transmembrane proteins by gene symbol?

Personally I like to run a local copy of SignalP, as it will yield consistent results. However bioinformatic approaches tend not to capture highly unusual transmembrane proteins or membrane associations (e.g.: only one leaflet etc.). http://www.cbs.dtu.dk/services/SignalP/

Depending on the organisms, and scope of your analysis, further possibility could be to visit supplemental tables of recent publications containing somewhat supervised computational predictions on secretion, e.g.: http://science.sciencemag.org/content/356/6340/eaal3321

> I am looking the secretome profile and the membrane receptor profile for a given cell type.

The above links don't inform about cell-type specific retention of potentially secreted proteins. - but their combination will yield cell-type specific RNA and protein levels of transmembrane and secreted proteins.

# Answer
> 0 votes
You can get expression data from databases and for membrane proteins, run a transmembrane protein prediction method like CCTOP (http://cctop.enzim.ttk.mta.hu/).

---
Tags: public-databases, networks
---


# Question
Title: Why is this makeblastdb command not working?
I am trying to make local databases for the ncbi-blast+ package (version 2.60). I am doing so for 4 T-cell receptor genes. 3 of the 4 (TRAV, TRAJ, TRBV) have worked fine, but I am having problems with TRBJ. I say this because they should be the same format and if one works the rest works.

They take the format below:

```
>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15| | |
NTEAFFGQGTRLTVV

```

I used the following command (like I say, it was fine for the rest, oddly enough):

```
makeblastdb -in TRBJ_AA.fa -parse_seqids -dbtype prot

```

and I get the following error message:

```
BLAST options error: TRBJ_AA.fa does not match input format type, default input type is FASTA

```

One issue I checked was the length of the sequences being too short, but I believe 11 is the minimum length required (for inputs, don't know if that is the same for databases)- but would someone be able to confirm, that strikes me as the most likely candidate?

I'm happy to use some other metric, but I was hoping on keeping BLAST for consistency's sake for the time being.

Thank you.

edit: I know it's to late, but for reference I'm on Ubuntu 16.04.2 LTS x86\_64 for all steps.

# Answer
This looks like a bug in `makeblastdb`. Removing the final `| | |` from your sequence's description makes it work:

```
>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48 nt|3| | | |15 AA|15+0=15

```

So does removing everything after the 1st space:

```
>K02545|TRBJ1-1*01|Homo sapiens|F|J-REGION|749..796|48

```

There is nothing in the FASTA definition that would explain this and, even more weirdly, the presence of the `| | |` alone isn't enough to trigger the issue. I played around with this a little and it looks like it doesn't like the combination of `+` or `-`along with empty `| |`, but there's no clear pattern I can see:

These fail:

```
>| | |15-| | |
>| | |15+| | |
>| | |15-1| | |
>| | |15+1| | |

```

These work:

```
>| |15+| |
>| | |151| | |

```

So yeah, sounds like a bug. I suggest you report it to the NCBI devs. In the meantime, a simple fix would be to just keep the first space-separated field of your fasta header. Pass your fasta file through this:

```
awk '{print $1}' fasta.fa > new.fa

```

That will print the 1st field of each line and, assuming you don't have spaces in your sequence, it should effectively sort this out for you.
> 8 votes

---
Tags: fasta, blast
---


title: hg38 GTF file with RefSeq annotations

text:
# Question
Title: hg38 GTF file with RefSeq annotations
I'm not sure what I'm missing, but I'm struggling to find an official `hg38` `GTF` file with RefSeq annotations. I'd like to provide the `GTF` to Salmon to get gene-level annotations.

Here's Salmon's help info for `--geneMap`:

> File containing a mapping of transcripts to genes. If this file is provided Salmon will output both quant.sf and quant.genes.sf files, where the latter contains aggregated gene-level abundance estimates. The transcript to gene mapping should be provided as either a GTF file, or a in a simple tab-delimited format where each line contains the name of a transcript and the gene to which it belongs separated by a tab. The extension of the file is used to determine how the file should be parsed. Files ending in '.gtf', '.gff' or '.gff3' are assumed to be in GTF format; files with any other extension are assumed to be in the simple format. In GTF / GFF format, the "transcript\_id" is assumed to contain the transcript identifier and the "gene\_id" is assumed to contain the corresponding gene identifier.

Here's what I've tried:

1. NCBI seems to only provide `GFF` files. I know you can convert between the two, but that hasn't worked well for me in the past. Both file formats allow a lot of freedom, which makes conversions sloppy. Salmon, specifically, is looking for the `gene_id` and `transcript_id` annotations.
2. UCSC provides `GTF` files from RefSeq, but the `gene_id` annotation is identical to the `transcript_id` annotation (i.e., it's the NM number). Or maybe there's an option I'm missing.
3. Ensembl provides exactly the format and all of the information I could hope for, except the `transcript_id` is the Ensembl ID (ENST#), not RefSeq (NM#).

Can anyone tell me how to get a `GTF` file with hg38 RefSeq annotations?

# Answer
You will probably be interested in the following UCSC wiki page, which explains how to go from most of the UCSC tables to GTF/GFF:
http://genomewiki.ucsc.edu/index.php/Genes\_in\_gtf\_or\_gff\_format

The basic gist is that UCSC doesn't store any data internally as GTF or GFF, and so you will need to use our genePredToGtf utility in order to convert from our internally stored genePred format into GTF/GFF3.

This is the command I usually give when we get this question on our mailing list:

```
$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -N -e "select * from refGene" hg19 | cut -f2- | genePredToGtf -source=hg19.refGene.ucsc file stdin stdout

```

Change stdout to the output filename you want in the last command to get an hg19 refGene GTF file:

```
chr1    hg19.refGene.ucsc   transcript  11869   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357";  gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    11869   12227   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "1"; exon_id "NR_148357.1"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    12613   12721   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "2"; exon_id "NR_148357.2"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    13221   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "3"; exon_id "NR_148357.3"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   transcript  11874   14409   .   +   .   gene_id "DDX11L1"; transcript_id "NR_046018";  gene_name "DDX11L1";
...

```

If you have further questions about UCSC data or tools feel free to send your question to one of the below mailing lists:

* General questions: genome@soe.ucsc.edu
* Questions involving private data: genome-www@soe.ucsc.edu
* Questions involving mirror sites: genome-mirror@ose.ucsc.edu
> 5 votes

---
Tags: ngs, public-databases, sequence-annotation, human-genome, gtf
---


# Question
Title: Sleuth: transcripts with beta close to 0 are considered differentially expressed in a likelihood-ratio test
I'm comparing the results that I obtain when doing a DE analysis with the Wald test and the likelihood-ratio test. One the thing that I've noticed is that there are many genes with 'beta' close to zero that are considered differentially expressed between the conditions.

I know that the likelihood-ratio test does't use the beta values to calculate the p-values, but I find it strange that transcripts with similar expression values between the conditions are being considered differentially expressed.

**Volcano plot (Wald test q-values):**

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

# Answer
> 1 votes
This question was answered here

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