# Question
Title: Software recommendation: find DNA sequence distribution over entire transcript
I would like to create a density/histogram of the distribution of a particular DNA sequence over the entire transcript using R and/or command line tools. From here, I would like to use the coordinates of the bins to map the intron-exon diagram below the plot. In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms.

I have already found the number of occurrences of the sequences I am interested in using the `seqinr` and `biostrings` packages. I have been using `biomaRt` to get the exon start and stop sequences as well as the full gene sequence. Between all this, I feel I have all the tools I need but I am missing how to put it all together.

The issue I am having is that `biomaRt` returns multiple exon sequences with identical start positions but different end positions for a total of 37 exons sequences. Meanwhile, Genatlas and NCBI are returning anywhere from 19 to 21 exons for Grin1, my gene of interest. How can I put all of these tools together to get the plot I need? Is there a better approach?

Finally, if there is an R package or command line tool that already has this function built in, that would be ideal.

Example code for one transcript. This code is actually from an earlier project where I used multiple transcripts but it works fine with one as well. I masked the sequences I'm searching for as "XXXX" at the request of my project manager.

```
library(magrittr)
library(biomaRt)
library(org.Mm.eg.db)
library(biostrings)
library(seqinr)

### Select only the longest sequences
getMaxLengthSequences <- function(x) {
  seqMax <- split(x, factor(x$mgi_symbol))
  seqMax <- lapply(seqMax, function(x) max(x$length)) %>%
    unlist() %>%
    as.data.frame()
  seqMax$mgi_symbol <- row.names(seqMax)
  row.names(seqMax) <- NULL
  names(seqMax) <- c("length", "mgi_symbol")
  seqUnique <- merge(seqMax, x)
  seqUnique <- seqUnique[order(seqUnique$mgi_symbol), ]

  return(seqUnique)
}

### Get sequence for gene of interest
getSequences <- function(genelist, seqType) {
  ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

  sequ <- biomaRt::getSequence(id = genelist,
                               type = "mgi_symbol",
                               seqType = seqType,
                               mart = ensembl)

  ### Remove unavailable sequences
  sequ <- sequ[!grepl("Sequence unavailable", sequ[ ,seqType], fixed = TRUE), ]

  ### Find the length of each sequence then pick the max length sequence if there are multiple sequences returned from BioMart
  sequ <- sequ[order(sequ$mgi_symbol), ]
  sequ$length <- sapply(sequ[ ,seqType], nchar)
  sequ <- getMaxLengthSequences(sequ)

  return(sequ)
}

### Get motif counts of interest (substrings are masked)
getmotifs <- function(x) {

  seq <- x %>%
    tolower() %>%
    s2c()

  output <- c("XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXX" = count(seq, 3)["xxx"] %>% sum())

  return(output)
}

### Using the above functions: Get sequence and then find motifs
seq <- getSequences("Grin1")
targetMotifs <- matrix(sapply(seq[ ,"coding"], getmotifs), nrow = nrow(targets), ncol = 3, byrow = TRUE)

targetMotifTable <- data.frame(Length = seq$length)
rownames(targetMotifTable) <- make.names(seq$mgi_symbol, unique = TRUE)
targetMotifTable <- cbind(targetMotifTable, targetMotifs) %>% as.data.frame()
names(targetMotifTable) <- c("Length", "XXXX", "XXXX", "XXX")

### Use biomaRt to get exons positions
ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

gb <- getBM(attributes = c('ensembl_exon_id', "exon_chrom_start","exon_chrom_end","gene_exon"),
            filters = "mgi_symbol",
            values="Grin1",
            mart = ensembl,
            bmHeader = TRUE)

```

# Answer
If I understand the question correctly, you'd like to plot the positions of the matches to you motif along with a gene model that shows the positions of introns and exons for the different transcripts. This can be accomplished fairly easily with ggbio:

```
library(EnsDb.Mmusculus.v79)
library(ggbio)
library(biomaRt)
library(stringr)
library(dplyr)

# plot gene model with ggbio
gene_model <- autoplot(EnsDb.Mmusculus.v79, ~ symbol == "Grin1")

ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

# You need the full length gene sequence to match the gene model
seq <- biomaRt::getSequence(id = "Grin1",
                         type = "mgi_symbol",
                         seqType = "gene_exon_intron",
                         mart = ensembl)

# retrieve position of gene in genome to match up with gene model
gene_start=biomaRt::getBM(attributes=c("start_position"), filters = 'external_gene_name', values = 'Grin1', mart=ensembl)

# use stringR to find positions of exact matches
motif <- "ggcc"
motif_positions <- as.data.frame(str_locate_all(seq$gene_exon_intron, toupper(motif))) %>%
  mutate(start=start+gene_start[[1]]-1)

# make plot similar to Llopis's answer
motif_plot <- ggplot(motif_positions, aes(x=start)) +
    geom_point(y=1) +
    scale_y_continuous(limits=c(0,1))

# combine plots
tracks(motif= motifs, gene_model=gene_model)

```

This gives all motif matches for the gene sequence, not just the transcripts. It would be quite straightforward to filter out positions that overlap introns
> 5 votes

# Answer
I assume you have the sequence in seq1:

```
# Set the input
seq <- getSequences("GRIN1")
motif <- "ggcc"
# Convert to lower case
seq2 <- tolower(seqinr::s2c(as.character(seq)))
istarts <- seq(from = 1 , to = length(seq2), by = 1)
# Create words of the same size as the motif
oligos <- seq2[istarts]
for (i in 2:nchar(motif)) {
        oligos <- paste(seq2, seq2[istarts + i - 1], sep = "")
}
# Remove exceeding words (words with NA, might need tweaking to motif size
oligos <- oligos[-c(length(oligos):(length(oligos)-nchar(motif)+2))]
# Find motif in the sequence
hits <- grep(motif, oligos, ignore.case = TRUE)
#Prepare the plot
positions <- ifelse(1:length(seq2) %in% hits, 1, 0)
plot(1:length(seq2), positions, col = c("white", "black")[position+1], pch = 19)

```

The plot is not the density, but you can choose how to plot it now, as you might be interested in windows sizes of 5 nucleotides or 20.

From here you can iterate through multiple sequences. If you want to find several motifs you can take advantage of the oligos object if they are of the same size.

There is the seqPattern package in Bioconductor for more advanced pattern matching. But surprisingly it doesn't have this simple motif search. I think it would be easier to do this with perl and pyhton than in R.

---

1: In the code posted the function getSequence returns an error.
> 1 votes

---
Tags: r, software-recommendation, sequence-annotation
---


title: snakemake rule all

text:
# Question
Title: snakemake rule all
I am trying snakemake for my mRNA-seq pipeline and I got stuck in a simple thing.

```
from os.path import join

RSEM_INDEX = "/app/ref/ensembl/human/idx/GRCh37"
SAMPLES, = glob_wildcards("Align/{sample}/Aligned.toTranscriptome.out.bam")

print("---Sample Names---")

for i in SAMPLES:
    print(i)

rule all:
    input:
        expand("Expression/{sample}/rsem.genes.results", sample=SAMPLES)

rule expression:
    input:
        bam = join("Align", "{wildcards.sample}}", "Aligned.toTranscriptome.out.bam")
    params:
        rsem_index = RSEM_INDEX
    output:
        join("Expression", "{wildcards.sample}", "rsem.genes.results")
    threads:
        8
    shell:
        """
        rsem-calculate-expression --estimate-rspd --bam --no-bam-output -p {threads} --paired-end Align/{sample}/Aligned.toTranscriptome.out.bam {params.rsem_index} Expression/{wildcards.sample}/rsem.genes.results >& Expression/{wildcards.sample}/rsem.log
        """

```

Here is an error I got

```
    MissingInputException in line 12 of /app/src/mRNAseq/TBD170049_20170320_RNAseq/quantify.snakemake:
Missing input files for rule all:
Expression/013_clamp_0m_buffy/rsem.genes.results
Expression/12_27_11_clamp_buffy_200/rsem.genes.results
Expression/010_200_buffy/rsem.genes.results
Expression/010_0_buffy/rsem.genes.results
Expression/12_27_11_clamp-p_buffy_0/rsem.genes.results
Expression/013_clamp_200m_buffy/rsem.genes.results
Expression/11_15_006_0m_buffy/rsem.genes.results
Expression/11_15_006_200m_buffy/rsem.genes.results

```

So, I thought rule all input for your final files you want to make, but it's looking for the final files that I want to make.

I would like to get some help. would you help me please? Thank you in advance!

My file trees.

```
drwxrwxr-x. 10 corea corea  4096 Sep 26 20:54 Align
-rw-rw-r--.  1 corea corea  1295 Sep 27 10:19 align.snakemake
-rw-rw-r--.  1 corea corea  1823 Sep 27 09:38 align.snakemake.test
-rw-rw-r--.  1 corea corea 31055 Sep 26 08:55 out
-rw-rw-r--.  1 corea corea   869 Sep 27 12:34 quantify.snakemake
drwxrwxr-x.  2 corea corea  4096 Sep 22 20:50 Rawdata
(mRNAseq) [corea@Biocore76 TBD170049_20170320_RNAseq]$ tree Align/
Align/
├── 010_0_buffy
│   ├── Aligned.out.bam
│   ├── Aligned.toTranscriptome.out.bam
│   ├── Log.final.out
│   ├── Log.out
│   ├── Log.progress.out
│   ├── SJ.out.tab
│   ├── _STARgenome
│   │   ├── sjdbInfo.txt
│   │   └── sjdbList.out.tab
│   └── _STARpass1
│       ├── Log.final.out
│       └── SJ.out.tab
├── 010_200_buffy
│   ├── Aligned.out.bam
│   ├── Aligned.toTranscriptome.out.bam
│   ├── Log.final.out
│   ├── Log.out
│   ├── Log.progress.out
│   ├── SJ.out.tab
│   ├── _STARgenome
│   │   ├── sjdbInfo.txt
│   │   └── sjdbList.out.tab
│   └── _STARpass1
│       ├── Log.final.out
│       └── SJ.out.tab
├── 013_clamp_0m_buffy

```

# Answer
The following (or something close) should work:

