title: How to paste RSIDs in CADD output

text:
# Question
Title: How to paste RSIDs in CADD output
I want to paste RSIDs in CADD output as CADD does not give RSID column in its output. For this purpose I am using bedtools intersect to compare two files and have RSID column in my CADD file.

This is the command I am using.

```
bedtools intersect -a CADD.bed -b Input.vcf| cut -f 3

```

But this command is giving me error

```
ERROR: file CADD.bed has non positional records, which are only valid for the groupBy tool.

```

Is there any way to solve this error or any other way to get RSID column in CADD output file.

# Answer
It would be useful if you show us some line of your `CADD.bed` and `Input.vcf` and how your desired output should look like.

`bedtools` is complaining about non positional records in `CADD.bed`, which means there are lines in this file where the first three columns aren't `CHROM`, `START`, `END`.

I guess you will have some header lines in the beginning that aren't introduced by an `#`. Just delete them (or add in `#` to the beginning of the line) and try again.

fin swimmer
> 1 votes

---
Tags: vcf, bedtools, variation
---


# Question
Title: Generating the reconstructed alignment from BAM
I have a (small) BAM file with CIGAR and MD fields.

Question 1: What tools exists in Python and/or R to reconstruct the alignment between the reference and the read in a BAM? Given that this is a very standard analysis, I feel that there should be R packages or Python libraries which have functions which do this task given a BAM input...

My goal was to create a \*tsv file with two columns, one with the reconstructed alignment and the other column with the corresponding sequence. This would be useful in my analysis. (I would either use R data.table, or pandas in Python.)

Here is read and reference:

```
CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
|||||||||||||||| |||||||||||||||||
CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG

```

Here is the desired output dataframe (with only the one example row). We could use an R data.table/data.frame, or a pandas DataFrame. Column1 are the alignments (i.e. the reconstruction of the alignment between the reference and the read based on the CIGAR and MD strings), Column 2 are the corresponding sequences. :

```
ReadAligned    Reference
CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG    CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
...    ....

```

Question 2: Using the CIGAR/MD tags, one could then create columns with the number of insertions and the number of deletions. This problem feels much more straightforward, e.g. the python library `cigar` or the cigar parsing functions of `pysam`

# Answer
If all you want are the number of insertions/deletions, you can use this to extract that info:

```
from cigar import Cigar
c = c=Cigar('20M5I10M3I5M')
num_insertions = len(list(filter(lambda x: x[1] == 'I', c.items())))
insertion_length = sum([x[0] for x in list(filter(lambda x: x[1] == 'I', c.items()))])

```

Deletions would be the same thing with `'D'` in place of `'I'`
> 5 votes

# Answer
I am answering question 1. However, as I am unfamiliar with python, I am providing a javascript solution only to show the basic logic. To try the code at the end, you need to install node.js and run it with `node this-script.js`.

> Given that this is a very standard analysis

I tend to believe this is a very rare analysis. I have only seen it done twice. Converting CIGAR+MD to padded alignment is complicated in coding, as you have to keep sequence, cigar and MD in sync all the time. The code below generates two strings in one go. It may be logically easier to generate the reference and the query strings separately.

In general, MD is ill defined. It is too difficult to work with. If possible, avoid this tag. You can do a lot of things with the combination of CIGAR and NM.

```
function md2pad(CIGAR, MD, SEQ) {
    var m, re_cigar = /(\d+)([MIDSHNX=])/g, re_MD = /(\d+)|(\^[A-Za-z]+)|([A-Za-z])/g;
    var cigar = [];
    while ((m = re_cigar.exec(CIGAR)) != null) // parse CIGAR into an array
        if (m[2] != 'H') // do nothing for hard clipping, as it has no effect
            cigar.push([parseInt(m[1]), m[2]]);
    var k = 0, sx = '', sy = ''; // k: k-th cigar operator; x: ref; y: query
    // cx/cy: start coordinate of cigar[k] on ref/query; mx/my: start of current MD op
    var cx = 0, cy = 0, mx = 0, my = 0;
    while ((m = re_MD.exec(MD)) != null) {
        if (m[2] != null) { // deletion from the reference
            var len = m[2].length - 1; // length of deletion
            sx += m[2].substr(1), sy += Array(len+1).join("-");
            mx += len, cx += len, ++k; // consume a deletion D from ref
        } else { // copy or mismatch
            var ml = m[1] != null? parseInt(m[1]) : 1; // length of this MD operation
            while (k < cigar.length && cigar[k][1] != 'D') {
                var cl = cigar[k][0], op = cigar[k][1];
                if (op == 'M' || op == 'X' || op == '=') {
                    if (my + ml < cy + cl) { // an MD ends in the middle of cigar M
                        if (ml > 0) {
                            if (m[3] != null) sx += m[3], sy += SEQ[my];
                            else sx += SEQ.substr(my, ml), sy += SEQ.substr(my, ml);
                        }
                        mx += ml, my += ml, ml = 0;
                        break;
                    } else {
                        var dl = cy + cl - my;
                        sx += SEQ.substr(my, dl), sy += SEQ.substr(my, dl);
                        cx += cl, cy += cl, ++k; // consume an M operator
                        mx += dl, my += dl, ml -= dl;
                    }
                } else if (op == 'I') {
                    sy += SEQ.substr(cy, cl), sx += Array(cl+1).join("-");
                    cy += cl, my += cl, ++k; // consume an insertion I to ref
                } else if (op == 'S') {
                    cy += cl, my += cl, ++k; // consume a soft-clipping S
                } else throw Error("inconsistent MD"); // not working with N
            }
            if (ml != 0) throw Error("inconsistent MD");
        }
    }
    if (cx != mx || cy != my) throw Error("inconsistent MD");
    return [sx, sy];
}

var cigar = '16025H3M2I55M1D108M306H';
var MD = '7T15C24T9^C0G2T6T16T3A1A7G21T17G6C1T6T5A4';
var seq = 'GATCACAGGCCTATCACCCTATTAATCACTCACGGGAGCTCTCCATGCATCTGGTATTTTTTCGGGGGGGGATGCACGCGATAGCATCGCGGGCCGCTGGAACCGGAGCACCCTATGTCGCAGGATCTGTCTTTGATTCCTACCTCATGCCATTATTAATCGCGCCTA';
var s = md2pad(cigar, MD, seq);
console.log(s[0]);
console.log(s[1]);

```
> 5 votes

# Answer
This program in C++ is able to perform the task https://github.com/mlafave/sam2pairwise
> 0 votes

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


# Question
Title: full visualisation of draft genomes alignment
I have two draft genomes (aprox. 500Mb each), one with ~10'000 scaffolds (`genome1.fasta`) and one with ~1'000 contigs (`genome2.fasta`). I would like to visualize their pairwise alignment, which a while ago I found it could easily be done by mummerplot in MUMmer.

In mummer, pairwise draft genomes alignment is supported:

> Many times it is necessary to align two genomes that have not yet been completed, or two genomes with multiple chromosomes. This can make things a little more complicated, since a separate alignment would have to be generated for each possible pairing of the sequences. However, both NUCmer and PROmer automate this process and accept multi-FastA inputs, thus simplifying the process of aligning two sets of contigs, scaffolds or chromosomes.

It is still ok if I do it aligning small portions (e.g. contig vs scaffold), as I could get here. However, if I do it all-vs-all, I end up with a "confetti"-plot (self-coined term) instead of a dot-plot :

I am not sure how to interpret this, and what I would like to see is more of a straight line as the one. Probably it's just a matter of filtering, so I'd like to ask for advice.

I used nucmer through delta-filter as follows:

```
# mummer version //4.0.0beta1

# nucmer
nucmer --prefix=ref_qry genome1.fasta genome2.fasta

# ref_qry.delta generated from previous step
delta-filter -r -q ref_qry.delta > ref_qry.filter

# plot
mummerplot ref_qry.filter -R genome1.fasta -Q genome2.fasta --filter --layout

```

Reasons for options `-r -q` in `delta-filter` step :

> To determine a one-to-one mapping of each reference and query, combine the options and use delta-filter -r -q.

I'll also accept other methods for visualisation (doesn't need to be withh MUMmer). Any input is greatly appreciated!

---

*EDIT*

I found some solutions which I reported as answers. I will accept anyone giving answer on how to tune MUMmer parameters to achieve the goal I describe in the original question.

# Answer
Found a solution, using D-Genies, worked great.

Some examples from their website:

Thanks to @user172818.
> 6 votes

# Answer
Also, Mashmap has an embedded visualisation tool generating dot-plots to visualize mappings, ideal for:

* highly similar species or to compare different assemblies of the same species.
* mapping genome assembly or long reads (PacBio/ONT) to reference genome

An example from their website:
> 1 votes

---
Tags: sequence-alignment, genome, software-recommendation, visualization
---


# Question
Title: Difference between genometools installed with anaconda vs apt-get?
I installed genometools with anaconda `conda install -c bioconda genometools`. Something is installed though apparently no binary called `gt`. If I type in gt Ubuntu suggests to install genometools via `apt`. What is the difference between both versions? Why is there no `gt` binary in the anaconda package?

# Answer
These appear to be two separate collections of tools. The conda version is a collection of python scripts whereas the apt version is based on a C library. So, as you might expect, the python version does not have a `gt` binary.

They also have separate GitHub accounts conda genometools vs apt genometools.
> 7 votes

---
Tags: bioconda
---


title: What is mate rescue in bwa mem?

text:
# Question
Title: What is mate rescue in bwa mem?
BWA mem has the `-S` and `-P` tags for *skip mate rescue* and *skip pairing; mate rescue performed unless -S also in use*.

What do these do? I presume `-P` aligns read pairs independently of each other. Is that correct?

And what does `-S` do?

# Answer
> 6 votes
Using `-SP` is equivalent to running `bwa mem` on each of the two mates as if they are single-end reads, but it formats the output as a proper paired-end output, i.e. with all pair-related flags added properly. Without `-SP`, by default `bwa mem` forces an alignment of a poorly aligned read if its mate is aligned somewhere. `-SP` turns off the forced alignment. We use `-SP` for mapping Hi-C reads.

# Answer
> 2 votes
If one read maps and the other doesn't, BWA attempts to rescue the other read by performing Smith-Waterman alignment with the unmapped mate, `-S` disables this mate rescue.

`-P` does indeed disables pairing, but mate rescues is still enabled - I think `-P` disables the setting of the proper pair flag. If you set `-P` and `-S`, BWA maps the paired reads essentially as single reads.

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


title: What are 2D reads in the Oxford MinIon?

text:
# Question
Title: What are 2D reads in the Oxford MinIon?
Reading the Mash: fast genome and metagenome distance estimation paper, I see the reference to "2d reads" notion under the minION context (a sequencing technique)
e.g

```
    In both cases Mash was able to
correctly differentiate these closely related species
(ANI ≈ 95 %) using 43,806 and 91,379 sequences collected
from single MinION R7.3 runs of B. anthracis
Ames and B. cereus ATCC 10987, respectively (combined
1D and 2D reads) <----------  HERE.

```

But couldn't find any reference in Wikipedia or online papers to what it is
Can someone shed some light?

# Answer
> 5 votes
Early MinION sequencing runs had forward and reverse DNA templates joined together by a hairpin adapter, so that the sequencer would read both strands from the same template. The consensus sequence that was generated by combining the base calls of opposite strands is referred to as a 2D read.

Due to a dispute with Pacific Biosciences, and for other reasons related to chemistry and basecalling, the equivalent mode in current sequencing runs is called 1D². There is no physical ligation of forward and reverse templates, but the kit and flow cell chemistry are modified to encourage the reverse strand to follow the forward strand very soon afterwards through the pore, with the joining and consensus generation carried out in software after the sequencing is done. In the output fast5 files, these 1D² reads can also be referred to as 2D.

---
Tags: phylogenetics, genome-sequencing
---


# Question
Title: how to filter a multi-VCF for sex-specific genotypes?
I have a single VCF file of variants for a population that includes males & females and I would like to pull out all variants that are sex-specific. For example, sites where all males are 0/1 and all females are 0/0. Is there a way to do this using bcftools or vcftools? I could enumerate each individual name, or ask for all individuals named \*.M and \*.F.

I first tried running smoove on the male and female populations separately then comparing the two VCFs, but I'm working with a draft genome with loads of contigs. As only contigs containing variants are output in each VCF, and these are different between the two files, vcftools couldn't handle the comparison.

