title: Salmon tximport

text:
# Question
Title: Salmon tximport
I ran bulk `RNA-seq` experiment and got `quant.sf` file. Now, I am struggling with understanding what `tximport` package does and how to use it correctly. My ultimate goal is to feed the data into `DESeq2` for differential expression analysis. In the documentation:

http://bioconductor.org/packages/release/bioc/manuals/tximport/man/tximport.pdf

We have:

```
samples <- read.table(file.path(dir,"samples.txt"), header=TRUE)
files <- file.path(dir,"salmon", samples$run, "quant.sf")

```

I do not understand what are these `samples.txt`. It is not generated by `salmon`, so I suppose that I need to write it myself. How then? What is the correct format for the file? And what if I want to feed just one sample in? That is actually what I want initially: process each sample separately. Could anyone explain step by step how we could utilize `tximport` when we have just one `quant.sf` file to get to raw count matrix to feed it into `DeSeq2`?

# Answer
You are right, `samples.txt` is not generated by Salmon (or any other transcript abundance quantifiers). From documentation, you can find a link to an example of how `samples.txt` should look like.

Note: At the moment, the server is unavailable (I got an `Error 503: Service Temporarily Unavailable`) - I'd suggest you get directly in contact with the owners / authors.

Luckily, there is a call of the function importing `samples.txt` file in the doc itself (variable `samples`), so you can get an idea of ow you have to structure it:

```
> samples

##   pop center                assay    sample experiment       run
## 1 TSI  UNIGE NA20503.1.M_111124_5 ERS185497  ERX163094 ERR188297
## 2 TSI  UNIGE NA20504.1.M_111124_7 ERS185242  ERX162972 ERR188088
## 3 TSI  UNIGE NA20505.1.M_111124_6 ERS185048  ERX163009 ERR188329
## 4 TSI  UNIGE NA20507.1.M_111124_7 ERS185412  ERX163158 ERR188288
## 5 TSI  UNIGE NA20508.1.M_111124_2 ERS185362  ERX163159 ERR188021
## 6 TSI  UNIGE NA20514.1.M_111124_4 ERS185217  ERX163062 ERR188356

```

As an alternate workaround, find the `samples.txt` format directly from tximportData Bioconductor page:

1. download source code
2. extract
3. navigate to: `[...]/tximportData 2/inst/extdata`
4. read content of `samples.txt`.
> 5 votes

---
Tags: rna-seq, deseq2, tximport
---


# Question
Title: Is there a Python/R package with the ability to convert an alignment and reference into a CIGAR?
I'm writing a python function from scratch to do this, but I feel like this must exist in some standard bioinformatics library already. In principle, this is a simply regex operation which many must have written previously. With the goal of having a centralized reference as opposed to everyone writing their own scripts to accomplish the same task (which is time-consuming and not necessarily reproducible), here is the "construct a CIGAR" function:

Here's an alignment:

```
ACGT-TGC

```

Here is a reference:

```
ACGTATGC

```

The reference will be required in order to record indels. One could in principle have a series of "queries", e.g. an R vector or a Python list.

This would generate the following CIGAR:

```
4=1D3=

```

Which R packages/ python library exists to do this operation? I suspect there is something related to Bioconductor as well.

# Answer
Your request is easy to implement. I wouldn't use any libraries in this case. As you haven't showed your implementation, I will provide one:

```
def gen_cigar(ref, qry):
    if len(ref) != len(qry):
        raise Exception('unequal length')
    cigar = []
    for i in range(len(ref)):
        r, q = ref[i], qry[i];
        if r is '-' and q is '-':
            raise Exception('both gaps')
        op = '=' if r is q else 'I' if r is '-' else 'D' if q is '-' else 'X';
        if len(cigar) > 0 and cigar[-1][1] is op: # add to the last operation
            cigar[-1][0] += 1
        else: cigar.append([1, op]);              # a new operation
    return "".join(map(lambda x: str(x[0]) + x[1], cigar)); # turn to string

print(gen_cigar('ACGTAT-CT', 'ACGT-TGGA'));

```

The output is:

```
4=1D1=1I2X

```

PS: It is interesting that you often approach problems from a different angle. The solutions that you think "must exist" actually rarely exist because those are not the ways we (me at least) work with data. This is why I prefer an asker to explain the bigger picture in addition to isolated intermediate problems – there could be simpler and more straightforward solutions without the need of solving the contrived problems.
> 5 votes

# Answer
Here is what I would do in Python. There are a few extra checks I've included, and my definitions of a few cases are different than Heng Li's above, @user172818. The performance is fairly similar.

```
def create_cigar(reference, alignment):
    '''Construct CIGAR from reference and alignment'''
    ## if reference and alignment have unequal lengths, stop
    if len(reference) != len(alignment):
        raise Exception('Error: reference and alignment are of unequal length')
    ## replace any '-' and '*' to '_'
    reference = reference.translate(bytes.maketrans(b"*-", b"__"))
    alignment = alignment.translate(bytes.maketrans(b"*-", b"__"))
    ### creates string with a CIGAR symbol for each position using zip()
    compare_strings = ''.join(['=', 'X'][ref == '_'] if ref == algn else 'X' if ref != '_' and algn != '_' else ['I', 'D'][algn == '_' and ref != '_'] for ref, algn in zip(reference, alignment))
    ## now, parse string into CIGAR format; faster than itertools
    total_match = 0         # total number of adjacent matches
    current_base = compare_strings[0]     # current character being processed
    result = ''     # result of function
    for i in range(len(compare_strings)):
        if compare_strings[i] == current_base:
            total_match += 1
        else:
            result += str(total_match) + current_base
            current_base = compare_strings[i]
            total_match = 1
    result += str(total_match) + current_base
    return result

```

Following the example:

```
print(create_cigar('ACGTAT-CT', 'ACGT-TGGA'))

```

which outputs:

```
4=1D1=1I2X

```

(1) I standard the inputs, such that gaps can be either `-`, `_`, or `*`. (Users could add their own.)

(2) There's no `N`, skipped regions, though based on the use case, I could see users writing this, e.g. a deletion \>= size some number could work.

(3) Hard and soft padding do not make sense in this context. These will be classified as `X`.

(4) By my definition above, if there's a gap in the reference and the alignment, this is a padding `P`.
> 3 votes

---
Tags: r, python, sequence-alignment, bioconductor, cigar
---


# Question
Title: What is a standard approach to binarize microarray gene expression data?
I am using gene expression data from a microarray dataset for the purpose of learning the network structure of genes in the dataset.

Now to do this I need to binarize my data, i.e convert the data for each gene to 0s and 1s. I am new to bioinformatics and am not sure what approach is considered standard for microarray data.

I appreciate your inputs regarding this and I am implementing this in R. So if you know any package that does this please let me know.

# Answer
Based on the info you provide, ArrayBin R package provides you the necessary tools:

* `binarize.array()` from ArrayBin, allowing:

> Implementation of an adaptive method for binarizing gene expression data on a per-probe basis and demonstrate the superior effectiveness of our method when compared with other, commonly used approaches. This adaptive binarization method can be applied to DNA methylation microarray data, which has implications for cross-platform integration, and can reduce batch effects in the data.

In complement, I'd advise you go through some literature and choose the approach which fits you best. Some suggestions:

Maybe you also want to look into more general concepts of transformation prior to microarray data analysis:
> 2 votes

---
Tags: r, microarray, data-preprocessing
---


title: What's a template switching site?

text:
# Question
Title: What's a template switching site?
Reading Islam et al. (2011):

> From each transcript, a single read is obtained, corresponding to a template-switching site located preferentially at the 59 end of the mRNA.

By reading this page I think that a template switching site is an additional set of nucleotides added at the 3' end of the newly synthesized cDNA, that serve to anchor a oligo-nucleotide used by the reverse transcriptase to synthesize the second cDNA strand. It is called template switching because reverse transcriptase switches strand upon reaching this set of nucleotides.

Is it correct?

1. Islam, S. et al. Characterization of the single-cell transcriptional landscape by highly multiplex RNA-seq. genome.cshlp.org 21, 1160–1167 (2011).

# Answer
From the page you cited:

> During first-strand synthesis, upon reaching the 5’ end of the RNA template, the terminal transferase activity of the MMLV reverse transcriptase adds a few additional nucleotide

That usually are CCC (1st is almost always C and the other two most of the times), see: Zajac et al., 2013

These CCCs are used for priming the Template Switching Oligo (TSO), that has GGG and whatever you want to add to the sequence (for instance UMIs and adapters).

So, answering to your question the **polimerase does not switch strands**, it **switches from mRNA as a template to TSO as template**. Because of that is called template switching.

You can see it in the graph from the page you cited. In the template switching step polimerase switch from the mRNA to the TSO and continue with the synthesis of the first cDNA strand using TSO as template.
> 5 votes

---
Tags: scrnaseq
---


# Question
Title: ComBat for batch effects removal on scRNA-seq data
Is it possible to use sva's ComBat for batch effects removal on scRNA-seq data?

Is there any difference between RNA-seq and scRNA-seq data which doesn't allow to use ComBat on single-cell data? (I am thinking about extreme sparsity, \> 80%). If so, are there alternatives?

# Answer
Look at this recent paper that uses ComBat on scRNA-seq data for batch effect removal and states that it "successfully does so".

I also suggest that you check out this publication on Distribution Matching Residual-Nets. Authors evaluated their method also on scRNA-seq data and thus it may be something you are looking for. I personally played a bit with their code and found it straightforward to use.
> 5 votes

# Answer
There are competing claims regarding how well ComBat works on scRNA-seq datasets. A recent paper from John Marioni's introduces a mutual nearest neighbors method that seems to outperform ComBat in at least some relevant scenarios. In general, it's best to look at tSNE or other diagnostic plots after batch correction to see if the results seem reasonable (this advice holds generally).
> 4 votes

---
Tags: rna-seq, scrnaseq, combat, batch-effects
---


# Question
Title: How to get unique somatic mutations for each individual patients
I have exome data 32 patients. My aim is to identify unique somatic mutations in each patients. For this I have completed my pipeline for Whole Exome Data analysis using various tools like BWA mem, picard, gatk etc. Now I have VCF and gVCF file of each patient.

For Variant Quality Score recalibration (VQSR), as per the doccumentation, we should have at least 30 patient data for good recalibration. So in order to perform VQSR, I combined gVCF of all 32 patients using genotypeGVCFs tool of gatk and done with VQSR. But now I have cohort of all patient data (a large combined vcf file).

I don't know how variant sub-class tool or other extraction tool work to extract all of the unique and related features to each patient. I know how to get significant mutations using mutect but here it is combined file but I want to get it for each one individually.

How can I get a vcf file for each patient (without missing a single variant associated to that patient) ?

# Answer
> 4 votes
If what you want is to split the main VCF file into 1 file per sample, you could use bcftools query and view commands.

A similar question was asked on biostars, adapting Jorge Amigo's solution to your situation, you might use:

```
for sample in $(bcftools query -l $combined_vcf)
do
    bcftools view -c1 -Ov -s $sample -o ${$sample.vcf} $file
done

```