```
rule expression:
    input:
        "Align/{sample}/Aligned.toTranscriptome.out.bam"
    params:
        rsem_index = RSEM_INDEX
    output:
        "Expression/{sample}/rsem.genes.results"
    threads:
        8
    shell:
        """
        rsem-calculate-expression --estimate-rspd --bam --no-bam-output -p {threads} --paired-end {input} {params.rsem_index} {output} >& Expression/{wildcards.sample}/rsem.log
        """

```

BTW, the original issue was probably due to the `}}` in your input line.
> 3 votes

---
Tags: snakemake
---


title: Why PHMM is used For Protein Analysis

text:
# Question
Title: Why PHMM is used For Protein Analysis
I am interested on how profile hidden markov model is used for protein analysis in bioinformatics. Can you guys explain why phmm is more preferable method comparing to simple MSA(Multiple Sequence Alignment) scoring or simple MSA probability scoring? For example we have three MSA sequences as based model and a target to be scored.

### MSA Sequences

A B C B C A
A B B C C A
A B C B A A

### Target Sequence

A B B B C A

By MSA scoring we can count score for target sequence comparing to a profile MSA. For example above target sequence protein column will counted as match(A), match(B), mismatch(B-majorityMSA(C)), match(B-majorityMSA(B)), match(C-majorityMSA(C)), match(A). Just give score for match and mismatch and we can get total score for a target sequence comparing to a profile MSA sequences.

By MSA probability scoring we can count probability of each protein in each column so for example above
**P(target\_seq similar to MSA profile)** = P(A in MSA sequences)\*P(B in MSA sequences)\*.... \*P(A in MSA sequences) = 1\*1\*(1/3)\*(2/3)\*(2/3)\*1 = 4/27

Isn't MSA scoring like above examples powerful enough?

Thanks a lot !!

# Answer
> 1 votes
The simple probabilistic scoring scheme you describe to query a target sequence against a (protein) family represented as a multiple sequence alignment (MSA) is actually similar in essence to a Position Weight Matrix (PWM), although PWMs are more refined. PWMs are powerful rapid tools often used for motif analyses.

However, your approach (and PWMs) don't handle insertions or deletions (indels), and you can't easily add these into your model (for example, consider that 2 1bp indels should not score the same as 1 2bp indel, and indels could occur within your MSA as well as in your target sequence). You could try and add some arbitrary gap opening and extension penalties somehow (like BLAST) but then you'd abandon your probabilistic model! By contrast, pHMMs have a fully fledged probabilistic indel model with insertion and deletion states.

Finally, of practical note, there are some excellent implementations of pHMMs for protein sequence analysis, such as HMMER3 and HHSuite.

# Answer
> 0 votes
I recommend you to read this paper in which the applications of HMM in protein sequence alignments are discussed.

---
Tags: sequence-analysis, sequence-alignment, hidden-markov-models
---


# Question
Title: Installing Sibelia for Ragout on Mac OSX
I am trying to use Ragout: https://github.com/fenderglass/Ragout

to fill the gaps in my *de novo* genome assembly. You can access the article freely here: https://www.ncbi.nlm.nih.gov/pubmed/24931998

For this, I first need to install Sibelia.

I tried doing so by cloning the GitHub repository into output/software/ragout/ and then running the following command:

```
git clone https://github.com/fenderglass/Ragout output/software/ragout

python2 output/software/ragout/scripts/install-sibelia.py

Installing Sibelia
Downloading source...
-- The CXX compiler identification is AppleClang 9.0.0.9000037
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- The C compiler identification is AppleClang 9.0.0.9000037
-- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/cc
-- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Performing Test HAVE_UNKNOWN_WALL
-- Performing Test HAVE_UNKNOWN_WALL - Success
-- Performing Test HAVE_UNKNOWN_FOMIT_FRAME_POINTER
-- Performing Test HAVE_UNKNOWN_FOMIT_FRAME_POINTER - Success
-- Looking for inttypes.h
-- Looking for inttypes.h - found
-- Looking for memory.h
-- Looking for memory.h - found
-- Looking for stddef.h
-- Looking for stddef.h - found
-- Looking for stdint.h
-- Looking for stdint.h - found
-- Looking for stdlib.h
-- Looking for stdlib.h - found
-- Looking for string.h
-- Looking for string.h - found
-- Looking for strings.h
-- Looking for strings.h - found
-- Looking for sys/types.h
-- Looking for sys/types.h - found
-- Performing Test HAVE_INLINE
-- Performing Test HAVE_INLINE - Success
-- Performing Test HAVE___INLINE
-- Performing Test HAVE___INLINE - Success
-- Performing Test HAVE___INLINE__
-- Performing Test HAVE___INLINE__ - Success
-- Performing Test HAVE___DECLSPEC_DLLEXPORT_
-- Performing Test HAVE___DECLSPEC_DLLEXPORT_ - Failed
-- Performing Test HAVE___DECLSPEC_DLLIMPORT_
-- Performing Test HAVE___DECLSPEC_DLLIMPORT_ - Failed
-- Check size of uint8_t
-- Check size of uint8_t - done
-- Check size of int32_t
-- Check size of int32_t - done
-- Looking for PRId32
-- Looking for PRId32 - found
-- Configuring done
-- Generating done
-- Build files have been written to: /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/build
Scanning dependencies of target divsufsort
[  4%] Building C object libdivsufsort-2.0.1/lib/CMakeFiles/divsufsort.dir/divsufsort.o
[  8%] Building C object libdivsufsort-2.0.1/lib/CMakeFiles/divsufsort.dir/sssort.o
[ 12%] Building C object libdivsufsort-2.0.1/lib/CMakeFiles/divsufsort.dir/trsort.o
[ 16%] Building C object libdivsufsort-2.0.1/lib/CMakeFiles/divsufsort.dir/utils.o
[ 20%] Linking C static library libdivsufsort.a
[ 20%] Built target divsufsort
Scanning dependencies of target Sibelia
[ 24%] Building CXX object CMakeFiles/Sibelia.dir/sibelia.cpp.o
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/sibelia.cpp:8:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:43:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence.h:103:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence/string_set_dependent_generous.h:38:9: warning:
      'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_' is used as a header
      guard here, followed by #define of a different macro [-Wheader-guard]
#ifndef SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_
        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence/string_set_dependent_generous.h:39:9: note:
      'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUSH_' is defined here; did you
      mean 'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_'?
#define SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUSH_
        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/sibelia.cpp:8:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:44:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score.h:48:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score/score_matrix.h:40:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/file.h:70:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system.h:82:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system/system_sema.h:120:27: warning:
      'sem_init' is deprecated [-Wdeprecated-declarations]
            SEQAN_DO_SYS(!sem_init(hSemaphore, 0, init));
                          ^
/usr/include/sys/semaphore.h:55:42: note: 'sem_init' has been explicitly marked
      deprecated here
int sem_init(sem_t *, int, unsigned int) __deprecated;
                                         ^
/usr/include/sys/cdefs.h:176:37: note: expanded from macro '__deprecated'
#define __deprecated    __attribute__((deprecated))
                                       ^
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/sibelia.cpp:8:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:44:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score.h:48:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score/score_matrix.h:40:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/file.h:70:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system.h:82:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system/system_sema.h:124:27: warning:
      'sem_destroy' is deprecated [-Wdeprecated-declarations]
            SEQAN_DO_SYS(!sem_destroy(hSemaphore));
                          ^
/usr/include/sys/semaphore.h:53:26: note: 'sem_destroy' has been explicitly
      marked deprecated here
int sem_destroy(sem_t *) __deprecated;
                         ^
/usr/include/sys/cdefs.h:176:37: note: expanded from macro '__deprecated'
#define __deprecated    __attribute__((deprecated))
                                       ^
3 warnings generated.
[ 28%] Building CXX object CMakeFiles/Sibelia.dir/postprocessor.cpp.o
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.cpp:9:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:43:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence.h:103:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence/string_set_dependent_generous.h:38:9: warning:
      'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_' is used as a header
      guard here, followed by #define of a different macro [-Wheader-guard]
#ifndef SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_
        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/sequence/string_set_dependent_generous.h:39:9: note:
      'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUSH_' is defined here; did you
      mean 'SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_'?
#define SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUSH_
        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        SEQAN_SEQUENCE_STRING_SET_DEPENDENT_GENEROUS_H_
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.cpp:9:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:44:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score.h:48:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score/score_matrix.h:40:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/file.h:70:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system.h:82:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system/system_sema.h:120:27: warning:
      'sem_init' is deprecated [-Wdeprecated-declarations]
            SEQAN_DO_SYS(!sem_init(hSemaphore, 0, init));
                          ^
/usr/include/sys/semaphore.h:55:42: note: 'sem_init' has been explicitly marked
      deprecated here
int sem_init(sem_t *, int, unsigned int) __deprecated;
                                         ^
/usr/include/sys/cdefs.h:176:37: note: expanded from macro '__deprecated'
#define __deprecated    __attribute__((deprecated))
                                       ^
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.cpp:9:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/postprocessor.h:10:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/align.h:44:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score.h:48:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/score/score_matrix.h:40:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/file.h:70:
In file included from /Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system.h:82:
/Users/cr517/Documents/phd/project/sibelia-build/Sibelia-master/src/include/seqan/system/system_sema.h:124:27: warning:
      'sem_destroy' is deprecated [-Wdeprecated-declarations]
            SEQAN_DO_SYS(!sem_destroy(hSemaphore));
                          ^
/usr/include/sys/semaphore.h:53:26: note: 'sem_destroy' has been explicitly
      marked deprecated here
int sem_destroy(sem_t *) __deprecated;
                         ^
/usr/include/sys/cdefs.h:176:37: note: expanded from macro '__deprecated'
#define __deprecated    __attribute__((deprecated))
                                       ^
3 warnings generated.
[ 32%] Building CXX object CMakeFiles/Sibelia.dir/indexedsequence.cpp.o
[ 36%] Building CXX object CMakeFiles/Sibelia.dir/util.cpp.o
[ 40%] Building CXX object CMakeFiles/Sibelia.dir/outputgenerator.cpp.o
[ 44%] Building CXX object CMakeFiles/Sibelia.dir/blockfinder.cpp.o
[ 48%] Building CXX object CMakeFiles/Sibelia.dir/blockinstance.cpp.o
[ 52%] Building CXX object CMakeFiles/Sibelia.dir/bifurcationstorage.cpp.o
[ 56%] Building CXX object CMakeFiles/Sibelia.dir/bulgeremoval.cpp.o
[ 60%] Building CXX object CMakeFiles/Sibelia.dir/dnasequence.cpp.o
[ 64%] Building CXX object CMakeFiles/Sibelia.dir/edge.cpp.o
[ 68%] Building CXX object CMakeFiles/Sibelia.dir/fasta.cpp.o
[ 72%] Building CXX object CMakeFiles/Sibelia.dir/serialization.cpp.o
[ 76%] Building CXX object CMakeFiles/Sibelia.dir/synteny.cpp.o
[ 80%] Building CXX object CMakeFiles/Sibelia.dir/test/unrolledlisttest.cpp.o
[ 84%] Building CXX object CMakeFiles/Sibelia.dir/platform.cpp.o
[ 88%] Building CXX object CMakeFiles/Sibelia.dir/stranditerator.cpp.o
[ 92%] Building CXX object CMakeFiles/Sibelia.dir/vertexenumeration.cpp.o
[ 96%] Building CXX object CMakeFiles/Sibelia.dir/resource.cpp.o
[100%] Linking CXX executable Sibelia
[100%] Built target Sibelia
[ 20%] Built target divsufsort
[100%] Built target Sibelia
Install the project...
-- Install configuration: "Release"
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/bin/Sibelia
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/NEWS.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/ANNOTATION.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/README.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/USAGE.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/INSTALL.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/SIBELIA.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/C-SIBELIA.md
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/LICENSE.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia/Staphylococcus_aureus
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia/Staphylococcus_aureus/NCTC8325.fasta
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia/Staphylococcus_aureus/README.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia/Staphylococcus_aureus/RN4220.fasta
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/C-Sibelia/Staphylococcus_aureus/variant.vcf
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/blocks_coords.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.conf
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.highlight.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.highlight1.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.highlight2.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.highlight3.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.highlight4.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.image.conf
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.png
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.segdup.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.sequences.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/circos/circos.svg
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/coverage_report.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/d3_blocks_diagram.html
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/genomes_permutations.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/Helicobacter_pylori.fasta
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Helicobacter_pylori/README.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/blocks_coords.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.conf
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.highlight.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.image.conf
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.png
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.segdup.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.sequences.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/circos/circos.svg
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/coverage_report.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/d3_blocks_diagram.html
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/genomes_permutations.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/README.txt
-- Installing: /Users/cr517/Documents/phd/project/sibelia-build/share/Sibelia/doc/examples/Sibelia/Staphylococcus_aureus/Staphylococcus.fasta
Traceback (most recent call last):
  File "output/software/ragout/scripts/install-sibelia.py", line 111, in <module>
    sys.exit(main())
  File "output/software/ragout/scripts/install-sibelia.py", line 107, in main
    return int(not install_deps(args.prefix))
  File "output/software/ragout/scripts/install-sibelia.py", line 31, in install_deps
    return install_sibelia(prefix)
  File "output/software/ragout/scripts/install-sibelia.py", line 64, in install_sibelia
    shutil.copy(sibelia_bin_src, sibelia_bin_dst)
  File "/usr/local/Cellar/python/2.7.13/Frameworks/Python.framework/Versions/2.7/lib/python2.7/shutil.py", line 119, in copy
    copyfile(src, dst)
  File "/usr/local/Cellar/python/2.7.13/Frameworks/Python.framework/Versions/2.7/lib/python2.7/shutil.py", line 83, in copyfile
    with open(dst, 'wb') as fdst:
IOError: [Errno 2] No such file or directory: '/Users/cr517/Documents/phd/project/lib/Sibelia'

```