Someone here https://github.com/samtools/bcftools/issues/118 managed to get it working ike this:

```
$bcftools view -f PASS -i 'GT~"0/0.*0/1.*0/1.*0/1.*0/1"' combined.vcf.gz

```

but the same syntax gives me no variants (yet I can see some that meet my criteria in IGV).

# Answer
> 2 votes
Found a great tool with SnpSift!

From the examples: I want to keep samples where the genotype for the first sample is homozygous variant and the genotype for the second sample is reference:

```
cat variants.vcf | java -jar SnpSift.jar filter "isHom( GEN[0] ) & isVariant( GEN[0] ) & isRef( GEN[1] )" > filtered.vcf

```

# Answer
> 1 votes
if you could provide seperate list of male and female samples names you can something like this:

```
$ bcftools view -e 'GT[@males.txt]!="0/1" && GT[@females.txt]!="0/0"' combined.vcf.gz

```

or if you don't want to take care of the exact genotype you can write it more generic:

```
$ bcftools view -e 'GT[@males.txt]!="alt" && GT[@females.txt]!="ref"' combined.vcf.gz

```

fin swimmer

# Answer
> 1 votes
Given an input file `input.vcf`, I would take this approach:

1. Create two text files: one with the male sample names (`males.txt`) and one with the female sample names (`females.txt`) (one sample per line)
2. Create a BCF (for speed) for all the samples in each list while dropping the `INFO` field so that `INFO/AC` is recalculated (I'm not sure dropping the `INFO` field is strictly necessary). Then only keep sites with a minor allele count greater than zero:

```
bcftools view --samples-file males.txt input.vcf -Ou | bcftools annotate -x INFO  -Ou | bcftools view --min-ac 1:minor -Ob -o male_variants.bcf
bcftools view --samples-file females.txt input.vcf -Ou | bcftools annotate -x INFO -Ou | bcftools view --min-ac 1:minor -Ob -o female_variants.bcf

```
3. Use `isec` to create intersects and complements of the two variant files and store them in `dir/*`:

```
bcftools isec male_variants.bcf female_variants.bcf -p dir -Oz

```

Double-check `dir/README.txt`, but the files you seek should be named `dir/0000.vcf.gz` and `dir/0001.vcf.gz`

**Disclaimer:** This code is untested.

# Answer
> 1 votes
using vcffilterjdk http://lindenb.github.io/jvarkit/VcfFilterJdk.html

> or ask for all individuals named \*.M and \*.F.

```
java -jar dist/vcffilterjdk.jar -e 'return variant.getGenotypes().stream().allMatch(G->(G.getSampleName().endsWith("M") && G.isHet()) || (G.getSampleName().endsWith("F") && G.isHomRef())); ' input.vcf

```

`variant.getGenotypes().stream()` get a stream of genotypes 'G' for the current variant.

`allMatch` all genotypes must match the statements below:

`G.getSampleName().endsWith("M") && G.isHet()` sample 's name ends with 'M' and is heterozygous

OR

`G.getSampleName().endsWith("F") && G.isHet()` sample 's name ends with 'F' and is Hom-ref

---
Tags: vcf
---


# Question
Title: How should I address batch effects in my experiment?
Let's say I have an RNA-Seq experiment, where I'm interested in the significantly differentiated genes between pre-treatment and post-treatment conditions. "rep" == biological replicate.

```
Sample PreA (3 reps **all** in 2017) vs Sample PostA (3 reps **all** in 2018)
Sample PreB (3 reps **all** in 2017) vs Sample PostB (3 reps **all** in 2018)
Sample PreC (3 reps **all** in 2017) vs Sample PostC (3 reps **all** in 2018)

```

I also have a sample:

```
Sample PreD (3 reps all in 2017) and Sample PreD (3 reps again in 2018)

```

It is safe to assume the samples in sequenced in 2017 were done differently to 2018.

**Q1:** Is there problem with batch effects in my study? When I do a test, there should't be any way to split variation from batch effects and biological variation?

**Q2:** If Q1 is correct. How should I reduce batch effects?

* Can I add a batch factor to my model in DESeq2/edgeR?
* Should I use my `PreD` control samples to derive normalisation factor as in RUV? (assume PreD as a negative control)
* Should I use `sva` package?

**EDIT**

Devon suggested for adding a batch variable, but how? That's how I'm doing (ignoring batch effects):

```
dds <- DESeqDataSetFromTximport(txi, data, ~Sample)

```

`data` holds a data frame with a `Sample` column. `Sample` would be "PreA" and "PostA" in the first test. Very simple.

I couldn't add `Date` to it, like:

```
dds <- DESeqDataSetFromTximport(txi, data, ~Date + Sample)

```

because the code crashed complaining of the 2017/2018 problem discussed above. I'm not sure how I should add `PreD` to the design.

# Answer
1. Yes, though thankfully Your 2018 PreD samples will help you resolve this.
2. Simply add the batch effect to the design (`~Batch + Treatment`) and DESeq2 (or edgeR or Limma) will handle this for you.

You do not need SVA or RUV, thankfully, since you quite cleverly sequenced one group in both batches.

To clarify, your coldata will be something like:

```
Group    Time    Batch
A        Pre     A
A        Pre     A
A        Pre     A
A        Post    B
A        Post    B
A        Post    B
...
D        Pre     A
D        Pre     A
D        Pre     A
D        Pre     B
D        Pre     B
D        Pre     B

```

And the design is actually `~Batch + Group + Time`. You can probably add an interaction term in if you then modify the model matrix manually to ensure it's full rank.
> 4 votes

---
Tags: rna-seq, bioconductor, batch-effects
---


# Question
Title: What is a simple command line tool for doing Needleman-Wunsch pair-wise alignment on the command line
I have two DNA strings:

> GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC

and

> AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

I want a tool that allows me to do something like this on the command line:

```
$ aligner  GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

```

and receive an ASCII visualization of a pairwise alignment. Something like this would work:

```
GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---------
---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

```

I prefer a tool that is packaged on bioconda. A tool that also does protein and RNA sequences is even more preferable.

# Answer
You are looking for the `needle` program from the `EMBOSS` suite. Available in bioconda.

http://emboss.sourceforge.net/apps/release/6.6/emboss/apps/needle.html

To read sequences from the commandline, you need so specify the format as `asis`. To get the output on the screen, you'll need `-stdout` and to use the default alignment parameters (gap penalty 10, extend\_penalty 0.5) you'll need `-auto`. Thus your query above would be:

```
$ needle -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG -auto -stdout
########################################
# Program: needle
# Rundate: Tue 27 Nov 2018 10:48:50
# Commandline: needle
#    -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC
#    -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG
#    -auto
#    -stdout
# Align_format: srspair
# Report_file: stdout
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: asis
# 2: asis
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 56
# Identity:      37/56 (66.1%)
# Similarity:    37/56 (66.1%)
# Gaps:          18/56 (32.1%)
# Score: 181.0
#
#
#=======================================

asis               1 GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---     47
                              |||||||||||||||||||||||||||||||||||||.
asis               1 ---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAG     41

asis              47 ------     47

asis              42 AGGAGG     47

#---------------------------------------
#---------------------------------------

```
> 11 votes

# Answer
Install biopython with conda and use the following script:

```
#!/usr/bin/env python
import sys
from Bio import Align
aligner = Align.PairwiseAligner()
aligner.mode = "local"
alignments = aligner.align(sys.argv[1], sys.argv[2])
print(alignments[0])

```

The output is then:

```
.GGA-GGAGGGAG--AAGGAGGGAGGGAAGA-GGAGGGAG--AAGGAGGGAGGC
.|-|-||||||||--|||-|||-|||||-||-||||||||--|||-|||-|||.
AG-AAGGAGGGAGGGAAG-AGG-AGGGA-GAAGGAGGGAGGGAAG-AGG-AGG.

```

There are a LOT of alignments with that same score and you can tweak the details of the mode by changing `aligner.mode`.
> 9 votes

# Answer
I'd just use something like `t_coffee` and write a little bash function to run it on the strings you give:

```
aligner(){
  tmp=$(mktemp);
  args=("$@")
  for ((i=0; i<$#; i++)); do
    printf '>%s\n%s\n' $i ${args[i]} >> "$tmp"
  done
  t_coffee "$tmp"  2>/dev/null | /bin/grep -Pv '^CLUSTAL|^\s*$'
  rm "$tmp"

}

```

If you add that to your shell's initialization file (e.g. `~/.bashrc` if using bash on Linux, `~/.profile` if using bash on macOS), you can run `align2seqs` with as many sequences as you want to align:

```
$ aligner 'GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC' 'AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG'
0               GG--AGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC
1               nAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGG-GAAGAGGAGG
                    *********      *******  *********  * ** *  *
$ aligner 'GAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC' 'GGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG' GAGGGAGGGAAGAGGAGGGAGAAGGATCCGAGGAAGAGGAGG AGAAGGAGGGAGGGAAGAGGCCCGAGATATATAAGGATCCGAGGAA
0               GA-GG----GAGAAGGAGGGAGGGAAGAGGAGGGAGAAGG--AGGGAGGC
1               GGAGG----GAGGGAAGAGGAGGGAGAAGGAGGGAG-GGA--AGAGGAGG
2               GA-GG----GAGGGAAGAGGAGGGAGAAGGATCCGA-GGA--AGAGGAGG
3               AGAAGGAGGGAGGGAAGAGGCCCGAGATATAT--AA-GGATCCGAG-GAA
                    *    ***      **   **     *       *    * *

```
> 4 votes

# Answer
Here's another option: https://github.com/noporpoise/seq-align; I've used it in the past as it is small and has no dependencies (just `make` and you're set).
> 3 votes

---
Tags: sequence-alignment
---


# Question
Title: How to separate a microbial strain of MTB resistant and susceptible drugwise?
Im trying to classify my list into three type: XDR - A strain must be resistant to:

```
MDR + CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN + AMYCACIN/CAPREOMYCIN/KANAMYCIN-> means ISONIAZID + RIFAMPICIN + any of these or two or three or all (CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN) + any of these or two or three or all (AMYCACIN/CAPREOMYCIN/KANAMYCIN)

```

PRE-XDR - A strain must be resistant to:

```
MDR + CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN -> means ISONIAZID + RIFAMPICIN + any of these or two or three or all (CIPROFLOXACIN/MOXIFLOXACIN/OFLOXACIN/LEVOFLOXACIN)

```

MDR - A strain must be resistant to:

```
ISONIAZID + RIFAMPICIN
and if it is resistant to any one/two/three of these or all of these in
addition to ISONIAZID + RIFAMPICIN than again it is MDR:-
PYRAZINAMIDE, ETHIONAMIDE, STREPTOMYCIN, ETHAMBUTOL

```

SUSCEPTIBLE - A strain must be Susceptible to all drugs

I tried the code as follows for MDR:

```
tail -n +2 drug.txt |\
awk '{a[$1] = a[$1]","$2};END{for(i in a)print i"\t"a[i]}'|\
grep -P "\t,ISONIAZID,RIFAMPICIN$"

```

It gave me MDR starins correctly but than i had to first remove all Susceptible strains and applied it to Resistant one only.

But when i applied the code for Pre-XDR:

```
tail -n +2 drug.txt |\
awk '{a[$1] = a[$1]","$2};END{for(i in a)print i"\t"a[i]}'|\
grep -P
"\t,ISONIAZID,RIFAMPICIN|CIPROFLOXACIN|MOXIFLOXACIN|OFLOXACIN|LEVOFLOXACIN"

```

It is giving me the result as below:

```
SRR671821       ,CAPREOMYCIN,ETHAMBUTOL,ETHIONAMIDE,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN,STREPTOMYCIN
SRR671863       ,CAPREOMYCIN,ETHAMBUTOL,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN
SRR671822       ,ETHAMBUTOL,ISONIAZID,OFLOXACIN,RIFAMPICIN,STREPTOMYCIN
SRR671864       ,CAPREOMYCIN,ETHAMBUTOL,ISONIAZID,KANAMYCIN,OFLOXACIN,RIFAMPICIN

```

Which is not correct for Pre-XDR.

So I'm looking for a the drugs if match than look for resistant and susceptible if these criteria match with any of the mention type (MDR,XDR...) it put those under that in a new file. Can anyone help? my file druglist.sxls

MDR eg from the file attached is:

```
ERR038736   ISONIAZID   RESISTANT
ERR038736   RIFAMPICIN  RESISTANT
ERR038737   ISONIAZID   RESISTANT
ERR038737   RIFAMPICIN  RESISTANT
SRR1270146  AMIKACIN    SUSCEPTIBLE
SRR1270146  CAPREOMYCIN SUSCEPTIBLE
SRR1270146  CYCLOSERINE SUSCEPTIBLE
SRR1270146  ETHAMBUTOL  RESISTANT
SRR1270146  ETHIONAMIDE SUSCEPTIBLE
SRR1270146  ISONIAZID   RESISTANT
SRR1270146  KANAMYCIN   SUSCEPTIBLE
SRR1270146  PARA-AMINOSALICYLIC ACID    SUSCEPTIBLE
SRR1270146  PYRAZINAMIDE    RESISTANT
SRR1270146  RIFAMPICIN  RESISTANT
SRR1270146  STREPTOMYCIN    RESISTANT

```

# Answer
> 2 votes
Since this involves combining multiple conditions, it is easier (and certainly more readable if you ever need to come back to it later) to write a script to do this. For example, in perl:

```
#!/usr/bin/perl

use strict;
use warnings;

## This will hold the information for each strain
my %strains;

## These are the functions that will test each strain
sub isMDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if ($thisStrain{"ISONIAZID"} && $thisStrain{"RIFAMPICIN"}){
    foreach my $drug (keys(%thisStrain)) {
      next if $drug eq "ISONIAZID";
      next if $drug eq "RIFAMPICIN";
      next if $drug eq "PYRAZINAMID";
      next if $drug eq "ETHAMBUTOL";
      next if $drug eq "STREPTOMYCIN";
      ## discard this strain if it is resistant to any other drug
      if ($thisStrain{$drug}) {
        return 0;
      }
    }

    return 1;
  }
  return 0;
}

sub isPreXDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if (
  ($thisStrain{"ISONIAZID"} && $thisStrain{"RIFAMPICIN"}) &&
  (
   $thisStrain{"CIPROFLOXACIN"}  ||
       $thisStrain{"MOXIFLOXACIN"} ||
   $thisStrain{"OFLOXACIN"}  ||
       $thisStrain{"LEVOFLOXACIN"}
      )
     ) {
    return 1;
  }
  return 0;
}

sub isXDR{
  my $_strain = shift;
  my %thisStrain = %{$strains{$_strain}};
  if (
  isPreXDR($_strain) &&
      (
       $thisStrain{"AMYCACIN"} ||
   $thisStrain{"CAPREOMYCIN"} ||
       $thisStrain{"KANAMYCIN"}
      )
     ){
    return 1;
  }
  return 0;
}

sub isSusceptible{
  my $_strain = shift;
  foreach my $drug (keys(%{$strains{$_strain}})) {
if ($strains{$_strain}{$drug}) {
      return 0
    }
  }
  return 1
}

## Iterate over each line in the input file.
## For each iteration, the special variable $_
# will hold the contents of the line.
while (<>) {
  ## Skip first line
  next if $. == 1;
  ## remove Windows-style line endings
  s/\r//;
  ## remove trailing newline
  chomp;
  ## split the line on commas
  my ($strain, $drug, $status) = split(/,/);
  if ($status eq "SUSCEPTIBLE"){
$status = 0
  }
  else {
$status = 1;
  }
  $strains{$strain}{$drug}=$status;
}

foreach my $strain (keys(%strains)) {
  print "$strain\t";
  if (isMDR($strain)) {
print "MDR\n";
  }
  elsif (isXDR($strain)) {
print "XDR\n"
  }
  elsif (isPreXDR($strain)) {
print "PRE-XDR\n"
  }
  elsif (isSusceptible($strain)) {
    print "SUSCEPTIBLE\n";
  }
  else {
    print "OTHER\n";
  }
}

```

If you save that script as `foo.pl`, you can run it on your file with:

```
perl foo.pl drug.csv

```

Selected output from running this on the `drug.csv` file you posted is:

```
$ foo.pl drug.csv | grep -E 'ERR038288|ERR038736|ERR038751|SRS553990|ERR133879'
ERR038751   OTHER
ERR038288   SUSCEPTIBLE
SRS553990   XDR
ERR038736   MDR
ERR133879   PRE-XDR

```

---

Various details about this question were added in chat: https://chat.stackexchange.com/transcript/message/47821646#47821646

---
Tags: python, perl, text-processing
---


title: What is Allelic Imbalance

text:
# Question
Title: What is Allelic Imbalance
Can anyone help me explain what allelic imbalance is, hopefully shortly? Surprisingly, we cannot find any introduction online.

# Answer
It might be used differently in different contexts, but generally speaking, in my world - allelic imbalance is when there's a difference in the level of gene expression from different alleles, usually through genetic (e.g. a variant in a promoter) or epigenetic mechanisms (e.g. one copy silenced, as in imprinted regions).

There are a few blog posts here with some more information, and plenty of papers (google "Allelic imbalance") if you want to dig deeper.
> 5 votes

---
Tags: phylogenetics
---


# Question
Title: Creating a computational table to show the number of available genes between all possible pairs of combinations in eight gene lists
I am looking to create a table (on RStudio) that will show the overlap between 8 different gene lists wherever possible. I have saved the 8 lists as .txt files all in one folder and set that particular folder as my working directory. Any idea of code that will work for this?

If anyone knows what I can alter to make the current code work or suggest an alternative way I would be very happy to hear any suggestions!

```
data.file1 <- "GSE108363_BCGdown_V.txt"
data.file2 <- "GSE108363_BCGdown_D.txt"
data.file3 <- "GSE108363_BCGup_V.txt"
data.file4 <- "GSE108363_BCGup_D.txt"
data.file5 <- "GSE108363_MTBdown_V.txt"
data.file6 <- "GSE108363_MTBdown_D.txt"
data.file7 <- "GSE108363_MTBup_V.txt"
data.file8 <- "GSE108363_MTBup_D.txt"

genevect1 <- scan(data.file1, what=character(), sep="\n")
genevect2 <- scan(data.file2, what=character(), sep="\n")
genevect3 <- scan(data.file3, what=character(), sep="\n")
genevect4 <- scan(data.file4, what=character(), sep="\n")
genevect5 <- scan(data.file5, what=character(), sep="\n")
genevect6 <- scan(data.file6, what=character(), sep="\n")
genevect7 <- scan(data.file7, what=character(), sep="\n")
genevect8 <- scan(data.file8, what=character(), sep="\n")

filelist <- list(data.file1, data.file2, data.file3, data.file4, data.file5, data.file6, data.file7, data.file8)
#-------------------------------------------------------------------------------------------------------
# "Lapply" is basically making a "for loop".
#-------------------------------------------------------------------------------------------------------
all_gene_vectors <- lapply(filelist, scan, what=character(), sep="\n")

#-------------------------------------------------------------------------------------------------------
# Then make the intersection of *all* genesets
#-------------------------------------------------------------------------------------------------------
final_inter <- all_gene_vectors[[1]]
for (next_genevect in all_gene_vectors[2:length(all_gene_vectors)]) {
  final_inter <- intersect(final_inter, next_genevect)
}
show(final_inter)

#-------------------------------------------------------------------------------------------------------
# Or compute all pairwise intersections (untested)
#-------------------------------------------------------------------------------------------------------
pairwise_inters <- list()
#-------------------------------------------------------------------------------------------------------
# Iterate over all possible pairs i,j (i ≠ j)
#-------------------------------------------------------------------------------------------------------
N <- length(all_gene_vectors)
for (i in 1:(N-1)) {
  ith_inters <- list()
  for (j in (i+1):N) {
    ith_inters[[ names(all_gene_vectors)[j ]]] <- intersect(all_gene_vectors[[i]], all_gene_vectors[[j]])
  }
  pairwise_inters[[ names(all_gene_vectors)[i] ]] <- ith_inters
}

#-------------------------------------------------------------------------------------------------------
#Generate table in wide format:
#-------------------------------------------------------------------------------------------------------
#Concatenate gene names of all tables to a separate file
#Extract unique values in this file
#Perform outer join between the separate file and each table

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

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

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

set.seed(11)
BCG_validation_Up <- sample(letters[1:429], 20)
BCG_discovery_Up <- sample(letters[1:250], 20)
MTB_validation_Up <- sample(letters[1:286], 20)
MTB_discovery_Up <- sample(letters[1:128], 20)
BCG_validation_Down <- sample(letters[1:267], 20)
BCG_discovery_Down  <- sample(letters[1:350], 20)
MTB_validation_Down <- sample(letters[1:244], 20)
MTB_discovery_Down <- sample(letters[1:86], 20)

cross_table <- matrix(, nrow = 4, ncol = 4)
rownames(cross_table) <- c("BCG_validation_Up", "BCG_discovery_Up", "MTB_validation_Up", "MTB_discovery_Up")
colnames(cross_table) <- c("BCG_validation_Down", "BCG_discovery_Down", "MTB_validation_Up", "MTB_discovery_Up")

for (i in 1:4){
  for(j in 1:4){
    cross_table[i,j] <- length(intersect(get(paste0("",i,"_Up")),(get(paste0("",j,"_Down")))))
  }
}

cross_table

```

# Answer
> 1 votes
To make a confusion matrix of 4 lists against 4 other lists, you need to get the length of all intersecting genes for each pair. Here an example with letters instead of gene names.

```
set.seed(11)
list1_up <- sample(letters[1:26], 20)
list1_down <- sample(letters[1:26], 20)
list2_up <- sample(letters[1:26], 20)
list2_down <- sample(letters[1:26], 20)
list3_up <- sample(letters[1:26], 20)
list3_down <- sample(letters[1:26], 20)
list4_up <- sample(letters[1:26], 20)
list4_down <- sample(letters[1:26], 20)

cross_table <- matrix(, nrow = 4, ncol = 4)
rownames(cross_table) <- c("list1_up", "list2_up", "list3_up", "list4_up")
colnames(cross_table) <- c("list1_down", "list2_down", "list3_down", "list4_down")

for (i in 1:4){
  for(j in 1:4){
    cross_table[i,j] <- length(intersect(get(paste0("list",i,"_up")),(get(paste0("list",j,"_down")))))
    }
}

cross_table
         list1_down list2_down list3_down list4_down
list1_up         15         14         15         15
list2_up         18         15         15         15
list3_up         17         15         14         16
list4_up         17         15         15         16

```

---
Tags: r, rstudio
---


# Question
Title: Score predefined ChIP-seq peaks with MACS2 or equivalent
I have performed peak calling on a number of separate ChIP-seq experiments and would like to harmonise the peaks from each of the experiments in order to convert my data into a matrix for further analysis. I would like to do this by taking the union of the called peaks (available in BED / narrowPeak / MACS2 format) and retrieving the enrichment score / q-value for all of the experiments on each peak. Is this possible without writing my own extension to MACS2?

# Answer
You could try setting a high p-value threshold when you're calling peaks to retain the "non-significant" values. Something like:

```
macs2 callpeak -t <ChIP>.bam -c <Control>.bam -f BED -g hs -n test -B -p 0.7

```

You probably don't want to set the threshold to 1, since you'll likely want to avoid some noise.

But you can then use these "non-significant" values in combination with a `bedtools map` command, to map them to your predefined union of "actual" peaks.
> 2 votes

---
Tags: chip-seq, peak-calling, macs2
---


title: Can't install newest Blast from Conda

text:
# Question
Title: Can't install newest Blast from Conda
I have run into an odd situation, where Anaconda repo page clearly states that the newest version of Blast in Bioconda channel is 2.7.1, however the version Conda wants to install (command `conda install -c bioconda blast`) is 2.6.0.

Even specifying the version with command `conda install -c bioconda blast=2.7.1` does not work and Conda returns a `PackagesNotFoundError`.

By observing the webpage mentioned above, one can see that there is a "v" before the version number, however adding this letter to the command with specified version did not help.

To make things even more odd, `conda search blast` returns a list of Blast versions available and there is a 2.7.1 version available in Bioconda channel.

I have added the Bioconda channel into my Conda.

My question is: How can I install the 2.7.1 (by the time of writing this question the newest version available in Anaconda repo) through Conda?

Edit: I have Ubuntu Server 18.04 LTS and Conda 4.5.11.

# Answer
I was able to reproduce your problem, even in a blank environment.

```
$ conda create -n test
$ conda activate test
$ conda install bioconda::blast=2.7.1
Solving environment: failed

PackagesNotFoundError: The following packages are not available from current channels:

  - bioconda::blast=2.7.1
  - boost=1.64
  - bioconda::blast=2.7.1
  - gnutls
  - bioconda::blast=2.7.1
  - nettle[version='3.3|3.3.*']
  - bioconda::blast=2.7.1

```

Note that `boost`, `gnutls`, and `nettle` are listed as requirements here.

If you do a search for those packages, you'll see that `boost` is available in a default channel, whereas `gnutls` and `nettle` are only in the `conda-forge` channel, which is not a default. Adding those packages to the install command resolved this issue, for me.

```
$ conda install bioconda::blast=2.7.1 conda-forge::gnutls conda-forge::nettle --dry-run
Solving environment: done

## Package Plan ##

environment location: /Users/hawleyj/anaconda3/envs/test

added / updated specs:
    - bioconda::blast=2.7.1
    - conda-forge::gnutls
    - conda-forge::nettle
...

```

As mentioned in the comment below, to avoid this issue entirely, before using Bioconda packages make sure you run the proper channel configuration first.

```
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge

```
> 9 votes

---
Tags: blast, bioconda
---


# Question
Title: Why the t-test for a specific gene shows different value compared to differential analysis?
I have RNA-Seq data for LUNG cancer. 370 tumor and 50 Normal. For differential analysis initially I did some filtering and kept approx. 19k genes for further analysis. I used edgeR.

With a FC \> or \< 2 and FDR \< 0.05
`tr <- glmTreat(fit, contrast=contrast.matrix, lfc = log2(2))`
I got approx 3000 differentially expressed genes.

Among the 3000 differentially expressed genes I'm interested in a gene `OIT3`.

```
GeneSymbol  logFC    unshrunk.logFC   logCPM          PValue      FDR
  OIT3  2.156487448  2.171599487    -0.437234244    0.000701037 0.005304312

```

I see the above gene `OIT3` is Upregulated. Why I see logCPM with negative value?

And I plotted the logCPM expression of the above gene and did a t-test which gave pvalue 0.1 which means not significant.

1) I see the gene is differentially expressed with Differential analysis using edgeR. Why I see logCPM with negative value?

2) With t-test why it is not significant?