---
Tags: vcf, mutations, exome, variation, somatic
---


# Question
Title: PCA on genotype matrix with multiple alleles
Consider an *m x n* genotype matrix of *m* haploid samples and *n* SNPs where each value is an allele encoded by an integer (0,1,2,3).

Is there a good/standard way to encode the alleles in order to perform a PCA on this matrix to investigate population structure ?

I have seen this with matrices of diploid samples where each SNP is encoded as 0/1/2 to represent the number of non-reference allele, but is there a way to consider more than 2 alleles ?

e.g. difference between alleles 1 and 3 should be equal to difference between alleles 1 and 2

# Answer
> 3 votes
You can make a 'dummy variable' for each allele. That means that you don't have info per SNP, but for SNPs with more alleles, the allele is present (1) or not (0).

# Answer
> 2 votes
I recommend the Hamming distance and doing a multidimensional scaling (a procedure similar to PCA but for distances), that way you don't create new variables for the same position.

The distance function can be defined as

```
hamming_dist <- function(x, y) {
    if (x != y) {
        1
    } else {
        0
    }

```

Or you can use a package:

```
library(e1071)
H <- hamming.distance(as.matrix(X))

```

Or as you yourself pointed other implementations.

You can use them to calculate the distance between all samples and then use it for the multidimensional analysis (`cmdscale` in R)

```
plot(cmdscale(H))

```

---
Tags: snp, pca
---


title: tx2gene file for tximport issue

text:
# Question
Title: tx2gene file for tximport issue
I am trying to use `tximport` to get raw count matrices of a bulk RNA-seq sequencing experiment. I followed the tutorial:

https://bioconductor.org/packages/3.7/bioc/vignettes/tximport/inst/doc/tximport.html

Where they generate `tx2gene` in the following way:

```
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

```

So, I reasoned to do it the same way. But since my data is mouse one, I used a different gene annotation file:

```
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

```

The transcript annotation file that I used when running `salmon` is `gencode.vM13.transcripts.fa.sidx`. Then when I try getting `txi` object:

```
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

```

It raises the error:

Seems like I am creating `tx2gene` in the wrong way. Maybe I need to use a different library, not `TxDb.Mmusculus.UCSC.mm10.knownGene`. Any suggestions would be greatly appreciated.

> Update

I looked into the `quant.sf` files that I got from `salmon` and it looked the following way:

A colleague adviced me to strip away everything that follows the dot in the `Name` column, so it started to look:

However, still no luck. The identifiers in the `k <- keys(txdb, keytype = "TXNAME")` look like that:

```
"uc007afg.1" "uc007afh.1" "uc007afi.2" "uc011wht.1" "uc011whu.1"
"uc057aty.1"

```

I suppose I need to convert these identifiers first to `ENS`-like ones, but do not know how to.

# Answer
> 2 votes
The solution was to use a different library as I expected. It should be `library(TxDb.Hsapiens.UCSC.hg19.ensGene)`. And also, as specified in the question, to remove everything after the dot in the `Name` column of the `quant.sf` files.

# Answer
> 1 votes
Do your ids in the expression file (from gencode) match up with the gene IDs in the transcript file (from UCSC?) I'm willing to bet that they don't. Stay consistent throughout your process. Either pull down the gencode transcript to gene translations, or do your pseudoalignments to UCSC transcripts.

---
Tags: rna-seq, salmon, tximport, tx2gene
---


# Question
Title: DeSeqDataSet experimental design: column with integer values
I am creating `DeSeqDataSet` and I am unsure whether I need to include the number of mice into the design formula. Basically, the samples look like this:

I have 4 different groups: `GF`, `B6`, `ET` and `Wld`. It is obvious that these should go into the design formula, but I am not sure about the mice number. Basically for each group, e.g. `GF`, there are 4 experiments with `1`, `2`, `3` and `4` mice each. When I am trying to use both `Groups` and `Mouse` columns as a design formula, `R` is telling that I need to convert `Mouse` to `factor`:

I am a bit confused about two things here: first, as I told already, should I include this `Mouse` column at all in the design formula of `DeSeqDataSet` and second, if I do need it, how to convert it to factor? Any suggestions would be greatly appreciated.

**Update**

It turned out that I was wrong in my initial assumptions: `Mouse#` really was an ‘id’ which I did not include in the design formula. There were 4 groups of 16 mice in each for 4 trials.

What concerns me also is that if I specify design formula as `~ Groups + Mouse`, the DeSeq() function starts comparing the `Mouse` but not the `Groups` which is the most interesting for me.

# Answer
Any IDs (e.g. replicate number) don't need to be included in the design formula if it's not a repeated factor in multiple samples.

In this case, you have said that mouse number is related to the number of mice used in the pool (and is a variable factor that is consistent across multiple samples), in which case it should be included in the design matrix. This would also be the case if the *same* mouse were being used for multiple RNASeq experiments. Converting an integer value to a factor is done by the `factor` method:

```
sampleTable$MiceUsed <- factor(sampleTable$MiceUsed);

```

Other examples of things that might be useful in the design formula would be sex, age, and sequencing run batch (assuming these are not consistent across samples).

With a design formula of something like `~Groups + Mouse`, it will also be comparing Groups, it's just hidden in the results. The DESeq2 vignette is really great at explaining what's going on, but as a quick answer you can pull out the `Groups` results with the `contrast` function. This will pull out group differences that are independent of mouse number:

```
results(dds, contrast=c("Groups", "GF", "Wld"))

```

\[first value is the variable name, second value is the numerator condition, third is the denominator condition\]
> 5 votes

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


# Question
Title: Highly variable genes analysis by DESeq2
In the following article. There are two ways of getting highly variable genes:

```
resSig <- subset(res, padj < 0.1)
head(resSig[ order(resSig$log2FoldChange, decreasing=TRUE), ])

```

Getting genes that are statistically significant, then ordering them by the `log2FoldChange` which is

> The column *log2FoldChange* is the effect size estimate. It tells us how much the gene’s expression seems to have changed due to treatment with dexamethasone in comparison to untreated samples.

And the second one:

```
topVarGenes <- head(order(rowVars(assay(rld)),decreasing=TRUE),20)

```

I am struggling to understand what are the differences between the two ways of getting the most variable genes. What is the difference?

# Answer
The first method (`resSig$log2FoldChange`) is for getting the most differentially expressed genes. These are the genes with the biggest differences between specific groups you pre-define.

The second method (`rowVars(assay(rld))`) is for getting the highly variable genes. These are the genes most variable across all samples regardless of which samples they are.
> 4 votes

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


title: Basic questions about GSEA

text:
# Question
Title: Basic questions about GSEA
I'm currently learning about Gene Set Enrichment Analysis (GSEA) in the hopes of using it in my analysis of differentially expressed genes, and I just had a few questions about the program, specifically about GSEAPreranked, which I need cleared up.

1. On the ranked list needed for GSEA input, should the list include all genes, or only those that pass a certain threshold of significance (i.e. fold change higher than 2, p value less than 0.05, etc.)? Ideally I'd like to sort the genes by fold change alone as I don't trust my p values as much, so should I only include genes with high fold changes?
2. I am comparing multiple conditions of disease with different treatments. Am I correct that GSEA only compares two conditions? If this is the case should I run GSEA for each control/treatment comparison? Would this be conventional?

# Answer
You seem to refer to the GSEA provided by the Broad institute, (there are other GSEA algorithms).

1) You can provide whatever you wish, but if you want to know if those gene sets in which side of the ordered list are they, then provide all the list (of genes) you have.

2) GSEA analyize if the order of a given list distributes in certain way the elements of another (unordered) list. In your case you would need to compute the GSEA for each treatment vs control comparison.
> 6 votes

# Answer
Use the following R package for Gene Set Enrichment analysis of RNA-seq data: seqGSEA

There is another R package (fgsea) recently published called "Fast Gene Set Enrichment Analysis" by Alexey Sergushichev.
> 6 votes

# Answer
GSEAPreranked is a limited version of GSEA. If your sample number is enough(more than 5,6 samples in each condition), you can use GSEA version. Moreover, you can ask this question to the gsea-help google group.
> 1 votes

# Answer
According to the GSEA documentation:

> we are recommending use of the GSEAPreranked tool for conducting gene set enrichment analysis of data derived from RNA-seq experiments

Then more specifically to your question:

> Based on your differential expression analysis, rank your features and capture your ranking in an RNK-formatted file. The ranking metric can be whatever measure of differential expression you choose from the output of your selected DE tool. For example, cuffdiff provides the (base 2) log of the fold change.
> 1 votes

---
Tags: gene, sequence-annotation, rna
---


# Question
Title: How to calculate call detection values in microarray datase if .cel file unavailable
I would like to do a differential gene expression analysis on a microarray data. From the literature, I understood, **call detection (Present, Absent, Marginal)** has to be done to minimize false positives. I found this tutorial for the call detection and also to find the DEGs, **from .cel files.**

For this particular dataset used for my study, **.cel file is unavailable** and only .txt file/ SOFT file is available for download.

1. Can call detection be done from these files?
2. Or is there any other method/ package that can be used for the same?

# Answer
> 0 votes
With things like MAS5, presence/absence calls were done using the P/M/A probe intensities and corresponded more to a "signal is reliable" p-value than literal presence/absence. This data is absent from the SOFT files, so you can't use the standard algorithms.

The only method left is to use "above X percentile of negative probes" as a cut-off for presence, which has been mentioned in a number of papersAFFX-BioB-3\_at (e.g., here). While there are a handful of control probes on that chip, I'm not aware of any of them being negative probes (or if there are any, there don't seem to be enough to be useful for this purpose). Thus I don't think you can even use this method.

Of course you can still do standard differential expression, even with GEO2R in the browser.

---
Tags: bioconductor, file-formats, differential-expression, microarray
---


# Question
Title: Filtering bam file with bedtools intersect using strand specificity
I have a small sample of records from a bam file. I obtained them by intersecting with a bed file not enforcing strand-specificity. So I get them all back upon performing the same intersection again:

```
$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 | samtools view
MN00481:17:000H2FH5N:1:11104:12480:15984    16  III 4499483 42  28M *   0   0   CAGGAACACGCACATCGTCGATTTTCAG    FFFFFAFFFFFFFFFFFFFFFFFFFFFF    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:28 YT:Z:UU
MN00481:17:000H2FH5N:1:11103:9925:11522 16  IV  2388448 42  30M *   0   0   CACTGGAACTGGTCCGATCTTCAAGCCGGT  FFFFFFFFFFA/FFFFAFFFFFF/FFAFAF  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:30 YT:Z:UU
MN00481:17:000H2FH5N:1:22109:11939:5243 16  IV  4749336 42  29M *   0   0   ATGTGGATGCGGTCAGCTTGGGCACGATG   FFFFFFFFFFFF/FFFF6FFFFFFFF/FF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:22104:17385:19431    16  V   4443064 42  29M *   0   0   CAAAATGTGAACAGGAACGTTTTCATGAG   F/FFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:23105:22058:12968    0   X   981524  42  29M *   0   0   CCATCGATTCTCATTTTAATGCTGGTGCT   FFFFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:12104:7867:7500  16  X   4214381 42  29M *   0   0   GCTGCAGGAATGACGGCTCTAGCAAACGT   FFFFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU

```