As you can see, I get a lot of errors. Could someone help me, please?

I also tried to install Sibelia following their instructions here: https://github.com/bioinf/Sibelia

```
cd Sibelia/build
➜  build git:(master) cmake ../src
-- The CXX compiler identification is AppleClang 9.0.0.9000037
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++
-- Check for working CXX compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- The C compiler identification is AppleClang 9.0.0.9000037
-- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/cc
-- Check for working C compiler: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Performing Test HAVE_UNKNOWN_WALL
-- Performing Test HAVE_UNKNOWN_WALL - Success
-- Performing Test HAVE_UNKNOWN_FOMIT_FRAME_POINTER
-- Performing Test HAVE_UNKNOWN_FOMIT_FRAME_POINTER - Success
-- Looking for inttypes.h
-- Looking for inttypes.h - found
-- Looking for memory.h
-- Looking for memory.h - found
-- Looking for stddef.h
-- Looking for stddef.h - found
-- Looking for stdint.h
-- Looking for stdint.h - found
-- Looking for stdlib.h
-- Looking for stdlib.h - found
-- Looking for string.h
-- Looking for string.h - found
-- Looking for strings.h
-- Looking for strings.h - found
-- Looking for sys/types.h
-- Looking for sys/types.h - found
-- Performing Test HAVE_INLINE
-- Performing Test HAVE_INLINE - Success
-- Performing Test HAVE___INLINE
-- Performing Test HAVE___INLINE - Success
-- Performing Test HAVE___INLINE__
-- Performing Test HAVE___INLINE__ - Success
-- Performing Test HAVE___DECLSPEC_DLLEXPORT_
-- Performing Test HAVE___DECLSPEC_DLLEXPORT_ - Failed
-- Performing Test HAVE___DECLSPEC_DLLIMPORT_
-- Performing Test HAVE___DECLSPEC_DLLIMPORT_ - Failed
-- Check size of uint8_t
-- Check size of uint8_t - done
-- Check size of int32_t
-- Check size of int32_t - done
-- Looking for PRId32
-- Looking for PRId32 - found
-- Configuring done
-- Generating done
-- Build files have been written to: /Users/cr517/Documents/phd/project/Sibelia/build
➜  build git:(master) make
Scanning dependencies of target lagan
[  3%] Creating directories for 'lagan'
[  6%] Performing download step (DIR copy) for 'lagan'
[  9%] No patch step for 'lagan'
[ 12%] No update step for 'lagan'
[ 15%] No configure step for 'lagan'
[ 18%] Performing build step for 'lagan'
rightinfluence.cpp:21:2: error: reference to 'end' is ambiguous
        end.score = -2;
        ^
rightinfluence.cpp:3:18: note: candidate found by name lookup is 'end'
Fragment origin, end;
                 ^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/iterator:1768:1: note:
      candidate found by name lookup is 'std::__1::end'
end(const _Cp& __c)
^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/iterator:1611:1: note:
      candidate found by name lookup is 'std::__1::end'
end(_Tp (&__array)[_Np])
^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/iterator:1760:1: note:
      candidate found by name lookup is 'std::__1::end'
end(_Cp& __c)
^
rightinfluence.cpp:22:22: error: reference to 'end' is ambiguous
        origin.totalScore = end.totalScore = 0;
                            ^
rightinfluence.cpp:3:18: note: candidate found by name lookup is 'end'
Fragment origin, end;
                 ^
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/include/c++/v1/iterator:1768:1: note:

```

...

# Answer
There is a link to compiled binaries for OSX on the project homepage
> 4 votes

---
Tags: genomics, scaffold
---


title: *very* unbalanced group sizes for DE

text:
# Question
Title: *very* unbalanced group sizes for DE
I downloaded some publicly available RNA-seq data and want to compare those samples carrying a mutation (~4) against the rest (~800!). I ran both EdgeR and DESeq2, and the first results in an asymmetric volcano plot: skewed in one of the sides, meaning that it results in very large logFC that have large (NS) adjusted p-values. So next, in order to see if this asymmetry is persistant, I ran several random partitions of the data as groups (but keeping them highly unbalanced), obtaining the same behavior. Here an example of what I mean (I run it with 4, 20 and 100 to show how this effect decreases as I increase the small class size)

\*I tried EdgeR-robust and I still get the same effect.

I do not see this effect when running DESeq2, only very few DE genes. What method would you recommend using when comparing such different class sizes?

# Answer
Random partitioning (as you have done) seems like a reasonable thing to do to work out what's going on. From the results you have suggested, it looks like DESeq2 is performing better than EdgeR, so use DESeq2. Can you show the volcano plots and the partition sizes?

Note: I'm a little bit biased in this, because DESeq2 is what I prefer to use (particularly where Salmon is not appropriate).
> 3 votes

# Answer
In general, you would expect roughly the same distribution on both sides of your volcano plot. The first two plots are concerning since there's so much on the right and not on the left. With such a small number of samples, it is hard to conclude if this is an artifact of not having enough samples, a technical artifact from the way the data was processed, or true biology. I remember reading that DESeq and other packages like a high number of samples for estimating dispersion or something, so this could also be contributing. I think you should find a better dataset (more samples) to prove or disprove your hypothesis.
> 3 votes

---
Tags: rna-seq, deseq2, edger
---


# Question
Title: How to apply upperquartile normalization on RSEM expected counts?
I see that TCGA RNASeq V2 RSEM data is normalized with upper-quartile normalization.

After doing Quantification with RSEM with the samples I have, I got "genes.results" as output which has gene id, transcript id(s), length, expected count, and FPKM. So, from all the sample output files I got the gene\_id and expected\_counts \[from all samples\]. for eg: it looks like following:

```
Gene_id S1  S2  S3  S4  S5
ENSG00000000003 1355    1121    242 709 1071
ENSG00000000005 5   0   0   0   0
ENSG00000000419 839 1345    48  976 1234
ENSG00000000457 429 1803.08 386 1156    628
ENSG00000000460 404 1444    523 479 1260
ENSG00000000938 294 312 93  189 683
ENSG00000000971 3911    4633.92 264 2863    5194
ENSG00000001036 1546    2276    271 1641    2141
ENSG00000001084 1945    2125    490 980 1533
ENSG00000001167 2054    4069.72 3457    2075    2207
ENSG00000001460 153 339 77  221 273

```

I want to apply "upper-quartile" normalization on this data. But doesn't know which R package I can use and which function?

Is anyone aware about this?

P.S. This is not for differential analysis

# Answer
You can use the `quantile` function in base R to get the value of a particular quantile (e.g. 0.75 for the upper quartile). This can then be used as a factor for normalisation: divide the observed expression by this number.

```
## set up simulated data
data.mat <- floor(10^data.frame(sapply(1:5,function(x){rnorm(100)})));
colnames(data.mat) <- paste0("S",1:5);
data.mat[data.mat > 0] <- floor(data.mat[data.mat > 0]);

```