# Answer
1. The $log(CPM)$ of any low-moderately expressed gene will be negative. There is nothing unexpected there.
2. Your statistics are inappropriate for a variety of reasons. Firstly, a CPM is not a robust value that's comparable between samples (this is why CPMs aren't used for statistics). `edgeR` performs more appropriate normalization and incorporates that into its model. Secondly, though you have a fair number of replicates, edgeR is pooling information from all of the genes to better estimate the variance.

I encourage you to read the edgeR paper and likely that of Limma too.
> 3 votes

# Answer
There is no reason your t-test should reproduce edgeR. In fact, edgeR exists because t-test is inappropriate.

edgeR does the tests by pooling information from all genes, because with the low number of replicates your t-test doesn't give sufficient statistical power.

What you need to do is: *check visually your gene and make a decision yourself based on graphical results and edgeR results if your gene shows biologically significant*.
> 4 votes

---
Tags: r, rna-seq, rna, cancer, gene-expression
---


# Question
Title: What command/invocation is used to generate NCBI 16SMicrobial blastdb
I'm looking for the exact invocation used to generate the 16SMicrobial database that you can download from here:

https://ftp.ncbi.nlm.nih.gov/blast/db/

I'm hoping to create the same type of blastdb with the same type of metadata with custom sequences.