Now I would like to enforce strand-specificity, but in either case, none of the records pass the filter:

Same strand (`-s`):

```
$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 -s | samtools view
$ # nothing here

```

Opposite strand (`-S`):

```
$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 -S | samtools view
$ # nothing here

```

**What is the correct way of doing this?**

### Edit

Bedtools version:

```
$ bedtools --version
bedtools v2.27.1

```

My bed file has strand information:

```
$ tail features.bed
X   17617826    17619458    "WBGene00015867"    +   .
X   17622777    17640743    "WBGene00019211"    -   .
X   17650177    17651594    "WBGene00044771"    +   .
X   17662457    17665357    "WBGene00003744"    +   .
X   17671028    17680790    "WBGene00020732"    +   .
X   17680910    17682294    "WBGene00044772"    -   .
X   17695231    17701744    "WBGene00017625"    +   .
X   17710485    17711674    "WBGene00019189"    -   .
X   17714455    17714832    "WBGene00007069"    +   .
X   17714967    17718717    "WBGene00007068"    +   .

```

# Answer
> 2 votes
Your bed file switch the score and strand columns. Strand should be column 6, score column 5. For example take this bed file with the 5th and 6th columns switched - zero results for strand specific intersection b/c strand is actually '.'

```
$ cat wrong.bed
chr1    30000000        31000000        noname  +       .
$ bedtools intersect -a mm10_500k_pe.bam -b wrong.bed -wa -u -f 1.0 \
    | samtools view | wc -l
166
$ bedtools intersect -a mm10_500k_pe.bam -b wrong.bed -wa -u -f 1.0 -s \
   | samtools view | wc -l
0

```

After fixing the bed file:

```
$ cat right.bed
chr1    30000000        31000000        noname  .       +
$ bedtools intersect -a mm10_500k_pe.bam -b right.bed -wa -u -f 1.0 \
      | samtools view | wc -l
166
$ bedtools intersect -a mm10_500k_pe.bam -b right.bed -wa -u -f 1.0 -s \
      | samtools view | wc -l
83

```

# Answer
> 1 votes
Your score and strand columns are swapped. If you correct that then you'll get stranded results.

---
Tags: sam, bedtools
---


# Question
Title: plotMA function issue: can not circle genes
I am trying to make `MA`-plot for my `bulk-RNA-seq` dataset and experiencing issues with it. I just copy-pasted the code that can be found here:

https://f1000research.com/articles/4-1070/v2

which is:

```
library("DESeq2")
plotMA(resLFC1, ylim=c(-5,5))
topGene <- rownames(resLFC1)[which.min(resLFC1$padj)]
with(resLFC1[topGene, ], {
  points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})

```

However, I am not getting the right plot with the top gene circled, but the following one:

Maybe it is a bug in the `plotMA` function, not sure. How can I reproduce the figure from the publication?