It's generally a good idea to first remove genes that have no counts in any samples. Rows can be removed using a logical selection based on the total expression in any given row.

```
## Only retain rows that have a positive total expression per gene
data.mat <- data.mat[rowSums(data.mat) > 0,];

```

Two methods of upper quartile determination are presented here, one which is based on the total counts, and one based on the expressed counts (i.e. excluding values with no expression). These both use the `quantile` function of R. The `apply` function calculates statistics that are based on a single specified dimension (or multiple dimensions) of the data. The first dimension of a 2D matrix in R corresponds to the row, while the second corresponds to the column. By using `2` in the apply function, data are collected along the second dimension (i.e. the column), and statistics generated from that:

```
data.quantileAll <- apply(data.mat, 2, function(x){quantile(x, 0.75)});
data.quantileExpressed <-
  apply(data.mat, 2, function(x){quantile(x[x>0], 0.75)});

```

Expressed counts may be preferred when read counts are low, so it's difficult to distinguish between a transcript with very low expression and a transcript with no expression.

Normalisation is applied by dividing the values by the column-based upper quartile statistic. R will repeat values along the last dimension first (which in this case is the column), so the matrix needs to be transposed first (with columns as genes) to get the arithmetic working properly, with a final transposition to return the orientation back to the original direction (with columns as samples).

```
data.norm <- t(t(data.mat) / data.quantileExpressed);

```
> 9 votes

---
Tags: r, rna-seq, normalization, rsem
---


# Question
Title: What are the algorithms for biometric identification in large database?
I have been thinking about this problem and reading some related papers for sometime.

The problem can be described as follows:

Suppose that each biometric template is a $d-$dimensional vector. A database is composed of $n$ such templates ( $n$ could be large, like 1 million). Given a query template $q$, the task is to identify the "closest" template according to some distance metric.

Now, since the dimension $d$ and the database size $n$ could be both large, it is not quite efficient to linearly scan the whole database.

So my question is, are there any algorithms to efficiently fulfill the above task? What are the state of the art algorithms for biometric identification?

Through some literature review, I can think of two approaches. But both of them seem to have some limitations.

1. Building space partitioning index like R-tree. The problem is that such index does not suit for high dimensional data. The final search efficiency would be close to linear search.
2. Locality sensitive hash. The problem is that this method can only provide **approximate** searching result. For biometic identification, it is not guaranteed that the matched template could be found.

# Answer
I think that a very nice example of what you are talking about saying the word bigger, is probably a microarray transciptional analysis of RNA extracted from a biological sample. Simplifying, in these analyis, each observation is analised on, up to 10000 genes expression levels of different mRNA coding genes. Thus, from these type of molecular biology experiments, we obtain matrices with few observations (d, the samples) respect to the dimensions observed (k, the target gene expression levels). Well, d \<\< k. In this situation, in order to perform a classification-like task, and to reduce the problem of the overfitting related with dimensionality maledition, the ways that we might follow are several:

1) PCA: it is a good tool to visualize multidimensional datasets reducing dimensionality and noise level in the consolidated and standardized data. Problem: it is an unsupervised machine learning algorithm. Alone, it haven't the capacity of class discrimination.

2) PCA-LDA is a good combination, where PCA eliminates hyperdimensionality and background in our data sheets, and than LDA (linear discriminant analysis) finds the directions onto the new subspace of characteristics, that are able to maximize variances between classes and decrease the variances within these. For this reason LDA is known as a kind of supervised machine learning technique. Furthermore, PCA and LDA are usually used together because LDA applied on high dimensional datasets could lead us toward overfitting problems. After the training of PCA-LDA pipeline, being LDA a classificator object trained and shaped on the principal component subspace, you can use this to examine new unknown samples or observation and predict its class by similarity. e.g. recently, Martin et al., on a Nature article talk about PCA-LDA approach to identify and than to predict cancer cells from those are normal, starting from ranman spectra of tissues.

3) PCA-LDA on transpose matrix: to use as better as possible the huge dimensionality, some authors used to transpose the initial matrix in order to force positioning of genes as rows and samples as columns. Then, the subsequent LDA is performed not on PC scores, but on loading factors of the eigenfactor matrix obtaining thus the list of the most important discriminant gene expression levels to segregate normal samples to cancer ones.
> 2 votes

---
Tags: database, identifiers
---


# Question
Title: How to check whether all BAM read contain defined read groups?
I'm trying to investigate whether there are errors within my BAM. After looking at the BAM header to see whether the read groups exists (using samtools, i.e.

```
samtools view -H file1.bam

```

Here, it appears that the header includes `@RG` tags. However, I suspect that the BAM contains reads that either lack the `@RG` tag or perhaps these read groups are not well-defined, i.e. it looks like some reads in the BAM may be missing an RG tag.

How do I explicitly check whether all reads have an RG tag, and whether these are well-defined?

# Answer
> 6 votes
According to the man page, running `samtools stats --split RG <file1.bam>` should produce summary statistics separated by read group. If it doesn't produce a list of ungrouped reads, counts/statistics can be compared to running without the `--split RG` argument.

Here's some example output from a BAM file with combined read groups:

```
$ samtools stats --split RG both.bam | grep '^SN'
SN  raw total sequences:    2392
SN  filtered sequences: 0
SN  sequences:  2392
SN  is sorted:  1
SN  1st fragments:  2392
SN  last fragments: 0
SN  reads mapped:   2341
SN  reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
SN  reads unmapped: 51
SN  reads properly paired:  0   # proper-pair bit set
SN  reads paired:   0   # paired-end technology bit set
SN  reads duplicated:   0   # PCR or optical duplicate bit set
SN  reads MQ0:  0   # mapped and MQ=0
SN  reads QC failed:    0
SN  non-primary alignments: 0
SN  total length:   2292231 # ignores clipping
SN  bases mapped:   2255361 # ignores clipping
SN  bases mapped (cigar):   1538560 # more accurate
SN  bases trimmed:  0
SN  bases duplicated:   0
SN  mismatches: 254720  # from NM fields
SN  error rate: 1.655574e-01    # mismatches / bases mapped (cigar)
SN  average length: 958
SN  maximum length: 9078
SN  average quality:    20.0
SN  insert size average:    0.0
SN  insert size standard deviation: 0.0
SN  inward oriented pairs:  0
SN  outward oriented pairs: 0
SN  pairs with other orientation:   0
SN  pairs on different chromosomes: 0

```

The `--split` command has also produced additional `.bamstat` files, which can be looked at for the summary statistics for individual read groups. In this case I'm only interested in the "raw total sequences" line:

```
$ grep 'raw total sequences' *.bamstat
both.bam_minusOnly.bamstat:SN   raw total sequences:    1901
both.bam_plusOnly.bamstat:SN    raw total sequences:    491

```

There are 1901 sequences from the "minusOnly" group, and 491 sequences from the "plusOnly" group, and I know from the previous statistics that there are 2392 total sequences. 1901 + 491 = 2392, so I know there are no ungrouped reads in this file.

# Answer
> 4 votes
using samjdk and invoking the function getReadGroup()

> getReadGroup() returns The SAMReadGroupRecord from the SAMFileHeader for this SAMRecord, or null if 1) this record has no RG tag, or 2) the header doesn't contain the read group with the given ID.or 3) this record has no SAMFileHeader

```
java -jar dist/samjdk.jar -e 'return record.getReadGroup()==null;'  input.bam

```

the output will be a sam/bam file with missing/bad RG .

# Answer
> 2 votes
A complete solution to validate a SAM/BAM file would be picard's ValidateSamFile. This will find many other issues with BAM files.

http://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile

---
Tags: sam, samtools, reads
---


# Question
Title: Uniprot: Determining protein and organism given a DNA sequence
I would like to find the identity of a protein, together with the organism name and kingdom.

The given DNA sequence is below:

```
ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTATGATATTCTTC
TTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATTCCAGATGAGTGTACCGGAGT
TAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTTAGTATGTGGTCTACTTTTGATTACGGGTTC
ATCTTGGATGGAGGAAGGTTCAGGAACGGCCTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCAT
AGCGGACTTGCTGTAGATACGTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCA
TCAACCTTATATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC
CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGCTATCACGATG
CTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGGCGGTGATCCTATGTTGTATC
AGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATATCATCATACTTCCAGTATTCGGTATAGTCAG

```

I've been asked to translate this DNA sequence into the 6 possible frames which are:

```
5'3' Frame 1
TVVR-SSVVQRPNNITCNAYDILLHHARHDGRTRQLTSAIPDECTGVSIPKD--HRYMIF
SMWSTFDYGFILDGGRFRNGLNRLSTTSAHCKS-RTCCRYVHYRIAHGRCKLPYRKHQPY
MYNRLCPPFTHGDAAVITLSLIHYNHCSVTHRSCACASRCYHDATH--KLEYQLL-QFGR
R-SYVVSALILGVWASRSLYHHTSSIRYSQ
5'3' Frame 2
LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI-FL
VCGLLLITGSSWMEEGSGTA-TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
CTIAYARRSLMAMLQSSLYP-SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
GDPMLYQHLFWVFGHPEVYIIILPVFGIV
5'3' Frame 3
CCSVIISCTTS-QHHMQCL-YSSSSCQARWQD-ATY-CHSR-VYRS-YSQRLITSVYDF-
YVVYF-LRVHLGWRKVQERPEPSIHH-RSLQVIADLL-IRSLSHCTWPVQAPLQEASTLY
VQSPMPAVHSWRCCSHHFIPDPLQSLQRYS-ELCLC-QVLSRCYSLIEVGVPASMTVRQA
VILCCISTYSGCLGIQKSISSYFQYSV-S
3'5' Frame 1
LTIPNTGSMMI-TSGCPNTQNKC-YNIGSPPAELS-KLVLQLLSVSSIVIAPASTGTTPM
SNAAVIVMDQG-SDDCSIAMSERRA-AIVHIRLMLPVRELAPAMCNAIMNVSTASPL-LA
VSASGG-TVQAVPEPSSIQDEPVIKSRPHTKNHIPMLLIFGNTNSGTLIWNGTSKLPSPA
IVPGMMKKNIISIACDVVRTLYN--SPNNS
3'5' Frame 2
-LYRILEV--YRLLDAQTPRISADTT-DHRLPNCHRSWYSNFYQ-VAS--HLLAQAQLL-
VTLQ-L-WIRDKVMTAASP-VNGGHRRLYI-G-CFL-GSLHRPCAMR--TYLQQVRYDLQ
-ALVVDRRFRPFLNLPPSKMNP-SKVDHILKIIYRCY-SLGILTPVHSSGMALVSCLVLP
SCLA--RRIS-ALHVMLLGRCTTDDHRTT
3'5' Frame 3
DYTEYWKYDDIDFWMPKHPE-VLIQHRITACRTVIEAGTPTSISE-HRDSTC-HRHNSYE
-RCSDCNGSGIK--LQHRHE-TAGIGDCTYKVDASCKGACTGHVQCDNERIYSKSAMTCS
ER-WWIDGSGRS-TFLHPR-TRNQK-TTY-KSYTDVINLWEY-LRYTHLEWH--VA-SCH
RAWHDEEEYHKHCM-CC-DVVQLMITEQQ

```