Platform isn't an issue, but let's say on ubuntu 14.04 or 16.04.

I would like to replicate the creation of the database as closely as possible. The most important feature is the taxonomic information as can be seen here:

> The databases on the FTP site contain taxonomic information for each sequence, include the identifier indices for lookups, and can be up to four times smaller than the FASTA. The original FASTA can be generated from the BLAST database using blastdbcmd

Creation of a blastdb using makeblastdb from a set of a fasta sequences is not an issue and can be achieved via:

```
makeblastdb -in <your_file.fasta> -dbtype nucl -out <database_name>

```

My question is specifically about the invocation NCBI uses to add the metadata that is present in the NCBI's 16SMicrobial blast database as I am keen to make sure I have replicated the process as closely as possible.

# Answer
> 1 votes
Running `blastdbcmd -db foo -info` provides a little information but I haven't seen anything which will report exactly how a blastdb was created.

A blastdb readme suggests that only the `parse_seqids` option has been added to the standard parameters.

> For those from NCBI, the following makeblastdb commands are recommended:
>
> * For nucleotide fasta file: makeblastdb -in input\_db -dbtype nucl -parse\_seqids
> * For protein fasta file: makeblastdb -in input\_db -dbtype prot -parse\_seqids

The `-parse_seqids` requires NCBI fasta headers in the correct format (spec can be found here).

# Answer
> 0 votes
I don't know what specific properties you are considering. As far as I know, that's just a normal blast database, like any other. These are produced using the `formatdb` command which is part of the `ncbi-blast` software. You haven't specified what system you are using so I can't help you find and install it, but the command will be something like:

```
formatdb -p F -i yourSeqs.fasta -n yourDBName

```
* `formatdb` : the command
* `-p F` : protein: false. This is needed to make databases of nucleotide sequences.
* `-i yourSeqs.fasta` : your input file with all your sequences in fasta format.
* `-n yourDBName` : the name (choose whatever you like) for your database.

You can then use `yourDBName` as a database to blast against. The above assumes yuou are

# Answer
> 0 votes
I am not sure whether you are looking for how to produce taxonomy file or rep\_set file (qiime reference file types). But you can produce fasta file (rep\_Set) like that:

```
$blastdbcmd -db 16SMicrobial/16SMicrobial -outfmt '%f' -entry 'all' > 16SMicrobial.fasta

```

Then you can create QIIME-compatible id\_to\_tax maps with this code: https://github.com/mtruglio/QIIME\_utilities/blob/master/id\_to\_tax\_mapmaker.sh

Hope this helps.

---
Tags: blast
---


# Question
Title: How to create a graph of multi column gene dataset?
well, I'm new to Bio. what I'm going to do is create a network/graph of a dataset related to the genes I think. here's the sample of the nodes dataset with CSV format named lung-nodes :

and this is a sample of edges I think with trrust name :

these data are stored in two separate file

what I must do is create a graph(network) of these nodes and edges by Python and a visualization library like networkX, but I do not how. I must make a relation between these to the dataset by creating a graph of them can you help me out?

# Answer
As far as I understand, it seems that the first table you provided contains biological/experimental information related to your genes, i.e. the nodes of the network you will build. The second table is the one you need to build the network: each line is an edge. You can use it as input and build the network using the IGraph library which works for both Python (as you requested) or R (as you suggested as an alternative).

Most probably, but I can't say for sure, the first table was just a dataset from which you (or someone for you) recovered a list of nodes that where then submitted to a database, i.e. String, commonly used to recover information related to the interactions between your nodes, i.e. genes. Once you have this info (in the second table) you can build the network!

edit:

if you're not familiar with IGraph, you can have a look at this question describing how to obtain an IGraph graph from a csv file or this other one which uses a different csv file format to build a graph
> 1 votes

---
Tags: gene, networks, visualization
---


title: Subset data frame based on ID

text:
# Question
Title: Subset data frame based on ID
I tried to subset the dataframe in R since there are empty space after each line in the sequence which shows 'NA' when i read the files into R ,as a result of which when i subset it just takes the sequence Protein\_IDS and sequence line by line not as whole

Im looking for this Protein\_ID in one column and respective sequence in other

```
Proteins_IDs                Sequence
  <chr>                       <chr>
1 fig|6666666.151743.peg.5284 MERQAVAVVFHDWRRAIQQRNKTAAADRGVDPVGVGHNSKIESLQPRSSVCRASLIQVPG
2 NA                          RLLNGEIFYSLAEAKVLIEAWRRHYNTVRPHSSLGYRPPAPETATPPYPASGSASLHLRP
3 NA                          DMAAMGLIH
4 NA                          NA
5 fig|6666666.151743.peg.5285 MIDRQTAQVRRWSFVLGSSIAAMAAASATPAWAQCSPDPTVADSVTNCTGTDSDGLRVTT
6 NA                          DGSMVEVAPSASVSGAGAAAITVDIAPVSGSYYARKATINVGGSVDGGTQSGISVLSGSV

```

So i have made another file with list of protein ID which i want to subset from this file..I tried removing empty space from the file but couldn;t as it the sequence should come in one line but instead there is line break

data files

I tried this `SEQ <- inner_join(seq,protein_id)` whose output is this its like each ID is associated with multiple sequence

```
fig|6666666.151743.peg.5044 MTKNTIHRTILVCASLIAMPAMAQGTTEAGQEGASLADIVVTAQRRSENLQNVPISVTAI
2 fig|6666666.151743.peg.5063 MTVEQALARMLEGTGLEAQRLSPTAYRLVQAAPPASRRDPVRPAPPLPPEGSVDIIVTGT
3 fig|6666666.151743.peg.5135 MLRVADGYVEEAADDAIVPSPAVPLSVPSAPVRFEAAPARYGARRGPNVPAILVIVLLHA

```

# Answer
> 5 votes
I want to advise you to use FASTA format, instead of obscure unknown formats. In R, there are many packages that can handle FASTA format and it is therefore more convenient. I will show you how to make a FASTA file from your data, and then I will show an example of how to import such a FASTA file and sub select your proteins of interest.

```
input <- scan("seq2.txt", what = character(), sep = "\t")
length(input)
new <- gsub("f", ">f", input[-1:-2])
write.table(new, "seq2.fasta", quote = F, row.names = F, col.names = F)

library(seqinr)
fasta.file <- read.fasta("seq2.fasta", seqtype = "AA", as.string = TRUE)
prot.id <- scan("protein_id.txt", what = character(), sep = "\t")
selection <- fasta.file[names(fasta.file) %in% prot.id]
length(selection)
[1] 87

head(selection, 2)
$`fig|6666666.151743.peg.5044`
[1] "MTKNTIHRTILVCASLIAMPAMAQGTTEAGQEGASLADIVVTAQRRSENLQNVPISVTAITSEQAIKLGITQPDVIAQMVPGLQLNHSAWGLTPFLRGVGTNSVAIGNEPAVALFVDDVYLPTGSAGVFQLNSISSIEVLKGPQGTLFGRNATGGVIHIRTRDPSFDPSASATASYGNYDTISGQLYGSTGLTDTVAINFALQGIDQRKGWGVNGTTGNEANSMSGWGGRAKLLWNMGEETSLLVAGSYDRRSSDQGQPGRSIAGTFTRDGYSADAAGLGFYDVMSNLDTVQKSTTRSLSGKLTHDFGSARLVLISAFSKTVVDPVVSDSDGSPSPFFQTREFAMGGRTFTQEVQLLSPEGSAVKWILGAFYLHDRSFMDATYEGQAFGTSKSHTDAVQRTNSIAGFGQVTTEILPRTNLTLGLRYTSDRRSLDAVARLLPATGAAVTFGPFEDAKTFNSVTGRASIDHHFSDDLMGYIAYNRGFKSGIYNFTGIVAGTTTLPPPVNPETLDAYTIGFKSEFLNRRLRVNAEAFLYEYSDLQVQIVLPTGTTLTNASKATIKGIDFEISAIPVRNLTVSASIGLLDGRYKDFRNAPSFFLLPPNAPTPYDPARCSGPNPYPAVQGSAPATQRSCDVSGNRLIQAPPFSSNLQLLYTLPTDLGDFDLSASWAHGGNYEYDVFSSPLTRQPRTDIVNGSVQWTAPGGTFNVRLWANNLTGEKYYTIATTGSTAGAKYAAAAPRTYGVTLGANF"
attr(,"name")
[1] "fig|6666666.151743.peg.5044"
attr(,"Annot")
[1] ">fig|6666666.151743.peg.5044"
attr(,"class")
[1] "SeqFastaAA"

$`fig|6666666.151743.peg.5063`
[1] "MTVEQALARMLEGTGLEAQRLSPTAYRLVQAAPPASRRDPVRPAPPLPPEGSVDIIVTGTKREEPLARLPISVAIVQAPGEQGGALLPGTEEVAQRLDGLTLTNLGPGRNRQFIRGVADSPFNGNSQSTVAILINDARATYNAPDPDLKLVDVERVELLKGPQGSLYGTGTLGGVYHIVTRPPDAGSFAAMAGMGASLATHGEMGGSAEAMVNLPVVADRLALRAVAYAVREPGWIDTAGRGDANRGELRGGRVALRALAGDWTVDLNGMAQMRDVADSQYVYARHALERPAQPAEPHDNDFALMALRVAGKLGGADLVYAGNWVDHDVSSTLDASDLSASAFGIAPPALYLEQRRNRLFDQEVRLSQTRDDGLSWLAGLSFVDVRNHLSGVVRGVRDETVGTGIQGVVEFSAFGDMRMPLTGALSLSVGGRLFHSRIEEEQGARHARIEASVARSGFTPSLALSWQAGAHDFLYARYAGALRPGGLSRTGAASAERFESDELQTVEMGWRHSASAWLTVNLDAFLTTWSHMQSDYLLADGLIGTRNVGDARIRGVEASVDWKPAAGWAISAGLVGQSALLTRTNEGLRLEDRHLPVVPDWKARFGVNRHFDLSGWSGSMGVNLAYLGPSRLSFDPGLDRRMGDYALANADLALRRKDWSIRLSGSNLLDCKADSFAFGNPFALARGTQYTPLRPRQWAIGFTKNW"
attr(,"name")
[1] "fig|6666666.151743.peg.5063"
attr(,"Annot")
[1] ">fig|6666666.151743.peg.5063"
attr(,"class")
[1] "SeqFastaAA"

```