# Answer
It's likely that the fold-change of the `topGene` is greater than 5. Consequently, the circle is being drawn and the text produced, but they're outside the bounds of the plot (the bounds don't update if you add a circle and text label). Look at the `log2FoldChange` column of `resLFC1[topGene, ]` and increase the bounds to encompass that.
> 4 votes

---
Tags: rna-seq, deseq2, figure-reproduction
---


# Question
Title: Performing t.test to obtain p.values for each gene
I have FPKM values for Normal(80 cases), Group1(135 cases), Group2(147 cases), and Group3(102 cases) for 160 genes (Groups were obtained by clustering tumor data).

I want to perform t-test (two-sided, unpaired, unequal variance) between Normal and Group1, Normal and Group2, and Normal and Group3 such that I get p.value for each gene.
The t.test function in R returns single p.value and not gene-wise (I don't want to do t-test in Excel).

Which package/function should I use to do this in R?

# Answer
You can simply `apply()` `t.test()` to your matrix. In general, though, I expect `lmFit()` from the `limma` package will be simpler to use.

Here is an example with `t.test()`:

```
m = matrix(runif(100, max=c(rep(1, 50), rep(1.5, 50))), nrow=10) # make some data
apply(m, 1, function(x) t.test(x[1:5], x[6:10])$p.value)

```

Neither `t.test()` nor `lmFit()` should be expected to produce the same results as excel, since excel is inferior in every way and should never be relied upon for any statistical analysis (plus `lmFit` is doing something much much more complicated than a simple t-test).

As an aside, FPKMs are highly problematic for use in statistics. If possible, I strongly encourage you to get either per-gene counts or per-gene estimated counts and to use those in `limma` or another differential expression package instead, as the results will be much more reliable. I presume, though, that you're using TCGA data, so you may be stuck with FPKMs.
> 4 votes

---
Tags: r, differential-expression
---


title: What's "signal" in gene expression?

text:
# Question
Title: What's "signal" in gene expression?
I'm exploring gene expression on public available databases. RefExa is one of the databases I've been studying.

When I searched something there, the website reported the expression value unit as "signal". See screenshot where I highlighted:

I'm not able to find anywhere what that really mean. TPM? RPKM? Raw counts? Normalizated for gene length?

My question is that how to interpret gene expression "signal"? What does signal 100 mean?

# Answer
Looking at this resource, these data are microarray data rather than RNAseq data. Thus the "signal" in question is the intensity of florescent light being emitted from the the probe associated with your gene. The relationship between intensity and quantity of RNA is different for every probe because each probe will have a different affinity for complementary fragments dependent on the sequence in question (not least the melting temperature).

Affy intensities are usually presented as log values going up to around 12, but the raw number is next to meaningless and usually only means something as a normalized (e.g. RMA) log fold change.

I've never come across RefExA before, and it appears to be undocumented and unmaintained. I wouldn't rely on it for usable data if I were you. A more modern database of expression data is the EBI Expression Atlas. It presents data from a wide range of different studies, processed and normalised on a unified pipeline and expressed as TPMs.
> 4 votes

---
Tags: gene, database
---


title: Information about control data

text:
# Question
Title: Information about control data
I am putting this question because I did not find any useful information from internet because of limited access. My question is related to control (or normal) data that we use for somatic mutation detection (with the help of mutect or SomaticSnipper etc) with tumor data.

I don't have control data and my aim is to find somatic mutations in Multiple myeloma WES caner data.

Here are the questions:

Can we use Hg19 as a control data? If not, what are the technical reasons of not using hg19 as a control data because it hg19 is consensus call.

Is there any way to infer (or predict) somatic mutations (either statistically or by something else) out of all possible mutations? For this reason, I want to know what information is given in control data.

# Answer
Without control data from your subjects, I don't think there's really no way to distinguish somatic mutations from germ-line mutations. The best you can do is to screen out common variants, which are germ-line mutations that are shared by large numbers of individuals using the population frequencies from something like the Exome Aggregation Consortium:

```
annotate_variation.pl -downdb -webfrom annovar -build hg19 exac03 humandb/
annotate_variation.pl -filter -build hg19 -dbtype exac03 example/ex1.avinput humandb/

```
> 2 votes

---
Tags: database, cancer, mutations, somatic
---


title: BioMart returns one more item

text:
# Question
Title: BioMart returns one more item
I am supplying 30 ensembl ids into BioMart, but somehow it returns 31 result. What is going wrong here? Here you can see the code:

```
library('biomaRt')
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
genes <- rownames(mat)
G_list <- getBM(filters= "ensembl_gene_id", attributes=
c("ensembl_gene_id","entrezgene",
"mgi_symbol","description"),values=genes,mart= mart)

```

genes list is the following one:

```
[1] "ENSMUSG00000097971" "ENSMUSG00000097582" "ENSMUSG00000061808" "ENSMUSG00000097206"
[5] "ENSMUSG00000022425" "ENSMUSG00000094708" "ENSMUSG00000095937" "ENSMUSG00000032484"
[9] "ENSMUSG00000060550" "ENSMUSG00000073409" "ENSMUSG00000026822" "ENSMUSG00000056071"
[13] "ENSMUSG00000035299" "ENSMUSG00000045932" "ENSMUSG00000006154" "ENSMUSG00000032496"
[17] "ENSMUSG00000040809" "ENSMUSG00000032690" "ENSMUSG00000055322" "ENSMUSG00000066443"
[21] "ENSMUSG00000005089" "ENSMUSG00000073489" "ENSMUSG00000005800" "ENSMUSG00000034855"
[25] "ENSMUSG00000029417" "ENSMUSG00000029561" "ENSMUSG00000021843" "ENSMUSG00000020108"
[29] "ENSMUSG00000037849" "ENSMUSG00000066861"

```

As you can see I have here 30 genes, but BioMart returns:

```
> G_list$mgi_symbol
[1] "Slc1a2"  "Mmp8"    "Eps8l1"  "Ddit4"   "Ktn1"    "Enpp2"   "Lcn2"    "Cxcl9"
[9] "Oasl2"   "Ngp"     "Ltf"     "Oas2"    "Cxcl10"  "Mid1"    "Ifi206"  "Chil3"
[17] "Ifit2"   "Tns1"    "S100a9"  "H2-Q7"   "H2-Q7"   "Ttr"     "Gm10163" "Oas1g"
[25] "H2-Q6"   "Ifi204"  "Gm10359" "Gm12671" "Gm26699" "Gm26527" "Gm26917"

```

# Answer
If you print `G_list`, you'll see this half way down:

```
1 ENSMUSG00000060550      15018      H2-Q7      Mus musculus histocompatibility 2, Q region locus 9 (H2-Q9), mRNA. [Source:RefSeq mRNA;Acc:NM_001201460]
2 ENSMUSG00000060550     110558      H2-Q7      Mus musculus histocompatibility 2, Q region locus 9 (H2-Q9), mRNA. [Source:RefSeq mRNA;Acc:NM_001201460]

```

So there are two `entrez` ids that map to the same `Ensembl` id. One of them maps to an alt locus.

If you exclude `"entrezgene"` from the query, you'll get 30 results.
> 6 votes

---
Tags: rna-seq, biomart
---


# Question
Title: Hisat2 : which option should mention for strand specific library read
I was trying hisat2 I get confused by the strand options... The question has been already asked on github here but didn't get any satisfactory answer. I'm surprise how much this strand information could be confusing in different tools.

Looking at the hisat2 manual I get confuse by the two different options:

`--rna-strandness`

> Specify strand-specific information: the default is unstranded. For single-end reads, use F or R. 'F' means a read corresponds to a transcript. 'R' means a read corresponds to the reverse complemented counterpart of a transcript. For paired-end reads, use either FR or RF. With this option being used, every read alignment will have an XS attribute tag: '+' means a read belongs to a transcript on '+' strand of genome. '-' means a read belongs to a transcript on '-' strand of genome.

and `--fr`/`--rf`/`--ff`

> The upstream/downstream mate orientations for a valid paired-end alignment against the forward reference strand. E.g., if --fr is specified and there is a candidate paired-end alignment where mate 1 appears upstream of the reverse complement of mate 2 and the fragment length constraints (-I and -X) are met, that alignment is valid. Also, if mate 2 appears upstream of the reverse complement of mate 1 and all other constraints are met, that too is valid. --rf likewise requires that an upstream mate1 be reverse-complemented and a downstream mate2 be forward-oriented. --ff requires both an upstream mate 1 and a downstream mate 2 to be forward-oriented. Default: --fr (appropriate for Illumina's Paired-end Sequencing Assay).

Should both need to be specified ? Are they redundant ? I didn't find so many person struggling with that, so I guess there is something obvious that I don't get. I hope someone may enlighten me about it !

# Answer
`--fr`/`--rf`/`--ff` should rarely be set, since they refer to the relative orientation of reads and basically everything these days is `--fr`. Mate-pairs were `--rf`, but you don't see those much any more.

`--rna-strandedness` is a very different option, since it sets how reads are expected to align against genes. Most stranded protocols in use these days follow the dUTP-method, where read #2 in a pair has the same orientation as the transcript from which it arose. So either `R` or `RF` would typically be appropriate, unless the library is unstranded. In practice, I expect this is more useful if you plan to run stringTie downstream, since then the `XS` auxiliary tag is set appropriately.
> 4 votes

---
Tags: read-mapping
---


# Question
Title: How to filter Ensembl cDNA and ncRNA FASTA files by primary assembly?
I'm currently performing differential expression analysis using alignment-free quantification using Kallisto. To do this, I need to create a Kallisto index using Ensembl's cDNA and ncRNA annotations available at the following two links:

ftp://ftp.ensembl.org/pub/current\_fasta/mus\_musculus/cdna/Mus\_musculus.GRCm38.cdna.all.fa.gz ftp://ftp.ensembl.org/pub/current\_fasta/mus\_musculus/ncrna/Mus\_musculus.GRCm38.ncrna.fa.gz

However, later on in my analysis I realized that I had multiple Ensembl gene IDs pointing to the same gene. For example, I was gettign results for both ENSMUSG00000026012 and ENSMUSG00000102412 which both map to Cd28. I realized that the reason for this error was due to the fact that the annotations I was using contained haplotypes.

I would like to remove these haplotypes, and to do this I believe that I need to rerun my analysis using a primary assembly. However, Ensembl doesn't have annotated cDNA and ncRNA FASTA files of the primary assembly. How can I filter these files to contain only genes present in the primary assembly? Is there a better way for me to solve this issue?

# Answer
The simplest method is to download the GTF file for GRCm38 and filter that. You can then use one of the many tools out there (bedtools getfasta (make sure to enforce strandedness and possibly use the `--split` option, if you happen to use a BED12 file), gtf\_to\_fasta from tophat2 (it's the only reason to keep tophat2 installed), etc.) to make a transcriptomic fasta file.

Alternatively, Ensembl BioMart can be used to directly download a transcriptomic fasta file using a filter to select only those transcripts on pre-defined chromosomes. Other filters can be applied as well. If you don't already have a GTF file downloaded and one of the previously mentioned packages installed, then this may be the simplest way to get what you want.
> 1 votes

---
Tags: ensembl, filtering
---


title: Deseq2 output issue

text:
# Question
Title: Deseq2 output issue
3 Control replicate and 2 test sample ,so Im comparing my Control which is my Progenitor cells to Monocyte which is mature stage

this is my code im using

```
condition <- factor(c(rep("Control", 3),rep("Test", 2)),levels=c("Control", "Test"))

(coldata <- data.frame(row.names=colnames(countdata), condition))
dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition)

```

This is my deseq2 output

```
gene baseMean log2FoldChange lfcSE  stat pvalue padj CB1 CB2 CB3 GD5_1  GD5_2
ENSG00000144655 2030.6146124999 -5.0812043868   0.1183868766    -42.9203348727  0   0   3327.4890036209 3267.1877800456 3362.7407570783 98.5079578403   97.1475639145

```

Im not sure how come the pvalue and padj is 0? even though the sample value there is a difference and there is also fold change ,something wrong with the result ?

# Answer
A p-value of 0 is expected in this case. The p-values are calculated by comparing the Wald stat to a standard normal, so you can confirm that this is correct like this:

```
> pnorm(-42.9203348727)
[1] 0

```

It really means p \< 2.225074e-308 See here for a more details
> 2 votes

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


# Question
Title: What could be the reason for the samples not clustering?
I'm performing RNA-seq analysis. I have used Hisat2 for aligning reads to the genome and stringtie for quantification and extracted read count information directly from the files generated by StringTie using "prepDE.py" mentioned in Stringtie manual. With the counts data I have used edgeR for Differential analysis.

Before that I have seen MDS plot showing in the following way:

I see that three samples are not clustered and thought that they are not sequenced well. So I removed the three samples and again made an MDS plot.

Can anyone tell me what could be the reason for this type of clustering? Is that because of sequencing? or using count data directly from stringtie or anything else?

# Answer
Samples will only cluster by experimental group if the experimental effect is large enough that it's the primary source of variance between your samples. If that's not the case then you'll get results like you're seeing. Whether that's a problem or not ends up depending on how large of an effect you expect to see. If your groups are different stages in differentiation or another large effect, then you expect that samples will cluster into groups. If you're looking at something like some large unaffected tissue in a disease vs. normal state then seeing a cloud of points is quit reasonable.
> 2 votes

# Answer
The largest variation differences in total expression will be seen in cell population differences. It's a good thing that your conditions aren't clustering on the first and second principal components when all genes are considered, because it means that cell population is probably not a factor in the variation of expression.

There may be some other principal component that *does* separate the conditions, but you might need to drill down fairly deep into the data. This will depend on how much the treatment conditions contribute to the total expression variation. We got lucky with our Th2 RNASeq study in that cell population was mostly on the first and second principal components, and treatment condition appeared on the third principal component.

If you want to force condition clustering, identify the genes that are differentially expressed, and only use those genes for the principal component analysis.
> 0 votes

---
Tags: rna-seq, r, rna, clustering, edger
---


title: BioMart does not return a result

text:
# Question
Title: BioMart does not return a result
I am supplying to `BioMart` `81` `ensemblID`s, but getting `80` results where `1` `ensembleID`: `ENSMUSG00000083840` is absent. What does it mean if `BioMart` can not find the gene? Does it mean that it is a non-coding one or that there is just no symbol attached to it? Here is my `BioMart` query:

```
library('biomaRt')
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
genes <- rownames(mat)
G_list <- getBM(filters= "ensembl_gene_id", attributes=
c("ensembl_gene_id","mgi_symbol","description"),values=genes,mart= mart)

```

# Answer
> 9 votes
Whenever you are wondering about things like this, just look up the identifier on the Ensembl web page. If you look up `ENSMUSG00000083840`, you will see:

> This identifier is not in the current EnsEMBL database

In this case, it doesn't seem to have been replaced by anything and it wasn't a protein coding gene, so it looks like it was simply a wrong annotation which has been removed.

---
Tags: r, biomart
---


# Question
Title: How are Principal Component analyses and Admixture analyses from a genetic alignment different?
How are Principal Component analyses and Admixture analyses from a genetic alignment different?

My understanding is that a PCA will take raw genetic differences across the entire alignment and plot them using dimensionality reduction techniques (single value decomposition for instance) whereas the Admixture analysis will plot likelihoods of sharing certain amounts of ancestry, based on allele frequency data. What I don't understand is how the two are fundamentally different. That is, isn't every variable position found in a PCA going to represent an allele used in Admixture?

# Answer
The fundamental difference lies mostly in the math.

An admixture analysis involves assuming that genotypes (or more likely, genotype likelihoods) in an unknown sample can be modeled with Hardy-Weinburg equilibrium as arising from the combination from two or more source pools of background variants. The output is then an admixture proportion that dictates the mostly percent contribution from each of the source populations (e.g., 23% Irish, 15% Norwegian, 30% Northern African, ...). There are versions of this specifically for handling pecularities of NGS datasets.

A PCA is a very very generic method that uses the same input allelic information, but simply finds **orthogonal** projections of that onto lower dimensions that maximize the variance on the first principal component. You've probably seen plots like the following (from this paper about 1000 genomes projects data):

That allows you to get an idea about how populations are structured, but they're not so useful for determining how individual samples are most likely related to the various populations.
> 3 votes

# Answer
Just to add to the great summary by Devon Ryan: admixture analysis tools are much more flexible than PCA (which is just a fixed mathematical operation), so they can be designed to incorporate LD patterns, phase information, different population evolution models etc.

Here's a paper from the creators of fineSTRUCTURE and Chromopainter admixture tools detailing the action of both - it seems they provide some comparisons with PCA as well: http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002453
> 4 votes

---
Tags: sequence-alignment, pca
---


title: Why chose CNN for a variant caller

text:
# Question
Title: Why chose CNN for a variant caller
Google released their variant caller DeepVariant which won the highest SNP performance award in the Precision FDA Truth challenge (99.999% accuracy).

From the linked github repo, we see that DeepVariant is a CNN, we provide images of the reads as an input and train the neural network with a stochastic gradient descent algorithm.

To me this is unintuitive. Why use images of reads when you already have the reads data?

# Answer
> 4 votes
It's more "images of reads" rather than actual images of reads. In reality, they're feeding in a pileup, which is just a matrix of letters or numbers (you can visualize this as an image, of course). This is more apparent if you look at one of their example python notebooks, CNNs are popular in the image processing field, but at the end of the day images are just large matrices of values so any time you can convert your input into that sort of format then a CNN can be useful (provided you have sufficient training data and some GPUs to throw at things).

# Answer
> 3 votes
> Why use images of reads when you already have the reads data?

Neural networks take one or multiple N-dimension arrays (n-d arrays) as input; they can't take arbitrary forms of data. To feed reads and their alignment into the network, you have to encode them as n-d arrays. If you ignore "colors" for now, a straightforward choice is a (windowLength, readDepth) 2D-array as you see in an alignment browser – it is an image. The "colors" of this image are often called "channels" or "features" in deep learning (DL). They represent base types and base quality in deepVariant. In the end, you get a (windowLength, readDepth, nFeatures) 3D array as input. This is the most natural encoding to DL researchers.

> Why chose CNN for a variant caller

That is because CNN has been repeatedly proved to be the best architecture for such (height, width, nChannels) 3D arrays as input. Read alignments are exactly that.

---

On a historical note, deepVariant was first built on top of google's now-deprecated DisBelief framework that only takes RGB images as input. They had to squeeze all types of information into three color channels. Now with TensorFlow, an image can have an arbitrary number of color channels. This simplifies implementation and improves accuracy. In addition to this image, deepVariant also takes other information as input.

---
Tags: variant-calling
---


# Question
Title: Scaling by linear regression against the number of reads
I am trying to build the preprocessing pipeline presented in The Tabula Muris Consortium et al. (pp).

It is a pipeline to preprocess single-cell sequencing data. There is one step that is not clear:

> Counts were log-normalized (log(1 + counts per N)), then scaled by linear regression against the number of reads (or UMIs), the percent of reads mapping to Rn45s, and the percent of reads to ribosomal genes.

I understand the first part (I assume that log in this context is log2), but I need help on understanding how to *scale by linear regression against the number of reads, the percent of reads mapping to Rn45s, and the percent of reads to ribosomal genes.*

# Answer
I don't know if this question has been solved already, but what they try to do is equalize the depth of sequencing for each cell. Therefore, they scale for the total number of reads. If you regress out (via linear or negative binomial regression) the differences in the number of reads per cell, you end up with cells that have been sequenced with the same depth.

In my opinion they use the ribosomal genes in the same way. They are considered to be a kind of housekeeping genes, which you can use to equalize the sequencing depth.

And yes, the log in this context is log2, used to get fold-change values instead of counts.
> 3 votes

---
Tags: reads, single-cell, scrnaseq, linear-regression, data-preprocessing
---


title: Count files using htseq-count?

text:
# Question
Title: Count files using htseq-count?
I am doing RNA-Seq data analysis using data available on SRA (SRR6047326). Samples are from the pig and they are paired. So, I did alignment using HISAT2.I used pig reference genome available on ENSEMBL. Then, I sorted .sam output files obtained from HISAT2. Now I am trying to get the count file using sorted sam file and it is not working. I am using annotated GTF file for pig which I downloaded from Ensembl.

I get count files using following commands, but all the counts are 0.

```
htseq-counts -i no myfile.sam genes.gtf>output.count

```

What is wrong here ? Are there any alternatives to htseq-count?

# Answer
> 2 votes
I m not sure about htseq-counts but I would suggest go for featureCounts, I will update my answer and share the script/command. Give it a try

the command i use

```
featureCounts -T 50 -C -p -t exon -g gene_id -a /run/media/data2/gencode.v21.annotation.gtf  -o PRIMARY_CELL.txt $(ls *.bam | sort )

```

---
Tags: rna-seq
---


# Question
Title: well identifiers in a 1536-well mtp plate
There's a commonly accepted convention for naming the wells in a 384-well-plate:

```
A01..A24
B01..B24
..
P01..P24

```

It's concise and unambiguous and thus very useful. For 1536-well plates (32x48) it's more difficult because there's more rows than letters in the alphabet. It seems there's no standard. Which nomenclatures are out in the field and which is the best?

# Answer
Here's one way I found:

```
A01..A48
..
Z01..Z48
AA01..AA48
..
AF01..AF48

```

Pros:

* Excel does it like that with columns
* it's clear what's the row and what's the column

Cons:

* alphabetic sorting of well ids does not what you want
* when reading A01 it's not clear whether the plate size is 384 or 1536
> 2 votes

# Answer
and here's a scheme I made up myself:

```
AA01..AA48
..
AZ01..AZ48
BA01..BA48
..
BF01..BF48

```

Pros:

* alphabetic sorting of well ids does what you expect
* it's clear what's the row and what's the column
* it's easy to distinguish well ids of 384-well plate from those of a 1536-well plate

Cons:

* it's not used anywhere AFAIK
> 2 votes

---
Tags: hts
---


# Question
Title: Script to run everything in a loop for extracting tar.gz files into fastq and to bam with alignment?
Original tar.gz files look like below:

```
UNCID_2668056.27d4c881-e906-437f-921d-467eb107798a.150108_UNC15-SN850_0402_AC5FTHACXX_2_GTTTCG.tar.gz
UNCID_2206882.58e94be5-2da2-4a96-944a-060f20216572.150108_UNC15-SN850_0402_AC5FTHACXX_3_CTTGTA.tar.gz
UNCID_2205121.2b3d596a-68ca-4c60-bb16-8c0cbd8fb190.140828_UNC11-SN627_0378_AC5FAWACXX_4_TAGCTT.tar.gz
UNCID_2209082.76a625e9-22b3-42b8-8813-27c133d0b248.110608_UNC11-SN627_0103_BD0CVVABXX_7.tar.gz

```

Before extracting the tar.gz files, I changed their names based on names.text file as below and using the command:

```
xargs -n2 mv -i <<<"$(sed '1d; s/$/.tar.gz/' names.txt)"

Tar.gz_filenames       Sample_Names
UNCID_2668056.27d4c881-e906-437f-921d-467eb107798a.150108_UNC15-SN850_0402_AC5FTHACXX_2_GTTTCG.tar.gz                Sample2
UNCID_2206882.58e94be5-2da2-4a96-944a-060f20216572.150108_UNC15-SN850_0402_AC5FTHACXX_3_CTTGTA.tar.gz                Sample1
UNCID_2205121.2b3d596a-68ca-4c60-bb16-8c0cbd8fb190.140828_UNC11-SN627_0378_AC5FAWACXX_4_TAGCTT.tar.gz                Sample3
UNCID_2209082.76a625e9-22b3-42b8-8813-27c133d0b248.110608_UNC11-SN627_0103_BD0CVVABXX_7.tar.gz                       Sample4

```

After changing the names, looks like below:

```
Sample1.tar.gz
Sample2.tar.gz
Sample3.tar.gz
Sample4.tar.gz

```

Using the tar.gz files after changing names I used the below code:

```
for i in /doc/*.tar.gz
do
    tar xvzf $i -C $TMPDIR
    for sample in $TMPDIR/*1.fastq
    do
        dir2="/destination/directory"
        base=$(basename $sample ".1.fastq")
        module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $TMPDIR/${base}.1.fastq -2 $TMPDIR/${base}.2.fastq | samtools view -Sb - > ${dir2}/${base}.bam
    done
done

```

When I run the above code, the output is like below:

```
150108_UNC15-SN850_0402_AC5FTHACXX_GTTTCG_L002_1.fastq
150108_UNC15-SN850_0402_AC5FTHACXX_GTTTCG_L002_2.fastq
150108_UNC15-SN850_0402_AC5FTHACXX_CTTGTA_L003_1.fastq
150108_UNC15-SN850_0402_AC5FTHACXX_CTTGTA_L003_2.fastq
140828_UNC11-SN627_0378_AC5FAWACXX_TAGCTT_L004_2.fastq
140828_UNC11-SN627_0378_AC5FAWACXX_TAGCTT_L004_1.fastq
110608_UNC11-SN627_0103_BD0CVVABXX.7_1.fastq
110608_UNC11-SN627_0103_BD0CVVABXX.7_2.fastq
Warning: Could not open read file "/scratch/username/slurm-job.36320423/140828_UNC11-SN627_0378_AC5FAWACXX_TAGCTT_L004_1.fastq.1.fastq" for reading; skipping...
Error: No input read files were valid
(ERR): hisat2-align exited with value 1
Warning: Could not open read file "/scratch/username/slurm-job.36320423/150108_UNC15-SN850_0402_AC5FTHACXX_CTTGTA_L003_1.fastq.1.fastq" for reading; skipping...
Error: No input read files were valid
(ERR): hisat2-align exited with value 1
Warning: Could not open read file "/scratch/username/slurm-job.36320423/150108_UNC15-SN850_0402_AC5FTHACXX_GTTTCG_L002_1.fastq.1.fastq" for reading; skipping...
Error: No input read files were valid
(ERR): hisat2-align exited with value 1

```

I dont mins whatever the name of fastq's in `$TMPDIR` but from `$TMPDIR` hisat2 should use those fastq's and give output bam files in /destination like below:

```
Sample1.bam
Sample2.bam
Sample3.bam
Sample4.bam

```

I would like to run a script doing tar.gz files -\> fastq -\> bam (everything in loop)

# Answer
This should work for you:

```
for i in *.tar.gz
do
    dir2="/destination/directory"
    tar xvzf $i -C $TMPDIR
    files=($TMPDIR/*)
    module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $files[0] -2 $files[1]| samtools view -Sb - > ${dir2}/$i.bam
done

```

(If this is the sort of thing you'll be doing often, it would be worth investing the time to learn a more robust approach such as Snakemake or Nextflow).

Edit:

If you want to run on all files at the same time, the easiest way is `find` / `xargs`

Change you script to this:

```
for i in $@
do
    dir2="/destination/directory"
    tar xvzf $i -C $TMPDIR
    files=($TMPDIR/*)
    module load HISAT2/2.0.4-goolf-1.7.20; module load SAMtools/1.3.1-goolf-1.7.20; hisat2 -p 8 --dta --rna-strandness RF -x /genome_index -1 $files[0] -2 $files[1]| samtools view -Sb - > ${dir2}/$i.bam
done

```

Then run it like this:

```
find . -d 1 -name '*.tar.gz' | xargs -n 1 bash NAMEOFSCRIPT

```

This works for SGE. I'm not sure about the details of SLURM. You might have to figure out how to run each job in the background or use `gnu parallel` or something. As `@Ian Sudbury` says, pipeline software can take care of all these details for you.
> 7 votes

# Answer
Your samples end with `_1.fastq` and `_2.fastq`, so you need to change this:

```
-1 $TMPDIR/${base}.1.fastq -2 $TMPDIR/${base}.2.fastq

```

to this:

```
-1 $TMPDIR/${base}_1.fastq -2 $TMPDIR/${base}_2.fastq

```
> 1 votes

---
Tags: rna-seq, file-formats, fastq, linux
---


# Question
Title: running a samtools command for multiple bam files from 1000 genomes project
I want to compute the depth of coverage only for specific intervals in phase 3, 1000 genomes project.I have not worked with 1000 genomes project before, so a bit unfamiliar with it. I do not want to download all of the bam files for the entire genomic region, just those intervals which I am interested in them. I guess I can use this (I found from biostars forum)

```
samtools -bu view 'http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam' "2:1000000-2000000" | (...)

```

However, I am just a bit confused. First of all, when I try the code above directly from the command line, I get an error message "No such file or directory"! Am I doing something wrong? or should I provide a specific path? Second, I want to run it for the mapped files for all individuals: `*.bam”.` through a loop.

Can anyone help me to solve my problem?

# Answer
The command you're looking for is:

```
samtools view -c http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam 2:1000000-2000000

```

For a series of intervals, you'd be better off scripting something with pysam:

```
import pysam
bam = pysam.AlignmentFile("some file.bam")
regions = open("your regions.bed")
for region in regions:
    chrom, start, end = region.split("\t")
     bam.count(chrom, int(start), int(end))  # Do something with the result

```

You might see if there are bigWig files available, they'd be faster to process.
> 1 votes

---
Tags: sam, samtools
---


# Question
Title: Obtaining identifier from plot of sequence GC%
I was looking at this script from https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc297 and wanted to use it to determine outliers for tentative Horizontal gene transfer.

The data set is provided in the example, under the section **20.2.2 Plot of sequence GC%**.

I want to yield a identifier for the x axis instead of just an arbitrary number of genes or cds' so I can actually determine which sequences have the outlying values, as I can not with the current graph.

I am not great at python or understanding list comprehensions, so if anyone could augment the code provided in the example to do so, it would be appreciated.

I provided the code below also.

```
#from : https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc297, data set under
#section 20.2.2 Plot of sequence GC%, ls_orchid.fasta

from Bio import SeqIO
from Bio.SeqUtils import GC

gc_values = sorted(GC(rec.seq) for rec in SeqIO.parse("ls_orchid.fasta", "fasta"))
import pylab
pylab.plot(gc_values)
pylab.title("%i orchid sequences\nGC%% %0.1f to %0.1f" \
            % (len(gc_values),min(gc_values),max(gc_values)))
pylab.xlabel("Genes")
pylab.ylabel("GC%")
pylab.show()

```

# Answer
> 2 votes
The tricky bit here is that you need a list of sequence records sorted in the same order as the CG values. You can accomplish this by making tuples of the seq records and GC values, then sorting the tuples specifying the GC values as the sorting key:

```
gc_w_records = sorted(((GC(rec.seq), rec) for rec in SeqIO.parse("ls_orchid.fasta", "fasta")), key=lambda tup: tup[0])
gc_values = [x[0] for x in gc_w_records]
records = [x[1] for x in gc_w_records]

```

I don't know much about `matplotlib` and it's going to be challenging to display 100 IDs, but you'll probably be interested in the extremes of the GC distribution, which you can obtain directly from the sorted list of records:

```
records[0:2]
# [SeqRecord(seq=Seq('CGTCACGAGGTCTCCGGATGTGACCCTGCGGAAGGATCATTGTTGAGATCACAT...CAT', SingleLetterAlphabet()), id='gi|2765587|emb|Z78462.1|PSZ78462', name='gi|2765587|emb|Z78462.1|PSZ78462', description='gi|2765587|emb|Z78462.1|PSZ78462 P.sukhakulii 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]),
# SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACAT...AGG', SingleLetterAlphabet()), id='gi|2765568|emb|Z78443.1|PLZ78443', name='gi|2765568|emb|Z78443.1|PLZ78443', description='gi|2765568|emb|Z78443.1|PLZ78443 P.lawrenceanum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])]

```

or:

```
records[-3:-1]
# [SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAGCAG...TAA', SingleLetterAlphabet()), id='gi|2765656|emb|Z78531.1|CFZ78531', name='gi|2765656|emb|Z78531.1|CFZ78531', description='gi|2765656|emb|Z78531.1|CFZ78531 C.fasciculatum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]),
# SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGCTCATTGTTGAGACCGCAA...AAG', SingleLetterAlphabet()), id='gi|2765625|emb|Z78500.1|PWZ78500', name='gi|2765625|emb|Z78500.1|PWZ78500', description='gi|2765625|emb|Z78500.1|PWZ78500 P.warszewiczianum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])]

```

# Answer
> 2 votes
There are various attempts to guess putative horizontal gene transfer using sequence composition (kmer-frequencies or GC content), you do not have to implement a new method by your own.

I co-authored one of them called SigHunt. It's kind of working, but the code is really messy and there are now bit more precise methods out there. The advantage is SigHunt is that it provides access to all the data calculated on the way (the tetranucleotide frequencies), so if you want to play with it, it's fun.

The method I would probably go for now if I would have searched for HGT is MTGIpick. The most of the tools detects HGT within a genome by a sliding window, but MTGIpick is refining boundaries where the change of kmer frequencies is rapid. I never tried it personally but the approach is sound to me.

Alternatively there are other paradigms for detection of HGT. Phylogeny based methods (there are many, but for instance) or by mapping reads to species boundaries. Phylogeny based methods usually provide stronger evident, but have a disadvantage to be database dependent (i.e. you can not build a tree if you don't have sufficient coverage of homologs across tree of life).

I would also like to mention that right now I do not follow the HGT literature closely, there might be some big new I just missed.

---
Tags: gene, genome, python, identifiers
---


# Question
Title: Why do NEB adapters have non-complementary sequence?
I am making some infographics of different library prep types and noticed something weird about NEB's Ultra II adapters.

`Molecule 1` is the desired product of the adapter ligation step. `*` indicates a phosphorothioate bond and `-` indicates a normal bond. `|` indicates complementary sequence through hydrogen bonding.

`Molecule 2` is the product after treatment with the USER enzyme.

`Molecule 3` is the result of melting molecule `2` and hybridizing the NEBNext Universal PCR Primer for Illumina.

My question is: Why does `Molecule 3` have a large non-complementary region on the 3' end of the original molecule? Is that just to force `Molecule 1` to have an open loop conformation so that the USER enzyme can access it better? Or maybe it is something to do with PCR kinetics?

I understand that PCR synthesizes in the `5'->3'` direction, so `Molecule 3` only forms with the original adapter-ligated molecules before PCR.

# Answer
> 1 votes
Module 3 does not exist during PCR amplification. Here is how NEBNext adaptor and multiplex PCR primers work.

During the 1st PCR cycle, only the P7 primers are used to bind and amplify ligated DNA, generating the complementary strands that contain the binding sites for P5 primers (Universal in this case). During the second cycle of amplification, the P5 primer hybridizes to the DNA strand generated during the first cycle and start to amplify. Please refer to the following figure for more details.

# Answer
> 1 votes
The original point of the non-complementary sequence is to allow formation of a Y-adapter after USER treatment, since that's what's needed for the Illumina chemistry to work. In reality the bottom strand may not be there in step 3, since USER has already digested one strand (creating the hairpin in step 2) and size selection has been performed. Regardless, it'll get removed once cluster formation is performed.

---
Tags: phylogenetics, adapter, pcr, library
---


# Question
Title: Select top 100 genes ranked by variance in read counts
I have a gene expression data (read counts) table with many samples. Showing some of them below:

```
Ensembl_ID          A   B   C   D   E   F   G
ENSG00000236601.2   0   0   0   1   0   0   1
ENSG00000237094.12  2   187 442 109 88  144 486
ENSG00000269732.1   0   0   0   0   0   0   0
ENSG00000284733.1   0   0   0   0   0   0   0
ENSG00000233653.3   0   0   0   0   0   0   0
ENSG00000250575.1   1   11  17  5   5   8   17
ENSG00000278757.1   0   0   2   0   1   0   7
ENSG00000230021.9   0   25  30  27  9   8   11
ENSG00000235146.2   0   0   0   0   0   0   0
ENSG00000225972.1   0   7   3016    11  5   14  5
ENSG00000225630.1   1   113 194 76  91  47  94
ENSG00000237973.1   25  1037    9767    1160    886 1321    5220

```

I want to select top 100 genes with reverse sort of variance. How to do that?

# Answer
> 4 votes
I assume by "reverse sort of variance" you mean "highest variance". Assuming you made a matrix out of that (set the row names to the first column and then remove it) and called it `m`:

```
sel = order(apply(m, 1, var), decreasing=TRUE)[1:100]

```

`sel` then contains the indices into your matrix (or the original dataframe).

BTW, I hope you're not filtering by variance before performing differential expression. That would be incredibly incorrect. Doing this sort of filtering for performing PCA or something like that can be useful, though you should really be using normalized counts (or a transformation of them).

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


# Question
Title: Regressing out unwanted sources of variation in single cell RNA-seq data
I have a dataset of single cell count data and I want to regress out the variation caused by the number of UMI's and the percentage of mitochondrial genes.

I know that count data is discrete data, and commonly follows a negative binomial distribution. To regress out these two confounders, should I use a linear model or GLM (poisson/negative binomial)? And how can I determine the optimal choice?

# Answer
Cleaning data before doing analysis can be more important than the method of data analysis. For UMIs, there's an obvious cleaning step that can be done: it would be better to filter out duplicates (prior to generating counts) than to try to incorporate that information in the count analysis.

I'm not familiar with single cell analysis, but I am aware that DESeq2 can be adapted to work with single cells, and incorporating additional variation is just a matter of adding an additional covariate into the DESeq2 formula model.
> 2 votes

# Answer
Seurat has as part of its protocol a step where you filter based on UMI counts and percent mitochondrial

http://satijalab.org/seurat/pbmc3k\_tutorial.html
> 2 votes

# Answer
If you just want to filter UMI related problems and cells with high mitochondrial to nuclear ratio, maybe there is no need to use modeling.

* **For UMIs**: You can **filter UMIs** that are **not supported by at least x reads** (for instance 3 reads) and that probably arise from sequencing errors.
* **For high mitochondrial to nuclear ratio**: You can compute (number of molecules of mitochondrial genes / number of molecules of nuclear genes) and plot this ratio for all the cells and try to put a **threshold value for cells that have high values for this ratio**. Thus, filter these cells that behave strange compared to the rest.
> 0 votes

---
Tags: rna-seq, scrnaseq, single-cell, linear-regression
---


# Question
Title: SAM format: Does the BAM "Integer or numeric array" field no longer exist? Why?
The SAM file format specification makes clear that optional fields must follow the format `TAG:TYPE:VALUE`

```
Type   Regexp matching                                       VALUE Description
A     [!-~]                                                  Printable character
i     [-+]?[0-9]+                                            Signed integer
f     [-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?                 Single-precision floating number
Z     [ !-~]*                                                Printable string, including space
H     ([0-9A-F][0-9A-F])*                                    Byte array in the Hex format10
B     [cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+    Integer or numeric array

```

However, looking at the samtools source code, I noticed this:

https://github.com/samtools/samtools/blob/develop/sam\_view.c#L920-L931

```
switch (aux_type) {
    case 'H':
    case 'Z':
        if (kputs(bam_aux2Z(s), linebuf) < 0) return false;
        break;
    case 'i': kputw(bam_aux2i(s), linebuf); break;
    case 'I': kputuw(bam_aux2i(s), linebuf); break;
    case 'A': kputc(bam_aux2A(s), linebuf); break;
    case 'f': kputd(bam_aux2f(s), linebuf); break;
    case 'B': kputs("*** Unhandled aux type ***", linebuf); return false;
    default:  kputs("*** Unknown aux type ***", linebuf); return false;
}

```

Questions:

(1) I have never seen an "Integer or numeric array" of `case "B"`. What does such an tag look like? Is there an example?

(2) Why is this currently `"*** Unhandled aux type ***"`? Was this the case with previous versions of samtools?

# Answer
> 3 votes
1. An example would be `ZC:B:i,271,271,1,0`, which is from TMAP. The `i` part denotes the values held in the array, which are comma-separated.
2. That code is only used in `samtools fastq`, where that tag type isn't output into fastq files.

---
Tags: sam, samtools, htslib
---


# Question
Title: Is there a publicly available tumor-normal sample?
I am looking for a publicly available matched tumor-normal sample. I need Illumina fastq reads (or an aligned bam file, since I could extract the reads from it) from a tumor and a matching, non-tumor control set from the same patient. Ideally, I would also have a truth set of known variants in the sample, but I can live without that.

I checked the Next Generation Cancer Knowledge Network but although they do seem to have appropriate bam files, which I could convert to fastq, all of the ones I found were "controlled access" and getting access seems very complicated.

Is there something like the public NA12878 sample for tumor-normal analyses?

# Answer
> 3 votes
Please take a look at

An open access pilot freely sharing cancer genomic data from participants in Texas

Although I have not worked on this paper but author claimed that both tumor and normal data is publicly available. In case if you still need something more, than you can work on virtual normal data. And there are several papers which provide machnine learning techniques to overcome the deficiency of normal data (not much efficient but still they give considerable results). Hope this information will help you. Thanks

---
Tags: ngs, cancer, public-dataset
---


# Question
Title: Official documentation for textual BLAST output?
This is pretty vexing: I got asked to explain the meaning of the match line in a blastp output, especially the subtle difference between a blank and plus. E.g.:

```
 Query  181  NHIEEL...
             N I+EL...
 Sbjct  181  NQIKEL...

```

For me this is under the "everyone knows this" category, but I wanted to have an authoritative reference I could cite.

Alas, none of the documents on the BLAST help page seems to have this information. My Google fu also failed as I was not able to find an authoritative source which I could cite.

Also, I'm not sure whether BLAST introduced this notation or it simply picked it up somewhere else.

Any pointers appreciated.

# Answer
The reference that is mentioned when I run a BLASTP search using NCBI's web service is this one. Here's some relevant text from the caption of Figure 3 of that paper:

> On the central lines of the alignment, identities are echoed and substitutions to which the BLOSUM-62 matrix (18) gives a positive score are indicated by a ‘+’ symbol.
> 1 votes

# Answer
I think gringer is correct to say this is because of the BLOSUM62 matrix and a ‘+’ symbol inidcating a positive score for a mismatch but here is a little more explanation.

BLASTP looks for (nearly) identical short protein matches pairwise, called seeds, and then extends the seeds in each direction one residue at a time using arbitrary gap penalties for indels and a matrix for substitutions. Based on some thresholds and more, the highest scoring alignments are kept and ultimately reported.

The default choice of matrix for BLASTP is BLOSUM62, which looks like this. Obviously matches have high positive scores, but also mismatches have varying positive or negative scores.

If you look up your changes:

`H -> Q`; ; score 0

`E -> K`; `+`; score 1

Note that this matrix is symmetric so the scores for `H -> Q == Q -> H` and `E -> K == K -> E`.

In terms of why `E -> K` has a positive score, well not all amino acid changes are viewed equally. BLOSUM matrices are not based on any biological theory but instead use the log-odds ratio, a likelihood approach, that estimates from real sequence alignments how frequently different amino acid changes occur. The more often `E -> K` substitutions are observed in these alignments, the higher the substitution will score in the derived matrix.

Hopefully these scores are intupretable based on the biochemical properties of the amino acids. For example, `R -> K` scores 2, a high score for a mismatch, but both amino acids are basic and so share a similar structure, presumably making a change more likely to occur and/or more likely to maintain the protein function and so not be harmful.

For your case, I don't have a good biological explanation for why `E -> K` scores positive, another author notes this is surprising, but ultimately the reason it is given a positive score is just because regardless of why, `E -> K` does apparently frequently occur, at least in the alignment data used to derive this scoring matrix. Incidentally, BLOSUM matrices are pretty ancient and based on few alignments by modern standards, but nonetheless appear to work quite well.
> 2 votes

---
Tags: blast
---


# Question
Title: Error while Importing fastq sequences to jMHC software
I am trying to analyze my .fastq files (MiSeq exon reads) with the jMHC software, but output files I have obtained have no results.

I think I might be doing something wrong during the first step: "connect to database". Because on the "import/extracting" step, I open all my .fastq files, click "Start", and it shows this message:

> "Success! Processed: 139, Imported: 0".

So I assume it is reading the files OK, but not importing them to the "database"? Could it be the extension of my files? The software says it operates on FASTA files. Are .fastq files not compatible?

Anybody that has successfully run this program could give me a hand?

# Answer
> 0 votes
Yes, use .fasta , you can try to convert your fastq files with NextGen Workbench for example.

---
Tags: sequence-analysis, errors
---


# Question
Title: How can I improve ESS values for ucldMean in beast?
Using starBeast I ran a 50 locus tree, with 34 individuals, using evolutionary models indicated by AICc in jModel test. I ran this tree with a MCMC chain length of 500 million, storing every 10,000 generations. I used a relaxed clock log normal.

I get low ESS values for posterior, prior, birthRate.t:Species, YuleModel.t:Species, popMean, and TreeHeight.Species. ucldMean.c:Locus\_name is low for each loci. The TreeHeight statisitcs also have low ESS values. Also rate.c:Locus\_name.mean for each loci is low as well.

Is there a way to change my parameters in Beauti to help increase the ESS values of these statistics? Is there a parameter that would likely cause the ESS values of these particular statistics to be low? Would increasing the MCMC chain length increase the ESS values, given that I've already run half a billion?

# Answer
The official BEAST website suggests a number of solutions to improve the ESS of a parameter

It is not uncommon to have several billions of generations so I don't imagine increasing your number of generations would hurt.
> 1 votes

---
Tags: genomics, phylogenetics
---


# Question
Title: Blast hits disappearing after changing -evalue
I was teaching an introduction to bioinformatics when the students and I noticed strange Blast behavior that we couldn't explain.

With the default `evalue` parameter, our best hit shows an evalue of `3e-06`, while that hit disappears with a evalue threshold of `1e-5` (which is bigger!)

*(after a few tries it seems that the evalue of the top hit starts changing with an evalue threshold of 2)*

**Input proteins:**

```
wget ftp://ftp.ncbi.nih.gov/refseq/B_taurus/mRNA_Prot/cow.1.protein.faa.gz
wget ftp://ftp.ncbi.nih.gov/refseq/H_sapiens/mRNA_Prot/human.1.protein.faa.gz
gzip -d *.gz
head -199 cow.1.protein.faa > medium_cow.faa

```

**Database building:**

```
makeblastdb -in human.1.protein.faa -dbtype prot

```

**Blast command 1 (default):**

```
blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast

```

output:

```
# BLASTP 2.2.31+
# Query: NP_001180091.1 SH2 domain-containing protein 1B [Bos taurus]
# Database: human.1.protein.faa
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 13 hits found
NP_001180091.1  XP_011539315.1  29.41   102     61      3       5       100     144     240     3e-06   42.0

```

**Blast command 2 (-evalue 1e-05):**

```
blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast -evalue 1e-05

```

output:

```
# BLASTP 2.2.31+
# Query: NP_001180091.1 SH2 domain-containing protein 1B [Bos taurus]
# Database: human.1.protein.faa
# 0 hits found

```

**Blast command 3 (-evalue 2):**

```
blastp -query medium_cow.faa -db human.1.protein.faa -outfmt 7 -out test_blast -evalue 2

```

output:

```
# BLASTP 2.2.31+
# Query: NP_001180091.1 SH2 domain-containing protein 1B [Bos taurus]
# Database: human.1.protein.faa
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 6 hits found
NP_001180091.1  XP_005244948.1  40.62   32      18      1       4       35      394     424     0.076   29.3

```

I reproduced this behavior with Blast 2.2.31+, 2.7.1+ and the web blast at https://blast.ncbi.nlm.nih.gov (screenshots attached)

Anything obvious that my students and I are missing? To my knowledge changing the evalue threshold should no change the values themselves, it should just change the reporting.

# Answer
> 3 votes
I am the author of the question / answer here Why does a very strong BLAST hit get lost when I change num\_alignments, num\_descriptions or max\_target\_seqs parameter?

> To my knowledge changing the evalue threshold should no change the values themselves, it should just change the reporting.

That is sadly not the case, though thats what I would say is typically expected.

Yours was the exact situation that lead me to question the implementation of NCBI blast. I got results with an evalue above a certain value with a relaxed threshold, but with a stricter threshold (still lower than the value of the results I got) those results disappeared.

From what I found all parameters imposing some threshold show this behavior and are applied in an early step of the program leading to prematurely removed results that later on might score higher than those thresholds. (Say your hit first has an evalue of 5: threshold 10 = ok and kept, threshold 1e-5 = gets removed, but that same hit gets improved further along in the program and then gets an evalue of 1e-6. Now this is only shown when you use an evalue threshold of 5 or higher, though 1e-6 is clearly also a result you would want to see when using a threshold of 3, 1, 0.01, ... etc.)

I was relatively shocked by this and the fact that its **basically nowhere officially documented and generally a not known fact** what influence these options have.

Sadly for now the stance of the developers is basically all is working as intended and the only feedback I got amounts to

> The safest approach is to relax (increase) the threshold and deal with a larger result set.

aka Suck it up and deal with it.

---
Tags: blast
---


# Question
Title: How to get all PDB homologs from Uniprot (mapping + BLAST)?
I'd like to create a dataset consisting of all sequences which are either present in the PDB, or whose homolog is present in the PDB. In other words, any sequence in the PDB or any sequence related to it. The similiarity margin is to be very wide, so anything above, for instance, 25% sequence identity is accepted (perhaps there are better criteria, but you get the point).

No further limitations, so any organism, with or without experimental data, not limited to the manually curated Swissprot set.

I have tried the SIFTS database, but that contains, as far as I can tell, only mappings for sequence identity \>85%. I would like to go far below that.

# Answer
For the correspondence between PDB and Uniprot entries you can use SIFTS -- a semi-automated mapping between PDB and UniProt maintained by PDBe.

The pipeline that creates the mappings uses BLAST and a few other criteria to decide which UniProt entry should be assigned to each PDB entry. The SIFTS website has all the data in CSV files.
> 4 votes

# Answer
You can download a BLAST db containing all sequences of proteins in Uniprot and in the PDB. The way I would go about this is first download the databases for uniprot and PDB, then query the PDB database for each sequence from Uniprot. If you get a BLAST hit above a certain threshold (whatever you define as homolog) then add that sequence to a file, and voila, you have a full set of sequences that fulfill your criteria. The only tools you would need for this are the BLAST executables (found here) and the databases (info on how to download can be found here).

I should note that you don't really need any Python to get this to work, only some shell scripting. Although you can run these commands through a Python wrapper if you really want to.

Good luck!
> 1 votes

---
Tags: python, public-databases, pdb, sequence-homology, uniprot
---


title: Seurat installation fails

text:
# Question
Title: Seurat installation fails
When I am trying to install `Seurat` in `R markdown file`:

```
install.packages("Seurat")

```

it fails with the following errors:

```
Installing package into ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3’
(as ‘lib’ is unspecified)
Warning in install.packages :
  dependency ‘proxy’ is not available
also installing the dependencies ‘igraph’, ‘dtw’, ‘diffusionMap’

trying URL 'https://cran.rstudio.com/src/contrib/igraph_1.2.1.tar.gz'

Content type 'unknown' length 2703426 bytes (2.6 MB)

downloaded 2.6 MB

trying URL 'https://cran.rstudio.com/src/contrib/dtw_1.18-1.tar.gz'

Content type 'unknown' length 849256 bytes (829 KB) downloaded 829 KB trying URL 'https://cran.rstudio.com/src/contrib/diffusionMap_1.1-0.tar.gz'

Content type 'unknown' length 39130 bytes (38 KB)
downloaded 38 KB

trying URL 'https://cran.rstudio.com/src/contrib/Seurat_2.3.0.tar.gz'

Content type 'unknown' length 846065 bytes (826 KB)

downloaded 826 KB
* installing *source* package ‘igraph’ ...
** package ‘igraph’ successfully unpacked and MD5 sums checked
installing to /sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph/libs
** R
** demo
** inst
** preparing package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded

Error in dyn.load(file, DLLpath = DLLpath, ...) :    unable to load shared object '/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph/libs/igraph.so': libicui18n.so.58: cannot open shared object file: No such file or directory

Error: loading failed Execution halted ERROR: loading failed
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph’ Warning in install.packages :   installation of package ‘igraph’ had non-zero exit status

ERROR: dependency ‘proxy’ is not available for package ‘dtw’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/dtw’ Warning in install.packages :   installation of package ‘dtw’ had non-zero exit status

ERROR: dependency ‘igraph’ is not available for package ‘diffusionMap’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/diffusionMap’ Warning in install.packages :   installation of package ‘diffusionMap’ had non-zero exit status

ERROR: dependencies ‘igraph’, ‘dtw’, ‘diffusionMap’ are not available for package ‘Seurat’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/Seurat’ Warning in install.packages :   installation of package ‘Seurat’ had non-zero exit status

```

How could it be fixed? What to try?

Operating system is `CentOS Linux release 7.3.1611 (Core)`

> Update

I also tried development version of `Seurat` but it failed too: http://satijalab.org/seurat/devinstall.html

> Update

I tried installing `Seurat` using `conda install -c bioconda r-seurat` and it is not working. It is giving the error:

> -bash-4.2$ conda install -c bioconda r-seurat Fetching package metadata ............. Solving package specifications:
>
> PackageNotFoundError: Packages missing in current channels:

* r-seurat -\> r-cowplot
* r-seurat -\> r-diffusionmap
* r-seurat -\> r-dtw
* r-seurat -\> r-ggridges
* r-seurat -\> r-ica
* r-seurat -\> r-pbapply
* r-seurat -\> r-ranger
* r-seurat -\> r-rcppprogress
* r-seurat -\> r-sdmtools
* r-seurat -\> r-tclust

> We have searched for the packages in the following channels:

# Answer
> 2 votes
As @terdon mentioned in the comments, you need to `sudo yum install libicu` on the relevant nodes (or have your admin do that for you). Alternatively, seurat can be installed via conda, which means you don't need root access.

---
Tags: r, seurat
---


title: Protein expression units?

text:
# Question
Title: Protein expression units?
I'm starting for protein expression analysis.

In gene expression, we have TPM, RPKM etc. However, all I can see for protein expression is: "High", "Medium", "Low" and "Not detected".

> https://www.proteinatlas.org/ENSG00000140326-CDAN1/tissue

Q: Is there a unit for measuring protein expression like gene expression? I've seen people talking about PE (protein existence). Can we use it for protein expression?

# Answer
I would guess that the real unit for measuring protein expression would be pM (pico-moles) or a count of protein molecules.

I'm not sure we have any methods to calculate such quantities though.
> 1 votes

# Answer
The unit of measurement in protein expression is more or less "abundance". In mass spec this is derived from the area under the curve in the spectrum. In things like an ELISA this is derived from a color intensity.

Protein existence is a vague summarization of multiple data sources to try to quantify the likelihood of a protein existing. You can't use that for expression, though expression would influence the likelihood of existence.
> 0 votes

---
Tags: proteins
---


title: visualisation of genome alignment

text:
# Question
Title: visualisation of genome alignment
I was asked to check the synteny of some genes in two genome assemblies of the same species (PacBio & Illumina). I was given two scaffolds couples (so 4 scaffolds total) in the Illumina genome and asked if they were assembled together in the PacBio assembly.

So, I proceeded with a DNA-DNA alignment with minimap2 of the two Illumina scaffolds on the PacBio whole genome. But I am not sure on how to directly use the sam file for the question I was asked. So I tried to visualize this alignment using Mauve. It worked for the first couple, but it failed (blank page) for the second couple of scaffolds.

Surely, Mauve is not the best option out, so can you suggest an alternative? This is a chordata genome.

# Answer
After a bit of digging, finally found exactly what I was looking for. I'll share here in case someone else was looking for the same thing.

I used a standard Nucleotide BLAST (blastn), coupled with MUMmer. Short tutorial.

Assuming that `xn` & `yn` are the couple of scaffolds I want to check if they are assembled together in the Pacbio assembly (a multi-fasta file):

1) first, I blasted the pair of Illumina scaffolds (as separate queries) to the PacBio genome (as subject):

`blastn -query ../illumina_asm/scaffold_x1.fa -subject ../pacbio_asm/genome.fa > out_scaffold_x1; blastn -query ../illumina_asm/scaffold_y1.fa -subject ../pacbio_asm/genome.fa > out_scaffold_y1`

The candidates PacBio scaffolds shall appear in the top blastn hits in both `out_scaffold_x1` & `out_scaffold_y1`.

Let's assume we found one common scaffold only - corresponding to the fasta entry `pacbio_scaffold` in the PacBio genome.

2) extract the sequence of the target scaffold from the PacBio assembly, using samtools faidx:

`samtools faidx ../pacbio_asm/genome.fa; samtools faidx ../pacbio_asm/genome.fa 'pacbio_scaffold' > pacbio_scaffold.fa`

3) concatenate the Illumina scaffolds to analyze:

`cat ../illumina_asm/scaffold_x1.fa ../illumina_asm/scaffold_y1.fa > illumina_scaffolds.fa`

4) compare the Illumina scaffolds to the PacBio assembly using nucmer from MUMmer:

`nucmer --prefix=pacbio_scaffold pacbio_scaffold.fa illumina_scaffolds.fa`

this will generate a `pacbio_scaffold.delta` file, used in the next step

5) filter the delta file to keep only the most meaningful alignments file using delta-filter from MUMmer (avoiding confusing messy plots later):

`delta-filter -m pacbio_scaffold.delta > pacbio_scaffold.delta.m`

6) use mummerplot from MUMmer to achieve a graph:

`mummerplot --png pacbio_scaffold.delta.m -R pacbio_scaffold.fa -Q illumina_scaffolds.fa --prefix=pacbio_scaffold -large -layout`

On the x-axis, the PacBio scaffold they have in common. On the y-axis, the two Illumina scaffolds (separated by the horizontal grey abline). Colors correspond to the orientation.

Example of Illumina scaffolds linked in the PacBio assembly:
> 4 votes

# Answer
I think you could use circos for this. The configuration is kind of tedious but it gives nice results. You would basically need to retrieve alignments coordinates and parse them into this format.

```
scaf1 19720 31917 scaf2 28307 34227
scaf1 28307 34227 scaf2 19720 31917

```

Example figure:
> 3 votes

# Answer
I remember a few options for pairwise genome visualisation (which may be outdated BTW):
> 2 votes

---
Tags: sequence-alignment, software-recommendation, mauve
---


# Question
Title: groupGO parameters explanation (ont and level)
I need to use `groupGO` in `clusterProfiler` to find functional profile of a list of genes, and I am having trouble finding out what some parameters of the function mean and which I should select for my specific case.

The function is defined as:

```
groupGO(gene, OrgDb, keyType = "ENTREZID", ont = "CC", level = 2,
  readable = FALSE)

```

I have mouse genes dataset. What I should select for `ont` and `level` parameters? What does `ont` types: `MF`, `BP`, `CC` mean? What is `level`? Where I could find this info? It is definitely not in the `vignette`...

# Answer
As mentioned in the comments `MF` is `molecular function`, `BP` is `biological process` and `CC` is `cellular compartment`. These are the 3 domains of the ontologies. The `level` refers to the level in the ontology graph. In the example in that link, `pigmentation` would be level 1, `pigmentation during development` level 2, `regulation of pigmentation during development` is level 3 and so on. You can find more discussion of this in this post on biostars from a year ago.

Regarding the exact settings you should use, I don't personally find `CC` to be an informative domain. The first couple levels are typically too generic to be useful and they'll include too many genes anyway for you to likely pick up any changes. I'd start around level 3 and see how things go (in other programs you would get all of the levels and domains at once).
> 1 votes

---
Tags: rna-seq, groupgo, clusterprofiler
---


# Question
Title: clusterProfiler, groupGO: meaning of GeneRatio
I am running `groupGO` function from `clusterProfiler`:

```
ggo <- groupGO(gene     = head(gene, 200),
           OrgDb    = org.Mm.eg.db,
           keyType = "ENSEMBL",
           ont      = "BP",
           level    = 2,
           readable = TRUE)

```

And getting the following output:

```
> head(ggo)
                   ID            Description Count GeneRatio
GO:0000003 GO:0000003           reproduction     5     5/200
GO:0008152 GO:0008152      metabolic process    40    40/200
GO:0001906 GO:0001906           cell killing     0     0/200
GO:0002376 GO:0002376  immune system process    20    20/200
GO:0006791 GO:0006791     sulfur utilization     0     0/200
GO:0006794 GO:0006794 phosphorus utilization     0     0/200

```

As can be seen `cell killing`, `sulfur utilization` and `phosphorus utilization` have `zero` genes. Then, why are they here at all? Should not the `groupGO` show only `biological processes` the fed in genes are related to? This seems very weird to me. What is the meaning of the `GeneRatio` column?

# Answer
As you seem to be learning, the clusterProfiler package is a bit "special". It doesn't filter the results you're seeing for relevance, it returns everything. The `GeneRatio` is the number of genes in your input list associated with the given GO term (e.g., `5`) and then a `/` and the total number of input genes. This is generally not that useful. Since many things aren't documented, you may find the source code helpful.
> 2 votes

---
Tags: rna-seq, groupgo, clusterprofiler
---


# Question
Title: What are negative FSC-A values in flow cytometry data?
I did flow cytometry analysis of *E. coli* cells 7 hours post inoculation (1% overnight culture in 3ml).

I am analysing the FCS3.0 files using the flowCore package in R. I observe that FSC-A has negative values in some events, although the corresponding FSC-H values were all positive. These points comprised just a minor fraction of the data (~250/10000 in this case). The FSC voltage was set at 788V and the threshold flow rate was 5000-10000 events/sec.

What do these negative values mean?

> **FSC-A vs FSC-H**: plotted only negative FSC-A values

---

> **SSC-H vs FSC-A**: red points indicate negative FSC-A values

# Answer
This is probably a matter of calibration. If this is not the raw data (actual voltage curves), look at that.

As an area below 0 is not possible, it has to be a calibration issue, either in the machine, or in the software.

the -H plots are about the maximum height of the voltage spike in the side or forward light detector, while the -W plots are about the width (using some definition of flank-detection, so variation prone). The area is also computed in various ways: using something like (H\*W)/c (with c something like 2) or with an actual computation of the area of the spike, but again this has all kinds of predefinitions in it (what's the beginning/end of a spike, what's the background-level, ...) Some of these predefinitions you can alter in software, some are hardcoded or even part of the functionality of the hardware, so it might be easy to impossible to get rid of the negative values.

The FSC-H values corresponding to the FSC-A values are in your plot, and positive, so the most obvious solution (a bad algorithm that calculated Area as Height times Width times Magic and did not account for negative Heights (which depending on spike definition, you might get)) is out.
> 3 votes

---
Tags: flow-cytometry
---