I ran jobs on Uniprot for all these sequences. Here are the results:

```
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082E0FZ
www.uniprot.org/blast/uniprot/B2017100383C3DD8CE55183C76102DC5D3A26728BEBF5BAE
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082E5FP
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082F08D
nothing
www.uniprot.org/blast/uniprot/B20171003F725F458AC8690F874DD868E4ED79B88CE32C4B

```

But how can I decide the identity of the protein and the name of the organism?

In the first one, the organism is *Pacific ridley turtle* but in the second one it is Diplonema sp. ATCC 50224. So I can only decide at a kingdom level, but not to the organism, right?

# Answer
Well, let's see. The first step is to simply run a nucleotide blast (blastn) against the nr (non-redundant) database. I did this using the NCBI blast tool. While, as expected, I got many hits, there was only one hit which matched the *entire length* of the query sequence perfectly:

```
Diplonema sp. ATCC 50224 cytochrome c oxidase subunit 1 (cox1) mRNA, complete cds; mitochondrial
Sequence ID: JF698657.1
Length: 1612Number of Matches: 1
Range 1: 141 to 770
Alignment statistics for match #1
Score   Expect  Identities  Gaps    Strand
1164 bits(630)  0.0 630/630(100%)   0/630(0%)   Plus/Plus
Query  1    ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTAT  60
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  141  ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTAT  200

Query  61   GATATTCTTCTTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATT  120
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  201  GATATTCTTCTTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATT  260

Query  121  CCAGATGAGTGTACCGGAGTTAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTT  180
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  261  CCAGATGAGTGTACCGGAGTTAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTT  320

Query  181  AGTATGTGGTCTACTTTTGATTACGGGTTCATCTTGGATGGAGGAAGGTTCAGGAACGGC  240
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  321  AGTATGTGGTCTACTTTTGATTACGGGTTCATCTTGGATGGAGGAAGGTTCAGGAACGGC  380

Query  241  CTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCATAGCGGACTTGCTGTAGATAC  300
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  381  CTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCATAGCGGACTTGCTGTAGATAC  440

Query  301  GTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCATCAACCTTAT  360
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  441  GTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCATCAACCTTAT  500

Query  361  ATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC  420
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  501  ATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC  560

Query  421  CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGC  480
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  561  CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGC  620

Query  481  TATCACGATGCTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGG  540
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  621  TATCACGATGCTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGG  680

Query  541  CGGTGATCCTATGTTGTATCAGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATAT  600
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  681  CGGTGATCCTATGTTGTATCAGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATAT  740

Query  601  CATCATACTTCCAGTATTCGGTATAGTCAG  630
            ||||||||||||||||||||||||||||||
Sbjct  741  CATCATACTTCCAGTATTCGGTATAGTCAG  770

```

So, we now know the sequence is from the gene for mitochondrial Cytochrome C Oxidase subunit 1 from *Diplonema sp.* ATCC 50224. That the sequence is mitochondrial, is very important since the mitochondrion of this species has a slightly different genetic code. If you go to the NCBI's page for this sequence (JF698657.1), you will see this line:

```
/transl_table=4

```

That tells you it is using translation table 4:

```
4. The Mold, Protozoan, and Coelenterate Mitochondrial Code and the Mycoplasma/Spiroplasma Code (transl_table=4)

    AAs  = FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
  Starts = --MM------**-------M------------MMMM---------------M------------
  Base1  = TTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGG
  Base2  = TTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGG
  Base3  = TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG

```

See how all of your translated reading frames have in-frame stop codons? That's because you translated them using the standard genetic code instead of the correct mitochondrial one. As you can see above, in this code, the codon TGA is not a STOP, as in the standard code, and instead codes for tyrosine (W).

So, if you look at the result of the protein blast you ran on UniProt, you will see:

```
Query1   LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI-FL   59
         LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI FL
F5CI61   LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGIWFL   107

Query60  VCGLLLITGSSWMEEGSGTA-TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI   118
         VCGLLLITGSSWMEEGSGTA TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
F5CI61   VCGLLLITGSSWMEEGSGTAWTVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI   167

Query119 CTIAYARRSLMAMLQSSLYP-SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG   177
         CTIAYARRSLMAMLQSSLYP SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
F5CI61   CTIAYARRSLMAMLQSSLYPWSITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG   227

Query178 GDPMLYQHLFWVFGHPEVYIIILPVFGIV                                  206
         GDPMLYQHLFWVFGHPEVYIIILPVFGIV
F5CI61   GDPMLYQHLFWVFGHPEVYIIILPVFGIV                                  256

```

The only locations where there isn't a 100% identical match is where a W is aligned to a STOP (`-` in the Uniprot output). So, your actual protein sequence is the 2nd frame of the 5' to 3' direction which, when translated using the correct code, is:

```
LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGIWFL
VCGLLLITGSSWMEEGSGTAWTVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
CTIAYARRSLMAMLQSSLYPWSITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
GDPMLYQHLFWVFGHPEVYIIILPVFGIV

```
> 4 votes

# Answer
Well if it is not a homework assignment, I would suggest to use blastx. Then you don't need to translate your DNA into protein first, but that will be covered by blastx.

If it is a homework assignment though, the teacher probably wants you to select the best blast hit. So which hit gives the highest score, and most matching (identical) amino acids between the alignment.
> 1 votes

---
Tags: database, sequence-alignment
---


title: Waterman-Smith-Beyer Sequence Alignment

text:
# Question
Title: Waterman-Smith-Beyer Sequence Alignment
I am having trouble understanding the affine gap penalty in the following example -

I am not sure where the 3 and 4 come from or the 4 and 5 in cells 1,3 and 2,4. I'm not understanding how the affine gap penalty works here. I understand how the linear gap penalty worked in the Needlemen\_Wunsch but not this.

# Answer
**The linear gap penalty** has a cost which is directly proportional to the length of the gap: the gap costs 1 unit per consecutive missing base.
The simplest case is assigning each unit (= each missing base) a penalty of 1:

> g(k) = k

Or more generally, if you want to increase the penalty per base to 2, 3 or n:

> g(k) = n * k
> with k = length of the gap, and n = cost of 1 missing base

**The affine gap penalty** has a cost that follows -as its name conveniently states- an affine function.
The idea behind it, is that any gap in an alignment has a cost with two components:

* **the "gap opening" penalty**: it is a cost payed only once per gap, basically just for existing.
* **the "gap extension" penalty:** equivalent to the linear gap penalty, it costs 1 unit for each consecutive missing base.

The general formula for the cost is then:

> g(k) = a + n\*k
> with a= opening penalty, n= extension penalty, k= length of the gap

Compared to the linear gap penalty, this will give a small bonus to longer consecutive gaps compared to skipping one nucleotide, matching one, skipping one again etc (because in the case of longer gaps, you only pay the opening penalty once, so the global cost per base decreases).

**Cells on your diagram:**
From what I understand, you're working with an affine penalty of g(k) = 1 + 1\*k. For the \[3,4\] you're "coming" from the score of 2 on the left (cell 1,2). You either:

* match the A (row 1) with G (column 3) and get a penalty of 1 for substitution -\> total of previous 2 + new 1 = 3
* open a gap to skip the G (column 3) and get a penalty of 2 (1 for opening + 1 because gap length of 1) -\> total of previous 2 + new 2 = 4.

So you choose the lesser penalty and go with 3, which is the number selected in the bigger box.
The explanation is the exact same one for the other cell you mentioned.
> 5 votes

# Answer
The formula your slide shows describes a dynamic programming with a general gap cost. Finding the optimal solution takes $O(n^3)$ time. No one is using the algorithm these days. This is the first time I see an attempt to fill a matrix with this algorithm. It is interesting, but I am not sure how it works. Although affine gap cost is a special case of general gap cost, the algorithm to find the best score is very different, so is the time complexity. With modern affine gap formulations, each cell holds three values:

1. The overall optimal score (or the optimal score if this cell is a match, which yields an equivalent but different formulation)
2. The optimal score when the cell ends at a deletion
3. The optimal score when the cell ends at an insertion

The values in each cell only depends on the up, left and up-left cells in the matrix. You fill the matrix from the top left corner to the bottom right. This is an $O(n^2)$ algorithm, because you traverse each cell once. For more details, you can read this note. I would recommend to learn affine gap cost the normal way.
> 2 votes

---
Tags: sequence-alignment
---


# Question
Title: How are centromeric/telomeric regions categorized in RepeatMasker? "Other" category?
I haven't been able to find the following in the documentation at http://www.repeatmasker.org/

(1) Why are certain regions classified as "Other"? Are these regions impossible to concretely classify due to alignment issues?

(2) I believe that centromeric and telomeric regions would be classified as tandem repeats. Are telomeric/centromeric regions classified as "DNA" within the hg19 version of RepeatMasker?

Are these documented somewhere?

# Answer
I'll leave question 1 to someone else since I don't know.

Regarding question 2, centromeric and telomeric regions are hard masked in hg19, so there's no sequence to repeat mask.

For hg38, the centromeres are labeled `ALR/Alpha` with a class of `centr`. Telomeres are still hard-masked, but they tend to have simple repeats next to them (as annotated by repeatmasker, which makes sense).

Apparently the `Other` class in hg19 is comprised of SVA elements. In hg38, these are reclassified as `Retrotransposon` class and `SVA` family, which didn't exist in hg19.
> 3 votes