---
Tags: r
---


# Question
Title: Python scripting with ete3 to query NCBI's Taxonomy: "sqlite3 Warning (can only execute one statement at a time)"
I am using this script:

```
import csv
import time
import sys
from ete3 import NCBITaxa

ncbi = NCBITaxa()

def get_desired_ranks(taxid, desired_ranks):
    lineage = ncbi.get_lineage(taxid)
    names = ncbi.get_taxid_translator(lineage)
    lineage2ranks = ncbi.get_rank(names)
    ranks2lineage = dict((rank,taxid) for (taxid, rank) in lineage2ranks.items())
    return{'{}_id'.format(rank): ranks2lineage.get(rank, '<not present>') for rank in desired_ranks}

if __name__ == '__main__':
    file = open(sys.argv[1], "r")
    taxids = []
    contigs = []
    for line in file:
        line = line.split("\n")[0]
        taxids.append(line.split(",")[0])
        contigs.append(line.split(",")[1])

    desired_ranks = ['superkingdom', 'phylum']
    results = list()
    for taxid in taxids:
        results.append(list())
        results[-1].append(str(taxid))
        ranks = get_desired_ranks(taxid, desired_ranks)
        for key, rank in ranks.items():
            if rank != '<not present>':
                results[-1].append(list(ncbi.get_taxid_translator([rank]).values())[0])
            else:
                results[-1].append(rank)

    i = 0
    for result in results:
        print(contigs[i] + ','),
        print(','.join(result))
        i += 1

    file.close()

```

The script takes taxids from a file and fetches their respective lineages from a local copy of NCBI's Taxonomy database. Strangely, this script works fine when I run it on small sets of taxids (~70, ~100), but most of my datasets are upwards of 280k taxids and these break the script.

I get this complete error:

```
Traceback (most recent call last):
  File "/data1/lstout/blast/scripts/getLineageByETE3.py", line 31, in <module>
    ranks = get_desired_ranks(taxid, desired_ranks)
  File "/data1/lstout/blast/scripts/getLineageByETE3.py", line 11, in get_desired_ranks
    lineage = ncbi.get_lineage(taxid)
  File "/data1/lstout/.local/lib/python2.7/site-packages/ete3/ncbi_taxonomy/ncbiquery.py", line 227, in get_lineage
    result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %taxid)
sqlite3.Warning: You can only execute one statement at a time.

```

The first two files from the traceback are simply the script I referenced above, the third file is one of ete3's. And as I stated, the script works fine with small datasets.

**What I have tried:**

* Importing the time module and sleeping for a few milliseconds/hundredths of a second before/after my offending lines of code on lines 11 and 31. No effect.
* Went to line 227 in ete3's code...

```
result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %merged_conversion[taxid])

```

and changed the "execute" function to "executescript" in order to be able to handle multiple queries at once (as that seems to be the problem). This produced a new error and led to a rabbit hole of me changing minor things in their script trying to fudge this to work. No result. This is the complete offending function:

```
    def get_lineage(self, taxid):
    """Given a valid taxid number, return its corresponding lineage track as a
    hierarchically sorted list of parent taxids.
    """
    if not taxid:
        return None
    result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %taxid)
    raw_track = result.fetchone()
    if not raw_track:
        #perhaps is an obsolete taxid
        _, merged_conversion = self._translate_merged([taxid])
        if taxid in merged_conversion:
            result = self.db.execute('SELECT track FROM species WHERE taxid=%s' %merged_conversion[taxid])
            raw_track = result.fetchone()
        # if not raise error
        if not raw_track:
            #raw_track = ["1"]
            raise ValueError("%s taxid not found" %taxid)
        else:
            warnings.warn("taxid %s was translated into %s" %(taxid, merged_conversion[taxid]))

    track = list(map(int, raw_track[0].split(",")))
    return list(reversed(track))

```

What bothers me so much is that this works on small amounts of data! I'm running these scripts from my school's high performance computer and have tried running on their head node and in an interactive moab scheduler. Nothing has helped.

# Answer
I overlooked the fact that (for my data) about 1:50,000 results returned from my blasting yielded a ";" within the hit denoting two taxids and did not account for these in my script.

How a given researcher will respond to their own results having ";"-hits will probably vary (you may want to do extra analysis, etc.). In my case, I am only interested in the first taxid before the ";". So I added a simply python .split() to my script before appending the taxid to the list that I query:

```
if __name__ == '__main__':
    file = open(sys.argv[1], "r")
    taxids = []
    contigs = []
    for line in file:
        line = line.split("\n")[0]
        tID = line.split(",")[0]
        if ';' in tID:
            tID = tID.split(";")[0]
        taxids.append(tID)
        contigs.append(line.split(",")[1])

    desired_ranks = ['superkingdom', 'phylum']
    results = list()
    for taxid in taxids:
        results.append(list())
        results[-1].append(str(taxid))
        ranks = get_desired_ranks(taxid, desired_ranks)
        for key, rank in ranks.items():
            if rank != '<not present>':
                results[-1].append(list(ncbi.get_taxid_translator([rank]).values())[0])
            else:
                results[-1].append(rank)

    i = 0
    for result in results:
        print(contigs[i] + ','),
        print(','.join(result))
        i += 1

    file.close()

```
> 1 votes

---
Tags: python, phylogenetics, database, taxonomy, ete3
---


# Question
Title: Question about research in bioinformatics for Computer Science student
I am a computer science student, who was interested in taking research in the field of bioinformatics.

The first research idea I had was to make a model that could predict the possibility of type 2 diabetes by using machine learning such as convolutional neural networks or recurrent neural networks by training the model with fasta data of protein that causes diabetes.

The research will focus more on the architecture of the machine learning itself.

Does this research possible to work on ? If yes, how and where can I get the data? All my lecturers have no background in bioinformatic, so I can't ask any of them.

Edit: I m sorry, not protein but gene like CAPN10, CDKAL1, CDKN2A, HHEX, HNF4A, IGF2BP2, KCNJ11, PPARG, SLC30A8, and TCF7L2 that shown have risk to be associated in Type 2 Diabetes (http://diabetes.diabetesjournals.org/content/56/12/3033).

# Answer
> 1 votes
I think , this kind of research never stop until find way for cure . I have some research in heart disease,and write this paperECG, you can have some research in improve learning method to find fast and more accurate diabetes , there are lot of database your can search in google and more famous is UCI .

---
Tags: fasta, machine-learning
---


# Question
Title: Does the ".full.aln" file produced by snippy-core contain all bases of my input sequences aligned to the reference genome?
I have a number of sequences and a reference genome.

I used `snippy` to align each individual sequence with the reference genome. I then used `snippy-core` on the `*.aln` files from `snippy` which produced a new `*.aln` file, and a `*.full.aln`. It is my understanding that the new `*.aln` file contains a core SNP alignment with only bases shared by all of my sequences, and that the `*.full.aln` file contains the core genome alignment with all bases of each sequence aligned to the reference genome - regardless of variation between sequences. Is my understanding correct? If not where am I going wrong?

Further, the program `snp-sites` uses both the new `*.aln` and `*.full.aln` file to produce a file which can be run through `FastTree` to produce a phylogenetic tree. Does this just process the files to give the phylogenetic distance between each of my sequences? I don't really understand what this step achieves at all...

# Answer
`snippy-core` is a tool that processes the `snippy_outdir_{1..n}` folders and produces the following files:

1. `core.aln` – a fasta format file that contains only the polymorphic sites (of {'A', 'C', 'T' and 'G'}) in your input set, whose length tends to decrease asymptotically as more isolates are added
2. `core.full.aln` – a fasta format file that contains all the core sites and base calls at any site where the reference has a site. The reference is not padded but the sequences mapped to it may be padded with gaps and 'N's to match the length of the reference. The alignment length does not change as more isolates are added.

`snippy{-core}` was written by Torsten Seemann et al. `snp-sites` was written by Andrew Page et al. `snp-sites` will reduce any alignment to the core polymorphic sites (columns) and will still output a fasta. Therefore the output of a `core.aln` processed through `snp-sites` should be the identical the input. All fasta alignments can then be processed downstream using your tree inference package of choice (e.g., `FastTree`, `IQ-Tree`, `MrBayes`, `BEAST` ...). If you want a distance matrix output from a fasta alignment, use something like `snp-dists`, also written by Torsten Seemann, or the `APE` package in `R`. Be sure to consider how gaps and Ns are treated during this step.
> 4 votes

---
Tags: sequence-alignment, snp, phylogenetics
---


# Question
Title: What does an FDR value of 1 in RNA-seq mean?
I am looking at the supplemental data from the paper "An allelic series of miR-17 ∼ 92-mutant mice uncovers functional specialization and cooperation among members of a microRNA polycistron" which lists the genes that are differentially expressed between a particular knock-out mouse and the wild-type control. In addition to showing the logFC changes for each of the genes, the table includes their FDR value. In many cases, the FDR value is 1. What does this mean?

# Answer
FDR stands for False Discovery Rate. It is a statistic tool used in multiple hypothesis testing.

As you may know, when you use a p-value cutoff (usually 0.05) for your experiments, it means strictly speaking that "is there was actually no signal, there would be a probability of 0.05 to observe this kind of extreme values". This can be understood as "there's a 5% chance that what you call a hit is actually not signal, just some extreme value that happened by chance, and really there's nothing to see here".

This definition has deep consequences if you run many experiments: imagine that you run 100 experiments (say, you test the effect of 100 different molecules), if none of them has any effect, on average you could still get 5% of false positives. So if you run these 100 experiments, and you get, say, 10 hits with a p-value below your thresold of 0.05; how can you tell which ones are real and which ones are just chance in the absence of biological signal?

This is where you can use FDR.

If you control for a false discovery rate of, say, 0.2 (20%), it means that after a few computations (the Benjamini-Hochberg procedure described on wikipedia is simple and commonly used) you will lower your p-value cutoff so as to be sure that, *among the N hits* (probably less than the 10 you initially selected above) *that you got above, a maximum of 20% are actually false positives* (i.e. no biological effect, but surprisingly high/low signal values by pure chance).

In the article you mention they have a different FDR value for each sample, which is unusual. The values are not round numbers, so it seems likely that they are numbers obtained after applying the FDR correction (i.e. they are "transformed p-values") on which they based their decision to include the samples or not. The FDR cutoff that they chose should be written somewhere in the main text.
> 9 votes

# Answer
In practice, you can interpret it just like a p-value.

An FDR value is a p-value adjusted for multiple tests (by the Benjamini-Hochberg procedure). It stands for the “false discovery rate” it corrects for multiple testing by giving the proportion of tests above threshold alpha that will be false positives (i.e., detected when the null hypothesis is true). Note this is less stringent than the Holm-Bonferroni FWER adjusted p-value (which considers the chance of any false positives among significant results).

So for a FDR p-value of 0.05, up to 5% of these tests will be false positives (note this is not the case for uncorrected raw p-values).

For an FDR p-value of 1, up to 100% of these tests will be false positives. This makes sense since if you take test with a p-value if 1, then it should include all of the negative results, where the null hypothesis cannot be rejected. Only tests with p-values lower than this, have any chance of being rejected under the null hypothesis as true positives. If you are getting FDR adjusted p-values, you have done so many tests that you do not have the power to detect true positives and cannot rule out these being false positives due to noise of multiple testing. So an FDR p-value of 1 is definitely not significant under any circumstances.
> 2 votes

---
Tags: differential-expression
---


# Question
Title: bcftools filtering all files in a directory
Probably a silly oversight on my part, but I'm trying to filter all the vcfs in a directory with bcftools using a simple loop. My basic command is working fine:

```
bcftools filter -i 'QUAL > 1000' -o filter/file1out.vcf file1in.vcf

```

but when I try to loop it, it echoes through each filename, but generates no output files or error messages.

```
for f in *.vcf; do echo "filtering $f"; bcftools filter -i 'QUAL > 1000' -o filter/$f_out.vcf $f; done

```

Am I not handling the output direction properly? Any help much appreciated!

# Answer
> 5 votes
I assume your shell is looking for the variable `$f_out`, but can't find it because you did not define it. It probably generated a hidden file with the name `.vcf` in your output directory.

Instead, use the following: `-o ${f}_out.vcf`. This will "protect/limit" the bash variable to just `$f` rather than `$f_out`.

When encountering issues like this it's worth putting `echo` in your loop right before your command, in this case bcftools (adding quotes whenever necessary). That would show you the command as how it is going to be executed and would show you, in this case, that the output name is incorrect.

# Answer
> 2 votes
Using a `for` loop for this is by far the safer option. However, you need to quote your variables and also make sure you're passing the right variable name. While you're at it, you can also remove the extension so that your filtered file name can be `foo_out.vcf` instead of `foo.vcf_out.vcf` (which is what your loop would create. Putting all that together:

```
for f in *.vcf; do
    echo "filtering $f";
bcftools filter -i 'QUAL > 1000' -o filter/"${f//.vcf}"_out.vcf "$f"; done
fi

```

# Answer
> 1 votes
For one-liners, `xargs` is often more convenient than the `for` loop:

```
ls *.vcf | sed s,.vcf$,, \
  | xargs -i echo 'bcftools filter -i "QUAL > 1000" -o "filter/{}_out.vcf" "{}.vcf"' \
  | parallel

```

Note that your command line (with @WouterDeCoster's fix) writes output to `foo.vcf_out.vcf`. My command line writes to `foo_out.vcf`.

---
Tags: vcf, bcftools
---


# Question
Title: ANY (technical) reason behind submitting sequences to GenBank versus ENA Sequence
The DNA sequence sections of the three INSDC databases (i.e., DDBJ, ENA Sequence and GenBank) are synchronized periodically and strive to keep their stored data as ubiquitously accessible as possible. Except for idiosyncrasies in their data submission routes, there should be little, if any, reason for preferentially submitting sequence data to one database over another. Yet many researchers display very noticeable preferences!

When asked point blank, many colleagues emphasize that these preferences are not linked to the database-specific data submission or retrieval interfaces available. Yet they cannot provide a technical reason for their choice either.

Maybe researchers base their preference on geographic vicinity ("buy local"), but maybe there are genuine technical reasons for their preference (even if they are unaware of these reasons themselves).

Thus, my question is:

Can you think of any technical reason - however minute or seemingly insignificant - in preferentially submitting DNA sequences to NCBI's *GenBank* over *ENA Sequence* (or vice versa)? For example, there may be differences in data storage or accessibility that are relevant to you.

Put differently: Short of flipping a coin, why would you (as a bioinformatics-prone end user) select one database over the other for your data submission?

--

**Edit 1**: Many SE posters attempt to reassure me that I can submit my DNA sequences to either database and that "the data will be just fine there". Having submitted data to both *GenBank* and *ENA Sequence* for 10+ years myself, this is not what this question is about. Instead, this question is about carving out genuine technical differences from a user perspective.

**Edit 2**: Since no answer has been posted yet despite a bounty, I figure this is a hard question as specified. Thus, I will broaden the question and now also accept non-technical reasons for preferentially submitting DNA sequence data to either *GenBank* or *ENA Sequence*. (That being said, I still hope for someone listing a technical reason.)

# Answer
> 1 votes
*Opinion*: I prefer to use EBI/ENA because the turnaround time for "customer service" on editing metadata about a submission is faster than on NCBI. Inquiries to change metadata on NCBI can go unanswered for up to six months in my experience.

*Answer to your question*: No, there are no technical reasons that make it preferable to submit data through one website or the other. **All entries end up being hosted by both databases anyway.** The preference is up to the user to select whichever upload interface they prefer.

---
Tags: database, genomics, phylogenetics, genbank, embl
---


title: How to find genomic overlaps faster?

text:
# Question
Title: How to find genomic overlaps faster?
I was trying to find a way to find overlaps between two genomic ranges. I found a post at Biostars but it couldn't consider the chromosome information. For example:

```
library(IRanges)
df1 = data.frame(chr = c("chr1", "chr12"), start = c(10000, 10000), end = c(20000, 20000))
df2 = data.frame(chr = rep("chr1", 4), posn = c(100, 12000, 15000, 250000), x = rep(1, 4), y = rep(2, 4), z = rep(3, 4))

df1.ir = IRanges(start = df1$start, end = df1$end, names = df1$chr)
df2.ir = IRanges(start = df2$posn, end = df2$posn, names = df2$chr)

df1.it = NCList(df1.ir)

overlap_it = findOverlaps(df2.ir, df1.it)
print(overlap_it)

```

It returns the overlap from both chromosomes. I am wondering if there is a way to keep the chromosome information intact.

# Answer
> 3 votes
Just wanted to follow up to a question asked at Biostars. My solution and suggestion when dealing with genomic ranges.

```
library(GenomicRanges)
df1 = data.frame(chr = c("chr1",  "chr12"), start = c(10000,  10000), end = c(20000, 20000))
df2 = data.frame(chr = rep("chr1", 4), posn = c(100, 12000, 15000, 250000), x = rep(1, 4), y = rep(2, 4), z = rep(3, 4))

df1.ir = IRanges(start = df1$start, end = df1$end, names = df1$chr)
df2.ir = IRanges(start = df2$posn, end = df2$posn, names = df2$chr)

df1.it = NCList(df1.ir)

overlap_it = findOverlaps(df2.ir, df1.it)
print(overlap_it)

# define a genomic interval tree (git) for faster search

df1.gr = GRanges (IRanges(start = df1$start, end = df1$end), seqnames=df1$chr)
df2.gr = GRanges(IRanges(start=df2$posn, end = df2$posn), seqnames = df2$chr)

df1.git  = GNCList(df1.gr)
t1 = Sys.time ()
overlap_git = findOverlaps(df2.gr, df1.git)
t2 = Sys.time()
sub1 = difftime(t2, t1, tz, units = c("auto"))
print (sub1)
print (overlap_git)

t1 = Sys.time ()
overlap_intersect = intersect(df1.gr, df2.gr)
t2 = Sys.time ()
print(overlap_intersect)
sub2 = difftime(t2, t1, tz, units = c("auto"))
print(sub2)

```
* One should be informed about differences between IRanges and GRanges (for example chromosome context); for example, `findOverlaps` with `IRanges` returns intersections from both `chr1` and `chr12`, but with `GRanges` it doesn't
* Interval Tree based search is way faster (almost 10x) than simple `intersect` function

---
Tags: bioconductor, genomicranges
---


# Question
Title: De novo motif discovery in protein sequences
I am trying to build a features matrix to be used for Random Forest based classification. I'd like to add, as features, short motifs which are common to all the protein sequences belonging to a specific gene family.

I tried to use MEME but I'd like to know if there are also other good tools (online or command line). Keep in mind I have more than 10.000 sequences and that I prefer short motifs (3/4 aa long).

# Answer
I think MEME is a good tool for your purpose, but there are others as well. I can think of InterProScan for example (although I am not sure if it is really *de novo*). Here a summary of available tools.
> 3 votes

---
Tags: proteins, sequence-analysis, motifs, machine-learning
---


title: Bash script error at paste command

text:
# Question
Title: Bash script error at paste command
I wrote script for pasting rsids on CADD output. Here is script.

```
#!/bin/bash
cd tmp
cut -f 1,2 CADD.tsv > fileA
paste fileA <(cut -f 2,125  CADD.tsv) > myNewFile
bedtools intersect -a  myNewFile -b New.vcf -wb |cut -f 1-4,7 > CADD.rsids.tsv

```

I have tested commands one by one and they work fine. But when I run script of these commands, it gives me this error.

```
$ sh cadd.rsids.sh
4: cadd.rsids.sh: Syntax error: "(" unexpected

```

Can anyone please tell me how to resolve this error? Here are few lines of CADD.tsv file.

```
"##CADD GRCh38-v1.4 (c) University of Washington, Hudson-Alpha Institute for Biotechnology and Berlin Institute of Health 2013-2018. All rights reserved."
#Chrom  Pos Ref Alt Type    Length  AnnoType    Consequence ConsScore   ConsDetail  GC  CpG motifECount motifEName  motifEHIPos motifEScoreChng oAA nAA GeneID  FeatureID   GeneName    CCDS    Intron  Exon    cDNApos relcDNApos  CDSpos  relCDSpos   protPos relProtPos  Domain  Dst2Splice  Dst2SplType minDistTSS  minDistTSE  SIFTcat SIFTval PolyPhenCat PolyPhenVal priPhCons   mamPhCons   verPhCons   priPhyloP   mamPhyloP   verPhyloP   bStatistic  targetScan  mirSVR-Score    mirSVR-E    mirSVR-Aln  cHmm_E1 cHmm_E2 cHmm_E3 cHmm_E4 cHmm_E5 cHmm_E6 cHmm_E7 cHmm_E8 cHmm_E9 cHmm_E10    cHmm_E11    cHmm_E12    cHmm_E13    cHmm_E14    cHmm_E15    cHmm_E16    cHmm_E17    cHmm_E18    cHmm_E19    cHmm_E20    cHmm_E21    cHmm_E22    cHmm_E23    cHmm_E24    cHmm_E25    GerpRS  GerpRSpval  GerpN   GerpS   tOverlapMotifs  motifDist   EncodeH3K4me1-sum   EncodeH3K4me1-max   EncodeH3K4me2-sum   EncodeH3K4me2-max   EncodeH3K4me3-sum   EncodeH3K4me3-max   EncodeH3K9ac-sum    EncodeH3K9ac-max    EncodeH3K9me3-sum   EncodeH3K9me3-max   EncodeH3K27ac-sum   EncodeH3K27ac-max   EncodeH3K27me3-sum  EncodeH3K27me3-max  EncodeH3K36me3-sum  EncodeH3K36me3-max  EncodeH3K79me2-sum  EncodeH3K79me2-max  EncodeH4K20me1-sum  EncodeH4K20me1-max  EncodeH2AFZ-sum EncodeH2AFZ-max EncodeDNase-sum EncodeDNase-max EncodetotalRNA-sum  EncodetotalRNA-max  Grantham    Dist2Mutation   Freq100bp   Rare100bp   Sngl100bp   Freq1000bp  Rare1000bp  Sngl1000bp  Freq10000bp Rare10000bp Sngl10000bp EnsembleRegulatoryFeature   dbscSNV-ada_score   dbscSNV-rf_score    RemapOverlapTF  RemapOverlapCL  RawScore    PHRED
1   3362704 T   A   SNV 0   Transcript  INTRONIC    2   intron  0.536423841 0.04    NA  NA  NA  NA  NA  NA  ENSG00000142611 ENST00000270722 PRDM16  CCDS41236.2 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  118572  24215   NA  NA  NA  NA  0   0   0   -0.405  -1.689  -1.593  964 NA  NA  NA  NA  0   1   1   5   0   0   1   0   0   0   0   1   0   1   12  0   0   0   0   0   8   11  5   2   0   NA  NA  2.46    -4.93   1   0.02    4.53362 1.5498  3.49082 1.00224 9.30658 3.84201 5.88121 2.95189 9.31556 5.27366 4.01294 1.08327 25.99   5.10373 4.47027 1.51567 1.55624 0.7046  4.95028 1.49632 8.57929 1.45764 1.19689 0.644106    0.01783 0.01476 NA  16  1   0   15  4   10  184 17  101 1741    NA  NA  NA  5   19  -0.173019   0.424
1   7785635 T   C   SNV 0   Intergenic  DOWNSTREAM  1   downstream  0.483443709 0.053333333 NA  NA  NA  NA  NA  NA  ENSG00000049245 ENST00000054666 VAMP3   CCDS88.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  112 2704    NA  NA  NA  NA  0.009   0   0   0.418   -0.242  -0.522  726 NA  NA  NA  NA  0   0   0   3   1   0   14  1   0   0   1   0   0   0   22  0   0   0   0   0   1   1   4   0   0   NA  NA  5.52    -11 NA  NA  20.4013 3.22673 13.621  2.38936 11.7271 2.25561 9.94895 2.15059 16.5529 2.79938 11.0724 2.29032 12.2962 2.07181 32.6588 8.63119 15.114  2.76883 7.39549 2.14311 14.1811 3.08082 0.603948    0.246907    0.55734 0.16905 NA  9   1   1   15  3   5   174 35  68  1667    NA  NA  NA  2   4   0.049024    1.923
1   7785635 T   C   SNV 0   Transcript  INTRONIC    2   intron  0.483443709 0.053333333 NA  NA  NA  NA  NA  NA  ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  112 2704    NA  NA  NA  NA  0.009   0   0   0.418   -0.242  -0.522  726 NA  NA  NA  NA  0   0   0   3   1   0   14  1   0   0   1   0   0   0   22  0   0   0   0   0   1   1   4   0   0   NA  NA  5.52    -11 NA  NA  20.4013 3.22673 13.621  2.38936 11.7271 2.25561 9.94895 2.15059 16.5529 2.79938 11.0724 2.29032 12.2962 2.07181 32.6588 8.63119 15.114  2.76883 7.39549 2.14311 14.1811 3.08082 0.603948    0.246907    0.55734 0.16905 NA  9   1   1   15  3   5   174 35  68  1667    NA  NA  NA  2   4   0.049024    1.923
1   7803233 C   T   SNV 0   Transcript  INTRONIC    2   intron  0.331125828 0.013333333 NA  NA  NA  NA  NA  NA  ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  182 726 NA  NA  NA  NA  0.014   0   0   -0.553  -1.498  -0.938  673 NA  NA  NA  NA  0   0   0   6   0   0   23  5   0   0   0   0   0   0   14  0   0   0   0   0   0   0   0   0   0   NA  NA  7.49    -15 NA  NA  4.36329 1.02663 4.55863 1.4118  4.72226 1.51608 6.80823 2.09997 7.44497 2.86816 5.71457 2.35876 3.06623 0.62198 18.3251 3.23263 2.1423  0.84413 8.63178 2.70111 6.29127 1.59229 0.198671    0.048037    0.0465  0.01705 NA  6   1   1   15  2   2   134 26  39  1314    NA  NA  NA  NA  NA  -0.395046   0.067
1   7808665 T   C   SNV 0   Intergenic  DOWNSTREAM  1   downstream  0.344370861 0   NA  NA  NA  NA  NA  NA  ENSG00000236266 ENST00000451646 Z98884.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  5614    1683    NA  NA  NA  NA  0.082   0   0   0.463   -0.365  -0.373  659 NA  NA  NA  NA  0   0   0   5   3   1   24  0   0   0   0   0   0   0   15  0   0   0   0   0   0   0   0   0   0   NA  NA  5.23    -10.5   NA  NA  7.16527 1.57568 5.63104 1.29085 5.96112 1.41329 8.31159 2.11434 4.41828 1.14843 5.03919 1.48679 3.91787 1.13788 4.97324 1.34708 5.22225 1.10759 7.38244 1.40721 11.4505 2.73652 0.977252    0.19269 0.09436 0.03409 NA  43  1   0   10  3   2   110 30  49  1288    NA  NA  NA  7   26  0.11875 2.796
1   7808665 T   C   SNV 0   Transcript  INTRONIC    2   intron  0.344370861 0   NA  NA  NA  NA  NA  NA  ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  5614    1683    NA  NA  NA  NA  0.082   0   0   0.463   -0.365  -0.373  659 NA  NA  NA  NA  0   0   0   5   3   1   24  0   0   0   0   0   0   0   15  0   0   0   0   0   0   0   0   0   0   NA  NA  5.23    -10.5   NA  NA  7.16527 1.57568 5.63104 1.29085 5.96112 1.41329 8.31159 2.11434 4.41828 1.14843 5.03919 1.48679 3.91787 1.13788 4.97324 1.34708 5.22225 1.10759 7.38244 1.40721 11.4505 2.73652 0.977252    0.19269 0.09436 0.03409 NA  43  1   0   10  3   2   110 30  49  1288    NA  NA  NA  7   26  0.11875 2.796
1   7827519 C   G   SNV 0   CodingTranscript    NON_SYNONYMOUS  7   missense    0.529801325 0.066666667 NA  NA  NA  NA  P   A   ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 NA  18/22   2854    0.45172523  2590    0.712909441 864 0.714049587 lcompl  NA  NA  177 17277   tolerated   0.05    benign  0.127   0.041   0   0   -0.553  -0.852  0.026   525 NA  NA  NA  NA  1   0   0   3   0   2   8   9   1   1   1   0   0   1   2   6   1   2   9   0   0   0   0   1   0   NA  NA  12.1    -24.2   NA  NA  8.89687 3.01407 58.9843 24.0334 58.3619 18.9077 55.4967 17.2149 10.5523 3.06178 87.6613 48.8498 8.66388 2.84705 31.4762 9.05454 11.8069 3.99886 8.7429  2.26128 72.4647 18.4371 1.08886 0.638979    2.10021 1.18017 27  2   2   1   19  4   8   204 31  59  1420    Promoter    NA  NA  86  94  0.346155    6.007
1   7827519 C   G   SNV 0   Intergenic  UPSTREAM    1   upstream    0.529801325 0.066666667 NA  NA  NA  NA  NA  NA  ENSG00000236266 ENST00000451646 Z98884.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  177 17277   NA  NA  NA  NA  0.041   0   0   -0.553  -0.852  0.026   525 NA  NA  NA  NA  1   0   0   3   0   2   8   9   1   1   1   0   0   1   2   6   1   2   9   0   0   0   0   1   0   NA  NA  12.1    -24.2   NA  NA  8.89687 3.01407 58.9843 24.0334 58.3619 18.9077 55.4967 17.2149 10.5523 3.06178 87.6613 48.8498 8.66388 2.84705 31.4762 9.05454 11.8069 3.99886 8.7429  2.26128 72.4647 18.4371 1.08886 0.638979    2.10021 1.18017 NA  2   2   1   19  4   8   204 31  59  1420    Promoter    NA  NA  86  94  0.346155    6.007
1   7827519 C   G   SNV 0   RegulatoryFeature   REGULATORY  4   regulatory  0.529801325 0.066666667 NA  NA  NA  NA  NA  NA  NA  ENSR00000000832 NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  177 17277   NA  NA  NA  NA  0.041   0   0   -0.553  -0.852  0.026   525 NA  NA  NA  NA  1   0   0   3   0   2   8   9   1   1   1   0   0   1   2   6   1   2   9   0   0   0   0   1   0   NA  NA  12.1    -24.2   NA  NA  8.89687 3.01407 58.9843 24.0334 58.3619 18.9077 55.4967 17.2149 10.5523 3.06178 87.6613 48.8498 8.66388 2.84705 31.4762 9.05454 11.8069 3.99886 8.7429  2.26128 72.4647 18.4371 1.08886 0.638979    2.10021 1.18017 NA  2   2   1   19  4   8   204 31  59  1420    Promoter    NA  NA  86  94  0.346155    6.007
1   7828155 A   G   SNV 0   Intergenic  UPSTREAM    1   upstream    0.364238411 0.066666667 NA  NA  NA  NA  NA  NA  ENSG00000236266 ENST00000451646 Z98884.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  813 17022   NA  NA  NA  NA  0.009   0   0   -0.389  -0.735  -0.731  522 NA  NA  NA  NA  0   0   0   4   0   0   20  4   1   0   0   0   0   0   14  1   0   0   3   0   0   0   0   1   0   NA  NA  4.56    4.56    1   0.12    10.0753 2.35506 38.5512 20.107  21.923  10.3608 21.217  11.8934 5.61317 1.38406 9.74608 6.16783 3.76905 1.10216 24.6436 5.18488 14.8421 7.93141 4.56698 1.13989 18.7058 4.31099 0.500167    0.211141    0.16854 0.06818 NA  20  2   0   14  4   4   159 34  57  1421    NA  NA  NA  11  14  0.225659    4.402

```

Here are lines of New.vcf

```
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO
1   3362704 rs11807862  T   A   923.01  PASS    "BaseCounts=31,0,0,30;BaseQRankSum=-0.108;DB;Dels=0;FS=3.561;GC=61.35;HaplotypeScore=1.7256;MQ=60;MQ0=0;MQRankSum=-0.368;QD=15.13;ReadPosRankSum=0.498;DP=282;AF=0.5;MLEAC=1;MLEAF=0.5;AN=14;AC=7"
1   7785635 rs228729    T   C   2294.01 PASS    "BaseCounts=0,44,0,0;DB;Dels=0;FS=0;GC=42.39;HaplotypeScore=0;MQ=60;MQ0=0;QD=28.13;BaseQRankSum=2.114;MQRankSum=-1.268;ReadPosRankSum=-0.548;DP=1430;AF=0.5;MLEAC=1;MLEAF=0.5;AN=70;AC=51"
1   7803233 rs228642    C   T   2082.01 PASS    "BaseCounts=0,19,0,19;BaseQRankSum=-0.963;DB;Dels=0;FS=2.884;GC=44.89;HaplotypeScore=0.9999;MQ=60;MQ0=0;MQRankSum=0.058;QD=14.92;ReadPosRankSum=0.409;DP=1468;AF=0.5;MLEAC=1;MLEAF=0.5;AN=62;AC=44"
1   7808665 rs228666    T   C   1925.01 PASS    "BaseCounts=0,36,0,26;BaseQRankSum=0.064;DB;Dels=0;FS=0.979;GC=30.42;HaplotypeScore=1.7333;MQ=60;MQ0=0;MQRankSum=-0.82;QD=18.03;ReadPosRankSum=-1.32;DP=939;AF=0.5;MLEAC=1;MLEAF=0.5;AN=42;AC=25"
1   7827519 rs228697    C   G   786.01  PASS    "BaseCounts=0,14,24,0;BaseQRankSum=0.016;DB;Dels=0;FS=4.925;GC=58.35;HaplotypeScore=0;MQ=60;MQ0=0;MQRankSum=-0.698;QD=18.63;ReadPosRankSum=0.99;DP=355;AF=0.5;MLEAC=1;MLEAF=0.5;AN=18;AC=9"
1   7828155 rs2859388   A   G   2385.01 PASS    "BaseCounts=0,0,56,0;DB;Dels=0;FS=0;GC=35.41;HaplotypeScore=0;MQ=60;MQ0=0;QD=24.91;BaseQRankSum=-0.231;MQRankSum=0.948;ReadPosRankSum=-0.128;DP=1136;AF=0.5;MLEAC=1;MLEAF=0.5;AN=46;AC=31"
1   7830057 rs2640909   T   C   1421.01 PASS    "BaseCounts=0,19,0,19;BaseQRankSum=1.051;DB;Dels=0;FS=5.041;GC=56.11;HaplotypeScore=0.734;MQ=60;MQ0=0;MQRankSum=0.496;QD=16.61;ReadPosRankSum=-0.526;DP=638;AF=0.5;MLEAC=1;MLEAF=0.5;AN=32;AC=19"
1   11107089    rs12139042  G   A   1668.01 PASS    "BaseCounts=43,0,0,0;DB;Dels=0;FS=0;GC=45.64;HaplotypeScore=0.734;MQ=60;MQ0=0;QD=29.7;BaseQRankSum=-1.652;MQRankSum=0.027;ReadPosRankSum=-1.192;DP=354;AF=0.5;MLEAC=1;MLEAF=0.5;AN=18;AC=10"
1   11787392    rs3737967   G   A   719.01  PASS    "BaseCounts=23,0,23,0;BaseQRankSum=-0.475;DB;Dels=0;FS=4.566;GC=61.85;HaplotypeScore=0.734;MQ=60;MQ0=0;MQRankSum=-0.84;QD=15.63;ReadPosRankSum=0.329;DP=378;AF=0.5;MLEAC=1;MLEAF=0.5;AN=20;AC=10"
1   11790870    rs2274976   C   T   777.01  PASS    "BaseCounts=0,18,0,24;BaseQRankSum=-0.042;DB;Dels=0;FS=4.677;GC=56.11;HaplotypeScore=0.6651;MQ=60;MQ0=0;MQRankSum=1.204;QD=15.05;ReadPosRankSum=0.836;DP=355;AF=0.5;MLEAC=1;MLEAF=0.5;AN=20;AC=10"
1   11792243    rs1476413   C   T   1751.01 PASS    "BaseCounts=1,18,0,21;BaseQRankSum=1.028;DB;Dels=0;FS=0;GC=56.11;HaplotypeScore=1.7287;MQ=60;MQ0=0;MQRankSum=-0.408;QD=16.65;ReadPosRankSum=-0.099;DP=954;AF=0.5;MLEAC=1;MLEAF=0.5;AN=50;AC=29"
1   11794400    rs4846051   G   A   2616.01 PASS    "BaseCounts=46,0,0,0;DB;Dels=0;FS=0;GC=59.35;HaplotypeScore=6.2573;MQ=60;MQ0=0;QD=31.75;DP=1654;AF=1;MLEAC=2;MLEAF=1;AN=78;AC=78"
1   11794419    rs1801131   T   G   2028.01 PASS    "BaseCounts=0,1,21,21;BaseQRankSum=-0.117;DB;Dels=0;FS=9.995;GC=59.35;HaplotypeScore=2.4328;MQ=60;MQ0=0;MQRankSum=1.03;QD=14.84;ReadPosRankSum=0.665;DP=1176;AF=0.5;MLEAC=1;MLEAF=0.5;AN=54;AC=33"

```

# Answer
Your problem is that you have written a script in `bash` but are then running it using `sh`. Bash, the Bourne-again shell, is not the same as the Bourne shell (`sh`). What is more, on Ubuntu, `sh` is actually not even the Bourne shell but another minimal shell called `dash`. The `/bin/sh` on Ubuntu is a symlink to `dash`:

```
$ ls -l /bin/sh
lrwxrwxrwx 1 root root 4 Mar  1  2018 /bin/sh -> dash

```

Dash is a POSIX compliant shell which means it has the features defined by the POSIX specification for the `sh` shell. In your case, the problem is caused by the fact that the `<()` input redirection construct is supported by `bash` but not by `dash`. So, in order to run your script, you can do one of two things:

1. Make the script executable (`chmod a+x cadd.rsids.sh`), and then run it directly:

```
/path/to/cadd.rsids.sh

```

This is the simplest approach since you already have the shebang line (`#!/bin/bash`) which will cause the script to be run by `bash` and not `sh`.
2. Run the script with `bash` instead of `sh`:

```
bash /path/to/cadd.rsids.sh

```

Instead of

```
sh /path/to/cadd.rsids.sh

```
> 5 votes

---
Tags: vcf, shell
---


# Question
Title: How to output all sequences with bwa mem, not `*`?
I've been running `bwa mem -a` for alignment, using the `-a` flag---this will

> output all alignments for SE or unpaired PE

I've noticed in the SAM that there are several alignments with `*` in the `SEQ` and `QUAL` fields. Based on the documentation:

> 10. SEQ: segment SEQuence. This field can be a ‘\*’ when the sequence is not stored. If not a ‘\*’, the length of the sequence must equal the sum of lengths of M/I/S/=/X operations in CIGAR. An ‘=’ denotes the base is identical to the reference base. No assumptions can be made on the letter cases.
> 11. QUAL: ASCII of base QUALity plus 33 (same as the quality string in the Sanger FASTQ format). A base quality is the phred-scaled base error probability which equals −10 log10 Pr{base is wrong}. This field can be a ‘\*’ when quality is not stored. If not a ‘\*’, SEQ must not be a ‘\*’ and the length of the quality string ought to equal the length of SEQ.

it appears the sequence isn't stored.

I would strongly prefer to have the sequence in this case. Is there any to direct bwa to output these sequences?

# Answer
I don't believe that is an option. Internally we created a patch and submitted it to BWA however that was 4 years ago and it has not been accepted.
> 5 votes

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


title: Install DESeq2 through anaconda

text:
# Question
Title: Install DESeq2 through anaconda
I tried running `conda install -c bioconda bioconductor-deseq2` in a `conda` environment, but when I run `R` console and try to import `DESeq2` I am getting the error message:

> Error: package or namespace load failed for ‘GenomeInfoDb’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI\[\[i\]\]): there is no package called ‘GenomeInfoDbData’ Error: package ‘GenomeInfoDb’ could not be loaded

When I try installing `GenomeInfoDbData` and `GenomeInfoDb` using `conda`, it says that the packages were already installed. If I run `install.packages(GenomeInfoDb)` from `R` console, it is telling me that the package is not available for `R` version `3.5.1`.

I want to import `DESeq2` in `Rstudio` though and it is giving me a different error when I try importing `DESeq2`:

> Error in library("DESeq2") : there is no package called ‘DESeq2’

I do not get why `Rstudio` does not see that I alrerady installed `DESeq2` with `conda`. If I try running `biocLite("DESeq2")` from `Rstudio` I get the following error:

> Error in dyn.load(file, DLLpath = DLLpath, ...) : unable to load shared object '/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so': dlopen(/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so, 6): Symbol not found: \_libintl\_dgettext Referenced from: /Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so Expected in: /Users/nikitavlasenko/miniconda3/envs/renv/lib/R/lib/libR.dylib in /Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so ERROR: lazy loading failed for package ‘Hmisc’ * removing ‘/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/Hmisc’ ERROR: dependency ‘Hmisc’ is not available for package ‘DESeq2’ * removing ‘/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/DESeq2’

So, supposedly the issue is with `Hmisc`. I installed the package successfully with `conda`, but `Rstudio` is apparently does not know about it. Running `install.packages("Hmisc")` gives the first part of the error with missing `cluster.so` library.

Any suggestions would be greatly appreciated.

# Answer
Ultimately my colleague helped me to solve the issue by following the steps:

1. Created environment: `conda create --name myenv`
2. Activated it: `source activate myenv`
3. Edited `.condarc` to set the priority of the channels to be the following (top to bottom): `conda-forge`, `bioconda`, `r`, `default`
4. Checked that the channels are set in the correct order: `conda config --get channels`
5. `conda install r=3.5.1`
6. `conda install r-essentials`
7. `conda install rstudio`
8. Installation of all of the `bioconductor` packages through `bioconda` and never install anything from inside `R` or `Rstudio` (I guess that is how it was broken)

Then launching `rstudio` from within the environment. Maybe it is not a kind of quick fix and it could have been solved easier, but I was not able to.
> 5 votes

# Answer
Sounds like you might have an issue with which `R` Rstudio is running. If you didn't install Rstudio through `conda`, Rstudio may be using its own version of `R` that is bundled with Rstudio, instead of the one you've installed via `conda`.

Assuming that your `conda` environment name is `renv`, try running this in the terminal:

```
$ conda activate renv
$ conda install r::rstudio -y

$ which rstudio
/Users/nikitavlasenko/miniconda3/envs/renv/bin/rstudio
# or something like this ^

# run Rstudio as a background shell process
rstudio &

```

This should open up the Rstudio interface, like normal, but using everything defined in your `renv` environment.

If you try loading the `DEseq2` library now, that might work.
> 2 votes

# Answer
As others have mentioned, this is likely an issue with having multiple versions of R or RStudio installed.

I highly recommend that any R/RStudio version not installed inside conda be removed. This includes any installed libraries.

I've previously written a blog post about this topic, and I can say from experience that installing R packages through conda can cause lots of headaches.

You'll have an easier time if you install R/RStudio via conda and then install R packages the traditional way, inside R/RStudio.

So, for example.

Bad: `conda install -c bioconda bioconductor-deseq2`

Good:

```
conda install rstudio
# Then, inside RStudio
source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2")

```
> 2 votes

---
Tags: deseq2
---


title: Real time transcript profiles

text:
# Question
Title: Real time transcript profiles
Do any methods exist (or are in the process of development) for investigating transcript data without lysing the cells, i.e, destroying the sample?

# Answer
> 3 votes
**Answer:** There are no methods that exist for investigating transcript data inside a cell without destroying the cell. Separating RNA from cells requires essentially destroying the cells and using a chemical cocktail to split its constituent components.

**Alternative 1:** One option is to investigate extracellular transcripts of living cells is to centrifuge cells in culture and isolating RNA from the growth media supernatant \[1\]. Cell-free DNA and cell-free RNA are hot topics in contemporary biology and bioinformatics as a potential source for disease biomarkers. \[2\]

**Alternative 2:** Another option to monitor a transcript from a living cell is to use a transgenic construct with a reporter molecule like luciferase or a fluorescent protein to produce light or become fluorescent when a promotor is activated \[, \]. This does not directly measure transcript quantity but answers a similar question: how does transcription change in a living cell?

Sources:

\[1\] Théry, C., Amigorena, S., Raposo, G., & Clayton, A. (2006). Isolation and characterization of exosomes from cell culture supernatants and biological fluids. Current protocols in cell biology, 30(1), 3-22.

\[2\] Yu, X. M., Wu, Y. C., Liu, X., Huang, X. C., Hou, X. X., Wang, J. L., ... & Ling, Z. Q. (2016). Cell-Free RNA content in peripheral blood as potential biomarkers for detecting circulating tumor cells in non-small cell lung carcinoma. International journal of molecular sciences, 17(11), 1845.

\[3\] Kaskova, Z. M., Tsarkova, A. S., & Yampolsky, I. V. (2016). 1001 lights: luciferins, luciferases, their mechanisms of action and applications in chemical analysis, biology and medicine. Chemical Society Reviews, 45(21), 6048-6077.

\[4\] Rao, J., Dragulescu-Andrasi, A., & Yao, H. (2007). Fluorescence imaging in vivo: recent advances. Current opinion in biotechnology, 18(1), 17-25.

---
Tags: rna-seq, transcriptome, rna
---


title: Linking of MeSH terms and genes

text:
# Question
Title: Linking of MeSH terms and genes
Is there an alternative to mesh2gene for linking MeSH terms to genes?

The http://www.ncibi.org/gene2mesh.html site appears to be down for months now (click on "Launch Gene2MeSH via gene2mesh.ncibi.org") and I've been unable to find an alternative so far...

# Answer
> 2 votes
I think you could use the MESH.db package in Bioconductor to annotate your genes. I have never used, but you should be able to do something like:

```
library("MeSH.db")
select(MeSH.db, keys = myGenes, keytype = "ENTREZ", columns = "MESH")

```

It depends on the MeSHDbi package from which you could retrieve other data form that database if needed

---
Tags: gene
---


# Question
Title: ncbi...eutils/esearch.fcgi: How to search now in the snp database?
According to the 2017(!) textbook, it worked at that time, for instance,

`https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=autism`

At this time, this does not work and for each conceivable term it will find 0 records.

What is the right way to search in `snp`, valid in December 2018?

# Answer
> 3 votes
eUtilities for dbSNP is in the process of being retired. See this post from NCBI. You may want to check out the new API.

From the NCBI Insights blog post:

> Entrez search is available only for human. We will retire all Eutils report formats (RS docsum XML and ASN.1, FASTA, and FLAT text) later this year. Only eSummary XML report will be supported.
>
> Batch Query Service for all organisms will be retired October 2018. Users of human dbSNP data should transition to the new API.
>
> We will place non-human FTP data in the archive directory. Please see the phase-out plan announced last year.

---
Tags: snp, public-databases, phylogenetics
---


# Question
Title: Understanding the initialization of the Needleman-Wunsch algorithm
I would like to ask why when we initialize 0 as the starting point, there are 3 gray boxes that are not used as shown (highlighted in the red box) in the pic?  Are these boxes to show that they are gaps or just to align them properly?

# Answer
> 7 votes
The first column is added in order to be able to align the sequences. They might be gaps if the alignment ends up beginning in different positions than the first. Without the first empty row and column it wouldn't be possible to have gaps at the beginning of the alignment

The dynamic programming is useful because you need first to calculate the score for each position, and then once you have the score, you can calculate the optimum path. So you solve your problem by first calculating what would happen if (there was a gap, there was a mismatch...) and then you find the optimum alignment.

---
Tags: sequence-alignment
---