# Answer
On (2), centromeres and telomeres primarily consist of satellite sequences. In human centromeres, the dominant type is alpha satellite (aka alphoid). GGAAT micro-satellite (aka SATIII centromeric micro-satellite) is also very common. There are less common types such as beta satellite. Because there are all kinds of centromeric/telomeric sequences of different lengths and properties, RepeatMasker doesn't put them into one class. Note that RepeatMasker has a "centr" class but it does not include all centromeric satellites (no SATIII, for example). In general, it is not advised to find centromeres with RepeatMasker output.

GRCh38 is an unusual reference genome in that its centromeric sequences are computationally generated modeled after Venter as I remember. They are mostly alpha satellite, but are not real sequences (which has always bugged me). GRC know exactly where they put these sequences. GRCh37 is a more typical reference genome. It lacks probably majority of alpha satellites, but still has many. When you map reads to GRCh37, you occasionally see ~10,000-fold read depth even though the average depth is ~30X. These regions are mostly alpha satellites. Mappers just pile all different copies of satellites onto the ones remained in GRCh37.
> 5 votes

---
Tags: repeat, centromere, telomere
---


# Question
Title: Error in creating a volcano plot in MATLAB
I am a complete newbie to MATLAB and the first task I have is to create a volcano plot. I have been using this link to understand about it and get started.

I tried to run it on dummy values -

```
a=[1 2 3]

b=[4.6 2.7 4.5]

c=[0.05 0.33 0.45]

```

And then I ran -

```
SigStructure = mavolcanoplot(a, b, c)

```

My understanding is that `a` represents the gene expression values for condition 1, `b` for condition 2, and `c` is the list of p-values for the 3 data points in `a` and `b`.

However running this code gives me the error -

```
Index exceeds matrix dimensions.

Error in mavolcanoplot (line 127)
appdata.effect = X(paramStruct.goodVals) - Y(paramStruct.goodVals);

Error in volc (line 4)
SigStructure = mavolcanoplot(a, b, c)

```

Can anyone explain where I am going wrong?

# Answer
> 3 votes
I'm not sure if nowadays matlab is still a good option for informatics. I would rather go for a R script to make a volcano plot. If you don't want to have fancy code but functional one, you can even make the foldchange calculations and log2 transformation in excel and just plot the log2 fold change in x axis and the -log10 transform of the p/q values on the y axis. If you need it, I could provide one of my own scripts.

Update: here is the actual script:

```
#load data into the 'res' variable
res <- read.table("genes.txt", header=TRUE)
head(res)

# Make a basic volcano plot with grey dots
with(res, plot(log2FoldChange, -log10(qvalue), pch=19, col="#4A4E4D", xlab="Log2 Fold Change", xlim=c(-8,8)))

# Add colored points: "indianred" if padj<0.0001, orange of log2FC>2, green if both)
with(subset(res, qvalue<.0001 ), points(log2FoldChange, -log10(qvalue), pch=19, col="#17B090"))
with(subset(res, abs(log2FoldChange)>2), points(log2FoldChange, -log10(qvalue), pch=19, col="#F89538"))
with(subset(res, qvalue<.0001 & abs(log2FoldChange)>2), points(log2FoldChange, -log10(qvalue), pch=19, col="#FA6367"))

```

# Answer
> 2 votes
so while I am not so familiar with MATLAB, in the definition of the volcano plot it specifies that *PValues* is supposed to be a column vector or a DataMatrix object. Your *c* vector is a row vector, and this confuses MATLAB.

# Answer
> 2 votes
Is there a reason why you are doing this in MATLAB instead of R?

In R with the ggplot2 package you can make a volcano plot pretty easily using something like this (may be errors as it is off the top of my head and untested):

```
# load libraries
library(ggplot2)
library(data.table)

# create a dataframe/data.table of the data
a = c(1,2,3)
b = c(4.6, 2.7, 4.5)
c = c(0.05, 0.33, 0.45)
dt = data.table(A=a,B=b,pval=c)

# add column to calculate log2 fold change between a/b
dt[,log2Diff := (log2(A)-log2(B))]
# add column to transform pvalue to -log scale
dt[,logPval := -log10(pval)]

#make the plot
ggplot(dt, aes(x=log2Diff,y=logPval))+geom_point()

```

---
Tags: differential-expression, matlab
---


# Question
Title: Dealing with indels when converting a position to chromStart and chromEnd for a bed file
Starting with a chromosome # and position, I am trying to get chromStart and chromEnd values for a .bed file, but I am not sure how to calculate chromEnd when I have a variant that is an insertion or deletion.

My file is not a VCF file but it gives the chromosome number, chromosomal position (pos), reference (ref) allele, and alternate (alt) allele for each genetic variant, like you would find in a VCF file. The tricky part is that not all of my variants are single nucleotide base changes--I also have insertions and deletions that can be multiple bases long.

For SNPs I am setting chromStart = pos - 1, and chromEnd = start + 1 = pos. (Note, .bed format requires 0-based indexing and the end position is exclusive. See https://genome.ucsc.edu/FAQ/FAQformat.html)

I am unsure of how to handle indels though (insertions and deletions). Currently I am adding the length of the reference allele to the starting position to get the end position (chromEnd = start + len(ref)). I am not sure if this is right though.

I will be using this .bed file with a .bed file from Gencode and bedtools. Does anyone know the standard way the research community (or at least Gencode) treats indels when determining chromosomal position for .bed files? Thanks!

# Answer
You can convert to VCF and use BEDOPS `vcf2bed --insertions`, `vcf2bed --deletions` and `vcf2bed --snvs` to segregate the VCF into three subsets.

For convenience, BEDOPS conversion tools deal with index shifts, including `vcf2bed`. The resulting BED file is 0-indexed, half-open.
> 1 votes

---
Tags: human-genome, bed, bedtools, gencode
---


title: Using large databases with BLAT

text:
# Question
Title: Using large databases with BLAT
I'm a computer scientist working with biologists at a small school that doesn't have dedicated bioinformatics staff. I apologize if I use incorrect terminology since I have limited bioinformatics background.

Our biology staff is studying a soil sample that has been sequenced. The resulting sequences (in the hundreds of millions) were analyzed using mg-rast, but they didn't seem to have a high confidence in the results.

I've been asked to help analyze these sequences. I setup a local instance of NCBI blast, but it was far too slow to analyze that many sequences. We moved on to using BLAT. They've requested we use the nt database from NCBI. I converted the nt FASTA file to 2bit file, splitting the files into small chunks.

gfServer seems to have limitations however. When I run gfServer, I get the following error while gfServer attempts to start.

> Exceeding 4 billion bases, sorry gfServer can't handle that.

Is it possible to use the full nt database using BLAT? Is there an alternative to BLAT that can handle that many sequences?

# Answer
I'd recommend having a go with either OneCodex (commercial / Freemium) or Centrifuge.

Chris Mason's lab has recently looked at metagenomic classification and OneCodex gave the most accurate result.

Centrifuge uses Burrows-Wheeler lookups, but hasn't yet been compared by the Mason Lab.
> 1 votes

---
Tags: sequence-alignment, blast, metagenome, blat
---


# Question
Title: Differences between NetAffx Hg-U133 Plus 2.0 Annotation file versions
From what I was able to observe, the annotation files (CSV) for NetAffx Hg-U133 Plus 2.0 is updated regularly.

However, from glancing at the CSV files (NA35: 2015 & NA22: 2007), I can't seem to notice any glaring difference between the two.

Does anyone know what information is actually being updated? I assume it may be gene identifiers for different databases, but wanted to make sure.

Thanks.

Below are the files I downloaded and examined from Affymetrix website (you need an account to download archived annotation files)

* HG-U133\_Plus\_2 Annotations, CSV format, Release 22 (20 MB, 5/31/07)
* HG-U133\_Plus\_2 Annotations, CSV format, Release 35 (33 MB, 4/16/15)

# Answer
After checking the files myself, it does seem like gene identifiers indeed have been updated. Some that were missing have been added with identifiers and some have additional identifiers now.

Among some other things that were updated from 2007 include:

* Genome version NCBI Version 36.1 =\> GRCh37
* Alignments
* Gene title
* Gene symbol

Thought this would be a good reference for those who may be looking into older researches involving microarrays. I will be sure to check the files more carefully before posting questions like this one.
> 2 votes

---
Tags: bioconductor, sequence-annotation, microarray
---


# Question
Title: How can I calculate the silent mutation for each position in codon?
I am trying to find the number of silent mutations in for each position in codon.

In the **homework**, I am asked:

> For the codons and for each position, count the number of silent mutation. There are 64 codons, for each position there are 3 mutations leading to a change at the DNA level. There are thus 192 possible mutations for each of the three positions of the codons. For each position, print the number of silent mutations.

I understand stand silent mutation is the mutation that doesn't lead to amino acid change. For example, **GCC** can become **GCT** but it is still **Thr**

I am not sure what does "calculating silent mutation for each position" mean? I mean it is easy to find silent mutation for each codon like the example I gave above, but couldn't really understand how to calculate the silent mutation for each position.

Should I change letters (nucleotides) one by one and see if it causing a silent mutation?

# Answer
> 3 votes
`Calculate` in this context would just mean `determine`. You don't actually have to change every base to determine whether a change is silent. Instead, try building a hash table of amino acids to codons, which you can then quickly iterate over to find all silent changes.

BTW, try to avoid loaded terms like "mutation", even if your professor incorrectly used it. Terming things like this as "silent polymorphism" or "silent change" or something along those lines is preferable.

---
Tags: phylogenetics, codon, homework
---


# Question
Title: How to obtain clusters of hierarchical heatmap when using Python?
Is there a good way of obtaining the labels (e.g. genes) within individual clusters that haven been clustered hierarchically in Python (preferentially, but not necessarily, by seaborn)?

I found these difficulties highlighted elsewhere (here, and here), and would hope for some simple option (e.g.: like in MATLAB where cutoff can be specified, and corresponding genes would be returned), without potential risk that the plotted hierarchical heatmap differed from the cluster assignment ( see here)

# Answer
> 6 votes
Clustering like this is typically done with scipy. Here's the code we use in deepTools (original context here):

```
from scipy.cluster.hierarchy import fcluster, linkage
Z = linkage(some_matrix)  # You might want to set `method` and `metric`
groups = fcluster(Z, nGroups)  # You might want to set `criterion`

```

`groups` is then a vector containing the group assignment of the input matrix according to the clustering. As an aside, this is exactly the same as what's shown in the last link in your post, so I'm not sure how you didn't find it.

BTW, at least in python, hierarchical clustering is painfully slow with large datasets.

---
Tags: python, clustering
---


# Question
Title: Best distance parameter for estimating physical interaction between residues in a PDB file
We can calculate the distance between residues in a PDB file regarding different parameters like closest atoms, alpha carbon, beta carbon, centroid and etc. Which one of these parameters are better to show physical interaction between residues in a PDB file?

# Answer
Use a distance cutoff of 12Å between Cᵦ atoms.
> 2 votes

# Answer
A popular (and I would say, respected and trusted) website is PDBsum https://www.ebi.ac.uk/pdbsum (which also has a Wikipedia article about it: https://en.wikipedia.org/wiki/PDBsum)

They measure protein-protein (prot-prot) contacts as any Nitrogen, Carbon, or Oxygen element atoms (N, C, O) within 4.0 Angstroms (Å) when measuring the 3d euclidean distance.

If your atomic contacts are not close enough, the amino acids may not actually interact at all, even though they theoretically could. Since in PPIs you generally want to use reliable ground-truth information, it is better to use cautious parameters. Often, in PPIs, the answer to "does it interact?" is not "yes" or "no" but rather "yes" or "don't know / unknown".
> 6 votes

# Answer
You should read my paper, where I did a comprehensive parameter sweep across different thresholds and definitions and benchmarked the resulting empirical potentials on a decoy ranking challenge:

The relevant phrase from the abstract is "compared using 90 different definitions of residue-residue contact". However, it boils down (as always) to the question, "what do you want to do".

Cite me if you use my paper for your research.
> 1 votes

---
Tags: protein-structure, pdb
---


# Question
Title: Tool for predicting interactions in the cell
What tools are available to predict, based on the structure of a certain protein, its interactions within a cell?

For example, I am considering the split GFP protein and I am trying to predict if there will be any interactions which might hinder its binding efficiency to the other half of the GFP protein.

Note that I am not only looking for protein interactions, but other molecules as well.

# Answer
> 7 votes
You could try one of these tools to predict protein-protein interactions:

1. Struct2Net

> Given two protein sequences, the structure-based interaction prediction technique threads these two sequences to all the protein complexes in the PDB and then chooses the best potential match. Based on this match, the method generates alignment scores, z-scores, and an interfacial energy for the sequence pair. Logistic regression is then used to evaluate whether a set of scores corresponds to an interaction or not. The algorithm is also extended to find all potential partners given a single protein sequence. Further details about the method are described here.
2. PRED

> Pred\_PPI is a web-based system that serves for predicting PPIs from different organisms. This server is freely available to any researcher wishing to use it for non-commercial purposes. Based on auto covariance (AC) and support vector machine (SVM), this tool is capable of predicting PPIs for any target protein pair only using their primary sequences, and assigning an interaction probability to each SVM prediction as well. So the user can use this tool to predict novel PPIs with high confidence.

For protein-RNA interactions, you can try:

1. *cat*Rapid

> Through the calculation of secondary structure, hydrogen bonding and van der Waals contributions, catRAPID is able predict protein-RNA interaction propensities with great accuracy (up to 89% on the ncRNA-protein interaction database, NPinter).

# Answer
> 2 votes
IBIS reports protein interaction with other biomolecules. Your protein or homolog must exist in the Protein Data Bank.

IBIS @ NCBI: https://www.ncbi.nlm.nih.gov/Structure/ibis/ibis.cgi

"For a given protein sequence or structure query, IBIS reports protein-protein, protein-small molecule, protein nucleic acids and protein-ion interactions observed in experimentally-determined structural biological assemblies. IBIS also infers/predicts interacting partners and binding sites by homology, by inspecting the protein complexes formed by close homologs of a given query. To ensure biological relevance of inferred binding sites, the IBIS algorithm clusters binding sites formed by homologs based on binding site sequence and structure conservation."

---
Tags: protein-structure, interactions, 3d-structure
---


# Question
Title: Gene expression signature (AUC) from gene expression data with 48 +48 samples
Does anyone have a good method to find a signature where you combine the expression from multiple genes to predict a specific condition. Assume you have performed sequencing of 48 cancer samples and 48 normal sample. Then find a combination of genes that predicts the cancer samples, using for instance AUC.

# Answer
> 1 votes
I have a method for doing this that uses a MCMC algorithm originally designed for population genetics (via the *structure* program) to generate a statistic between 0 and 1 indicating genetic risk. This needs the use of reference populations for cases and controls that the query individuals are added to for classification. I use T1D as an example. My methods paper is not peer-reviewed, and probably doesn't go into sufficient detail for what you are looking for, but at least indicates one way this can be done:

https://www.biorxiv.org/content/early/2017/01/31/104497

Based on the work I've done with this, using 48 cases and 48 controls with a very specifically defined type of cancer should be sufficient to create a generalised signature from human populations. If a generic cancer profile is needed, the populations are not representative of a general human population, or there are some cases in the control group, then things get trickier and require more individuals. On the other hand, 1500 per group is probably overkill.

# Answer
> 1 votes
There are many different ways of doing it! I would recomend the well established R pacakge CARET which have a whole chapter and many build in functions for it.

In addition you probably want to do some low-level filtering based on variance and fold changes (or the combination - ad DE analysis) to get the number of features down before progressing to the more advanced methods.

---
Tags: r, rna-seq
---


# Question
Title: Downloading CopyNumber450k package from Bioconductor
I am trying to download CopyNumber450k package in the R processor on my laptop in order to see if I could do a step of protocol for my research proffessor. I was unable to successfully download the package I think (because it was too outdated for the newest version of R?) The version of R I have downloaded is 3.4.1. This package appears to work on the version 1.8.0 of R. I thought this was strange, as that seems really old. I then tried to see if this download was possible using some sort of e-version of an R-processor : https://rdrr.io/bioc/CopyNumber450k/ Even the e-version appeared to be unable to find and download this package of R. I have a feeling that I may be misinterpreting something? I tried to find an equivalent package that could run on the newer version of R and didn’t seem to find anything. Then, I tried to see if it was possible to download an old version of R. I couldn’t find a clear way to do that either. Does anyone on bioinformatics stack exchange have insight regarging what I may be doing wrong, something i may be misinterpreting, or how I may overcome this challenge? Thanks

# Answer
I think you need to explain better what you have tried, including some code. You say that you were unsuccessful to download it, because of an outdated version of R. How did you do that?

I am not sure what you have tried, but if you read the instructions for installation here, you need to type the following in R console:

```
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("CopyNumber450k")

```

Show us what kind of error you get from this. Then we can maybe help you further.
> 2 votes

---
Tags: bioconductor, illumina
---


title: How to interconvert InChI and InChIKey?

text:
# Question
Title: How to interconvert InChI and InChIKey?
I would like to retrieve IDs from several databases using `InChI` as an input, e.g.

```
InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1

```

One can use unichem from `bioservices` for this, however, these functions all require the `InChIKey` as input, e.g.

```
KDXKERNSBIXSRK-YFKPBYRVSA-N

```

Is it possible to interconvert the two using `bioservices` and if not is it possible to somehow use the functions in `unichem` with `InChI` rather than the `InChIKey`?

I tried:

```
from bioservices import *
u = UniChem()
u.get_src_compound_ids_from_inchikey('KDXKERNSBIXSRK-YFKPBYRVSA-N')

```

which works fine, however,

```
u.get_src_compound_ids_from_inchikey('InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1')

```

does not work and returns `400`.

If it is not passible via bioservices, is there any other python API available for this?

# Answer
Recently I've developed a python lib for converting InChI to InChIKey: https://github.com/liwt31/chembl\_ikey. Maybe it will help.
> 4 votes

# Answer
Not sure if directly possible in `bioservices` but one can do the following workaround using chemspider:

```
import requests

host = "http://www.chemspider.com"
getstring = "/InChI.asmx/InChIToInChIKey?inchi="
inchi = 'InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1'

r = requests.get('{}{}{}'.format(host, getstring, inchi))
if r.ok:
    res = str(r.text.replace('<?xml version="1.0" encoding="utf-8"?>\r\n<string xmlns="http://www.chemspider.com/">', '').replace('</string>', '').strip())
else:
    print "provide a valid inchi!"

```

This will give the desired InChIKey

```
'KDXKERNSBIXSRK-YFKPBYRVSA-N'

```

which can be used in `unichem`.
> 1 votes

---
Tags: python, kegg, data-retrieval
---


# Question
Title: What sequencing data metrics should I record?
I'm new to bioinformatics and Im starting a new microbial sequencing project and I'd like to record all of the qc data correctly. This is research and I'm not sure what I'll need later. Is there any consensus on useful qc metrics to record in a database for each sequencing run?

No rnaseq, running FastQC, assembling with a defined reference for a single species.

# Answer
> 3 votes
I would say that the following metrics are useful:

* Sample source (Extraction place, person, date)
* Details about pre-HTS preparation (e.g. amplicon PCR, rRNA removal)
* Reads per sample (+ per barcode / lane, if available)
* Sequencing kit (e.g. TruSeq v3), useful for adapter trimming
* Sequence length (e.g. 125bp)
* Paired end / single end
* If paired end, target fragment length in sample prep (e.g. 650bp)
* Reference sequence metadata
    + Scientific name
    + Location of reference genome (preferably a stable web address)
    + Version of reference genome
    + Read mapping rate (per sample / barcode / lane), as a percentage of total reads
    + Mapping program + version
* If available, potential contaminating organisms
    + Name / description
    + Location of contaminant sequence (preferably a stable web address)
    + Estimated mapping rate (bearing in mind that there may be some overlap with the target)
    + Mapping program + version

More generally, ENA / EBI has an excellent reference in sample checklists, which indicate the minimum acceptable metadata for attaching to sequence archive data:

https://www.ebi.ac.uk/ena/submit/checklists

---
Tags: database, ngs, qc
---


# Question
Title: How can I calculate gene_length for RPKM calculation from counts data?
I have read counts data and I want to convert them into RPKM values. For this conversion I need the gene length.

Does the gene length need to be calculated based on the sum of coding exonic lengths? Or are there any different ways for that?

I know that gene length can be taken from the Gencode GTF v19 file. Could you please tell me how that `Gene_length` is calculated?

The Data I'm having is RNA-Seq data. I don't have any idea whether I need to include UTR's in this calculation or only exons?

In Github I have seen RPKM calculation from Counts data with the Gene\_length from Gencode GTF file. Do you think this is the right way of calculation?

And why RPKM is - Its not for differential analysis. For TNBC subtyping they use microarray data. I would like to give a try with RNA-Seq data. So for this I'm trying out different and the right way.

# Answer
> 7 votes
Here you can find some example R code to compute the gene length given a GTF file (it computes GC content too, which you don't need). This uses one of a number of ways of computing gene length, in this case the length of the "union gene model". In this method, the non-duplicated exons for each gene are simply summed up ("non-duplicated" in that no genomic base is double counted). This is a very simple way of getting a gene length.

There are alternative methods that you should be aware of, among which are:

1. Median transcript length: That is, the exonic lengths in each transcript are summed and the median across transcripts is used. This is probably a little more valid than the code that I linked to.
2. Per-sample effective gene lengths: the optimal method, though it requires using something like RSEM, which will give you an effective gene length.
3. Using the length of the "major isoform" in your tissue of interest. This isn't as good as method 2, but is more accurate than all of the others.

At the end of the day, you're just coming up with a scale factor for each gene, so unless you intend to compare values across genes (this is problematic to begin with) then it's questionable if using some of the more correct but also more time-involved methods are really getting you anything.

**Edit**: Note that if you want to plug these values into some sort of subtyping tool (TNBC in your case), you should first start with some samples for which you know the subtype. Then you can at least see if you're getting reasonable results. After that, do read up on how the method works and see if there's anything about RNAseq that makes it incompatible.

# Answer
> 1 votes
I assume you are mapping against the genome rather the transcriptome, since for the later the length would be trivial.

Assuming the first, I think not only the coding sections should be included but also the UTR, since reads can map against them which is what we ultimately care about.

In general, I found gene annotation files (e.g. gff or gtf) can be inconsistent in terms of naming, so it's good practice to inspect and double check. Below is some R code to import the annotation and calculate isoform lengths:

```
library(rtracklayer)
# Reading into a GRanges object
anno <- import.gff3("annotation.gff")
# Filtering exons and UTRs
exons <- anno[anno@elementMetadata$type %in% c("exon","five_prime_UTR","three_prime_UTR"),]

```

Depending on the annotation at hand, the most sensible is probably best to count the length of each isoform which are often contained in the "Parent" column of the annotation file:

```
# splitting up isoforms as preparation for the next step
tmp <- split(exons,as.character(exons$Parent))
# for each isoform, calculate the sum of all reduced exons
Gene_length <- sum(width(reduce(tmp)))

```

Note, reduce merges overlapping intervals together, since UTRs can "contain" bits of exons which would be otherwise double counted.

This code can of course be adapted mainly by changing the "Parent", "exon" etc.

---
Tags: gene, gtf, gencode, fpkm
---


# Question
Title: How can I only get the species name for fasta sequences from blast results?
I am trying to make a phylogenetic tree from sequences obtained with blast. I have files that contains hundreds of fasta protein sequences which are all named like `>NP_567483.1 transcriptional regulator [Pseudomonas aeruginosa PAO1]`. But if I use these files to do multiple alignments or phylogenetic trees using either ClustalX or MEGA7, every sequence are named after the first ID in the sequence names (for this example it would be `NP_567483.1`), and what I want is to have the name of the strain instead (here `Pseudomonas aeruginosa PAO1`) which is part of the sequence description.

So I don't know whether the softwares I'm using understand that the information between the \[ \] is the species and that I just didn't find where to set this as the name or if I have to write some kind of a script to change my fasta files before. I'm really new to this and don't know much about it. Thanks for your help in advance. Julian

# Answer
There might be a "cleaner" way to do this, but what I would do is first remove all of the spaces from the fasta sequence names and replace them with underscores using something like sed.

```
sed 's/ /_/g'

```

Run your BLAST after you remove the spaces and you will get the whole ID string in the result instead of only that first part (e.g. `NP_567483.1`). From here you can parse out the species name, or whatever other ID you would like.

Hope this is helpful! :)
> 3 votes

# Answer
You can replace fasta header of each sequence with the string that is inside the brackets using awk:

```
awk -F '[][]' '{if ($0~/>/) {print ">"$2} else {print $0}}' input.fasta > output.fasta

```

For lines starting with `>` (fasta header), it replaces the line with the string inside the brackets and for other lines (which contain the sequence), it prints the line without making any changes.
> 1 votes

# Answer
With `sed`, first test it:

```
echo '>NP_567483.1 transcriptional regulator [Pseudomonas aeruginosa PAO1]' | sed "s/>.*\[\(.*\)\]/>\1/"

```

Then modify your fasta:

```
sed -i.bak "s/>.*\[\(.*\)\]/>\1/" file.fasta

```
> 1 votes

---
Tags: fasta, blast, sequence-alignment
---


# Question
Title: Obtain CD (cluster of differentiation) marker expression information
There are 394 CD marker genes. Various resources (BD handbook, Abcam guide) contain information about which cell types express which CD markers. What is the best way to obtain or download this data so that it can be used with programs such as `R`?

How can we download similar annotation data on expression for any gene or cell type?

# Answer
> 5 votes
I can offer this as a way to scrape the data from the Abcam resource:

```
library(rvest)
library(dplyr)

antigens <- read_html("http://www.abcam.com/primary-antibodies/human-cd-antigen-guide") %>%
  html_node("table") %>%
  html_table(header = TRUE)

# what targets expressed in "Langerhans" ?
antigens %>%
  filter(grepl("Langerhans", `Cellular Expression`)) %>%
  select(`Target Name`)

   Target Name
1         CD1a
2         CD1b
3         CD1c
4         CD1e
5         CD14
6         CD39
7         CD50
8         CD69
9         CD74
10        CD83
11       CD101
12       CD206
13       CD207

```

The data could be further processed by, for example, splitting the `Cellular Expression` column on either `and` or `,` (comma) to try and make one row per cell type.

---
Tags: r, sequence-annotation, gene-ontology, flow-cytometry
---


# Question
Title: Why do reverse radtags have different start points in radtag sequencing?
Why is it that when doing paired end Radseq that R1 and R2 have different start points? For instance, while doing Radseq you do a restriction digestion of your DNA using enzyme X- shouldn't the sequencing reaction take place from either end of where the enzyme cut? If it doesn't- but happens randomly where do your multiplex sequences go?

For instance- if you are sequencing the sequence on an illumina machine with the "\<\>" representing the cut site-

Where does the sequencing reaction start from in each direction?

Where does the multiplexing sequence go? - does it go on each?

`5-GGCTAGC<CCTGCA*GG>GCGCTT-3
 3-CCGATCG<GG*ACGTCC>CGCGAA-5`

# Answer
Assume that after digestion you're left with the following sequence:

```
5' AGCTAGCTAGCTAGCTGACGTCGTAC 3'
3' GTACGACGTCAGCTAGCTAGCTAGCT 5'

```

Then one read would start with `AGCTAG` and the other would start with `TCGATC`.

The barcodes (what are used for multiplexing) are outside of this and part of the illumina adapters. They're sequences in completely different reads. This is nicely illustrated below:

Index 1 and 2 don't actually need to be present. Even if they are, whether the sequencer actually sequences them is a user setting.

For those curious, the order of the reads as they're sequenced are:

1. read #1
2. Index 1 (optional)
3. Index 2 (optional)
4. read #2 (optional)

If you ever get more than 2 reads from a sequencing provider, then it's likely that they included the index sequence (I have no clue why, but some sequencing centers like to do this...even though it confuses the heck out of their users).
> 1 votes

---
Tags: ngs, radseq
---


# Question
Title: Standard Cutoff for Moderated T-statistics
I'm looking at some microarray data. For the first time I've calculated a moderated T statistic from limma.

Is there any standard practice for where to cut off that value? For log2 fold change I usually cut off at +/- 1.5 and adjust accordingly; for adjusted p value I cut off at p\<0.01 and adjust accordingly.

These moderated T statistics seem a little higher than log2 fold change and I'm not sure what is generally accepted as significant.

Thanks!

# Answer
> 3 votes
You're misinterpreting the moderated T-statistic, it's basically the fold-change divided by its variance. The p-value comes directly from that, so if you filter by moderated fold-change you're just setting an unknown (unless you go through the trouble of figuring it out) p-value threshold. Instead, do exactly as you've been doing and set appropriate (adjusted) p-value and fold-change thresholds.

---
Tags: bioconductor, statistics, differential-expression, microarray
---


# Question
Title: How can I compare structures of proteins responsible for dehydration tolerance in plants?
For my undergraduate research project, I'm trying to compare the structures of some proteins, specifically transcription factors crucial for dehydration tolerance in plants, in the hope of finding differences that make each protein more or less dehydration tolerant than the others. i want to compare the selected transcription factors at their structural level. structure of one transcription factor from different plant species such as rice,sorghum,maize etc. are considered. ultimately i'm hoping to find a relationship between these structures of proteins and their respective level of dehydration tolerance in each plant.

Is there any database or a site where I can find information regarding this?

# Answer
Have you explored in PDB?. If you already have a starting sequence, you can use it as input to search with phmmer using PDB as database. Once you have a curated set to analyse, you can visualize and align the structures with VMD, for example.
> 1 votes

---
Tags: public-databases, protein-structure
---


# Question
Title: Why are inversions defined as the reverse complement and not just the reverse of the reference?
I can’t quite understand the way inversions are defined. In particular I expect an inversion to be only a reversed version of the reference and not its reverse complement. Most sources use diagrams such as below which totally obscure the fact that it is also complemented.

Could someone provide an explanation on what causes an inversion and why it is defined the way it is?

# Answer
If an inversion were just reversed then there would be 3' -\> 3' bonds and 5' -\> 5' bonds. That's why inversions are reverse complemented, you then maintain the normal 5'-\>3' direction. This is nicely covered over on our biology sister site, which also includes this image:
> 5 votes

---
Tags: ngs, structural-variation
---