title: FASTQ file trimming

text:
# Question
Title: FASTQ file trimming
I have a FASTQ file:

```
@NS500455:80:HG7TNBGXB:1:11101:17723:1055 1:N:0:ATCACG
ACTTANGTGTATGTAAACTTCCGACTTCAACTGTATAGGGATCCNAGCTCCAATTCGCCCTATAGTGAGTCGTAT
+
/AAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@NS500455:80:HG7TNBGXB:1:11101:8821:1057 1:N:0:ATCACG
ACTTANGTGATGTAAACTTCCGACTTCAACTGTATAATAAATATCTAGATCGGAAGAGCACACGTCCGAACTCCA

```

I want to remove the string ACTTAAGTGTATGTAAACTTCCGACTTCAACTG from it while retaining the sequences starting from TA onwards. I wrote `grep "ACTTAAGTGTATGTAAACTTCCGACTTCAACTGTA" SRR_1.fastq | sed "s/ACTTAAGTGTATGTAAACTTCCGACTTCAACTGTA/TA/g"`, but it removes the headers of the FASTQ file as well.

The desired output would be

```
@NS500455:80:HG7TNBGXB:1:11101:17723:1055 1:N:0:ATCACG
TATAGGGATCCNAGCTCCAATTCGCCCTATAGTGAGTCGTAT
+
/AAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@NS500455:80:HG7TNBGXB:1:11101:8821:1057 1:N:0:ATCACG
TATAATAAATATCTAGATCGGAAGAGCACACGTCCGAACTCCA

```

How to do this. I dont want to use cutadapt.

Kindly help.

# Answer
> 4 votes
What application are you using that will be okay with a mismatching quality string?

Use trimmomatic or cutadapt, or bcl2fastq could probably do the job too. Those will also have the benefit of dealing properly with one-off errors, which awk and sed won't.

# Answer
> 3 votes
Skip the `grep`, use just the `sed`. Unless the header contains the long sequence as well (which I doubt it will), your command will work.

```
sed 's/^ACTTAAGTGTATGTAAACTTCCGACTTCAACTGTA/TA/' SRR_1.fastq > SRR_1.edited.fastq

```

# Answer
> 2 votes
If you need the quality vector for each read, then your code will create a shift because it does not cut the quality vector. You should rather use a trimming tool a systematically cut all 33 first characters. **fastp** will allow you to do it easily with the `-f, --trim_front1` argument.

---
Tags: fastq
---


title: Files for paired end RNA sequencing

text:
# Question
Title: Files for paired end RNA sequencing
I am looking at the videos at a DIY Transcriptomics course and the speaker mentions that to run Kallisto for read alignment with paired end sequencing, one would enter the following:

where the sample files at the end refer to the two different end reads (forward and reverse read files). However, browsing the NCBI for paired end data, I noticed that only 1 Fastq file is provided per replicate for a paired end procedure (here's an example). So I'm a bit confused about this.

# Answer
> 0 votes
Building off of the answer by @An Ignorant Wanderer:

Essentially there are indeed the paired end data files available, but it takes a little bit of extra navigation to reach it.

Using the specific example you provided of BioProject PRJNA574273: First select a specific sequencing "Run," aka the sequencing of a specific sample / experimental condition. For example, click on the link to run: SRR10189234

This brings you to the default tab: "Metadata." As @An Ignorant Wanderer pointed out, then you can navigate to the "Data Access" tab.

Then look under the "Original format" subheading. Here you'll find links to two files, corresponding to the two orientations of the paired end sequence files.

Edit: Also please see the answer by @burger, the website linked in that answer is much more user friendly in my opinion.

# Answer
> 1 votes
The problem is that the SRA website can be a little confusing. View the same project on ENA where the organization is more human-readable: https://www.ebi.ac.uk/ena/browser/view/PRJNA574273

# Answer
> 0 votes
Actually after clicking on one of the "Runs" and going to "data access" I noticed two files are provided per run

# Answer
> 0 votes
Paired ends is a configuration in sequencing platforms. Illumina is the most popular at the moment and most single-cell sequencing is done using paired-end Illumina sequencing. It varies between single-cell technologies but most 3’ scRNA-Seq protocols use paired ends.

Kallisto is a tool for handling raw data. Many databases supported processed data in BAM format or UMI/read counts. You should not need to run kallisto unless you are involved in handling raw data from sequencing performed by your laboratory. Some databases such as SRA and EGA support raw data in FASTQ format, others do not for ethical reasons or the data is too large.

I recommend downloading from SRA via their command line interface (this is documented on their website). It is possible to extract paired ends from an aligned BAM file using bedtools and samtools. I’ve done this with the Macosko et al (Cell 2015) DropSeq data for example.

Note that in some cases, the barcode is stored in databases as the headers of the FASTQ, BAM, BED, or CTSS file rather than a separate read. Single cell transcriptomics is an emerging field so there isn’t a standard way to do this, it varies between labs and data centres. For single-cell specifically, some databases are being set up, for example: https://data.humancellatlas.org/

---
Tags: rna-seq, kallisto
---


# Question
Title: Understanding the "relativeness" in RNA seq experiments
Let's say we have two samples, sample 1 and sample 2, and we're running an RNA seq experiment. Consider genes A and B for both. For simplicity, let's say that when we run the RNA seq experiment, we get the following results:

gene A sample 1 count: 2

gene B sample 1 count: 4

gene A sample 2 count: 2000

gene B sample 2 count: 4000

Based on my understanding, RNA seq gives relative results. So I guess in this guess, how can you tell whether either one of the following are true:

1- You sequenced sample 2 a thousand fold more compared to sample 1

2- Sample 2 really has a thousand fold more expression in both genes A and B

# Answer
> 2 votes
1. We cannot be sure about this. For all we know, they could have been sequenced to the same depth and there is a biological difference that drives the 1000x (ridiculously large) fold change.
2. See above

The only thing you can say for sure is that gene B is expressed twice as much as gene A in samples 1 and 2. That's what "relative" means - relative expression of genes within a sample.

# Answer
> 1 votes
You obviously have to normalize for differences in depth followed by a proper statistical test to get fold changes and significances. Sequencing depth is a technical confounder that has to be removed, and that is what normalization does. You have been posting a lot of very basic RNA-seq questions, I suggest you read the RNA-seq Hitchhikers guide: https://peerj.com/preprints/27283/ before posting new questions first.

# Answer
> 0 votes
The process of RNA-seq is a sampling process of all the RNA molecules. We cannot make reliable comparisons before we know how deep we sample. And since we know it is a sample from a population, the absolute value does not make sense.

---
Tags: rna-seq
---


# Question
Title: How to use the DI-tector software for defective genomes analysis?
I am trying to use a free .py script from this site

```
  http://www.di-tector.cyame.eu/

```

while unfortunately there's neither manual nor any documentation presented.

I have downloaded the .py script however I have absolutely no idea how it can be started. Can anyone give me an idea of how that can be done?

The corresponding article (https://pubmed.ncbi.nlm.nih.gov/30012569/) only displays parameters for executing it and there is no way to understand what steps must be performed firstly. Thank you.

# Answer
> 1 votes
See the help:

```
python DI-tector_06.py -h
usage: DI-tector_06.py [-h] [-g HOST_REF] [-s MIN_SEGMENT] [-m MIN_MAPQ]
                       [-n NB_READS] [-o OUTPUT_DIRECTORY] [-t TAG] [-d]
                       [-l INDEL_LENGTH] [-f] [-p POLARITY] [-q] [-k]
                       [-x NB_THREADS]
                       Virus_Ref Input_Data

positional arguments:
  Virus_Ref             Virus genome reference sequence in FASTA format.
  Input_Data            File containing single reads in FASTQ format.

optional arguments:
  -h, --help            show this help message and exit
  -g HOST_REF, --Host_Ref HOST_REF
                        Host genome reference sequence in FASTA format.
  -s MIN_SEGMENT, --Min_Segment MIN_SEGMENT
                        Minimum segment length. Default is 15.
  -m MIN_MAPQ, --Min_MAPQ MIN_MAPQ
                        Skip alignments with MAPQ smaller than INT. Default is
                        25.
  -n NB_READS, --Nb_Reads NB_READS
                        Show only DVGs with counts reads > or egal to INT.
                        Default is 2.))
  -o OUTPUT_DIRECTORY, --Output_Directory OUTPUT_DIRECTORY
                        Enter a directory name that all compiled output files
                        will be saved in.
  -t TAG, --Tag TAG     Enter a tag name that will be appended before each
                        output file. Default is 'DI-tector'.
  -d, --DVG_sequences   Generate multi-fasta file with DVG sequences. Default
                        is (OFF).
  -l INDEL_LENGTH, --InDel_Length INDEL_LENGTH
                        Skip alignments with size of InDels smaller or egal to
                        INT. Default size is 1.
  -f, --Fasta           Select '-f' if data is in FASTA format fasta. Default
                        is FASTQ.
  -p POLARITY, --Polarity POLARITY
                        [0] Positive strand genome / [1] Negative strand
                        genome. Default is 0.
  -q, --No_Quantification
                        Inactive percentage quantification. Quantification
                        need bedtools Default is (ON).
  -k, --Keep_files      Keep intermediaite files (i.e. alignment etc...).
                        Default is (OFF).
  -x NB_THREADS, --Nb_threads NB_THREADS
                        Number of threads. Default is 1.

```

---
Tags: python, genome
---


# Question
Title: Parsing a GenBank file with multiple gene entries
Basically a GenBank file consists of gene entries (announced by 'gene') followed by its corresponding 'CDS' entry (only one per gene) like the two shown here below. I would like to extract part of the data from the input file shown below according to the following rules and print it in the terminal.

* There are two blocks of gene data shown below. I want to extract part of both blocks.
* The extracted text for each block starts with a line that contains spaces at the beginning of the line followed by `gene`
* The extracted text for each block ends with a line that contains `/db_xref="GeneID`
* The perl and awk tags are just suggestions. An answer can use a different program(s).

Input file (source: biotech):

```
     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
                     /translation="MYYKALAHFLPTLSTMQNILSKSPLSLDFRLLFLAFIDKR"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"
                     /translation="MQTVNVDVAIVGAGGGGLRAAIAAAEANPNLKIALISKVYPMRS
                     HTVAAEGGAAAVAKEEDSYDKHFHDTVAGGDWLCEQDVVEYFVEHSPVEMTQLERWGC
                     PWSRKADGDVNVRRFGGMKIERTWFAADKTGFHLLHTLFQTSIKYPQIIRFDEHFVVD
                     ILVDDGQVRGCVAMNMMEGTFVQINANAVVIATGGGCRAYRFNTNGGIVTGDGLSMAY
                     RHGVPLRDMEFVQYHPTGLPNTGILMTEGCRGEGGILVNKDGYRYLQDYGLGPETPVG
                     KPENKYMELGPRDKVSQAFWQEWRKGNTLKTAKGVDVVHLDLRHLGEKYLHERLPFIC
                     ELAQAYEGVDPAKAPIPVRPVVHYTMGGIEVDQHAETCIKGLFAVGECASSGLHGANR
                     LGSNSLAELVVFGKVAGEMAAKRAVEATARNQAVIDAQAKDVLERVYALARQEGEESW
                     SQIRNEMGDSMEEGCGIYRTQESMEKTVAKIAELKERYKRIKVKDSSSVFNTDLLYKI
                     ELGYILDVAQSISSSAVERKESRGAHQRLDYVERDDVNYLKHTLAFYNADGTPTIKYS
                     DVKITKSQPAKRVYGAEAEAQEAAAKKE"

```

Desired output:

```
     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

```

I tried using `pcregrep --multiline .*'START-SEARCH-TERM.*(\n|.)*END-SEARCH-TERM' path/to/SOURCE-FILE`. I installed pcregrep (grep utility that uses Perl-style regexps) in Ubuntu with `sudo apt install pcregrep`. I also installed Biopython with `sudo apt install python3-biopython` and ran the Simple GenBank parsing example from *Biopython Tutorial and Cookbook*.

# Answer
> 1 votes
You can simply use grep for this purpose as shown below,

```
grep -v /translation bio.txt | grep -B100000000 /db_xref= > output_file.txt

```

Just make sure that you keep the number with B bigger than the number of lines of your file.

If you print the contents of the above file you get your desired output as given below,

```
:~$ cat output_file.txt
gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

```

# Answer
> 1 votes
You can do this very easily with `awk`:

```
$ awk '/^ {5}\w/{a=1} /\/translation/{a=0}a' file.gb
     gene            complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /db_xref="GeneID:7278619"
     CDS             complement(8972..9094)
                     /locus_tag="HAPS_0004"
                     /codon_start=1
                     /transl_table=11
                     /product="hypothetical protein"
                     /protein_id="YP_002474657.1"
                     /db_xref="GI:219870282"
                     /db_xref="GeneID:7278619"
     gene            9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /db_xref="GeneID:7278620"
     CDS             9632..11416
                     /gene="frdA"
                     /locus_tag="HAPS_0005"
                     /note="part of four member fumarate reductase enzyme
                     complex FrdABCD which catalyzes the reduction of fumarate
                     to succinate during anaerobic respiration; FrdAB are the
                     catalytic subcomplex consisting of a flavoprotein subunit
                     and an iron-sulfur subunit, respectively; FrdCD are the
                     membrane components which interact with quinone and are
                     involved in electron transfer; the catalytic subunits are
                     similar to succinate dehydrogenase SdhAB"
                     /codon_start=1
                     /transl_table=11
                     /product="fumarate reductase flavoprotein subunit"
                     /protein_id="YP_002474658.1"
                     /db_xref="GI:219870283"
                     /db_xref="GeneID:7278620"

```

The idea here is to set `a` to `1` if this line starts with 5 spaces followed by a word character. Then, we set `a` back to `0` if this line matches `/translation`. The default action for `awk` when an expression evaluates to true (not 0) is to print, therefore the final `a` will cause all lines read while `a` is not 0 to be printed, effectively removing everything after each `/translation` line.

A more easily understandable version of the same code would be:

```
awk '{
        if(/^ {5}\w/){
            a=1
        }
        if(/\/translation/){
            a=0
        }
        if(a==1){
            print
        }
    }' file.gb

```

---
Tags: python, text-processing, genbank
---


# Question
Title: Non-exome variants called from whole-exome sequencing data
I'm working on rather standard whole-exome sequencing data and treat it the same as whole-genome sequencing data (aligning it to the full GRCh38 reference assembly and calling variants with no exome-related parameters).

I recently noticed that my variant calling results contain numerous variants mapped to a non-exomic region, incidentally, some of them were significantly associated with the phenotype I study.

My questions are as follows:

1. Are the non-exome variants in my exome data invalid? or could it be that the exome kit captured more than it was designed to capture?
2. Should I treat raw exome sequencing data any different than genome sequencing data, and if so, what parts of my pipeline should be modified to accommodate exome data?

# Answer
> 1 votes
These are called off-targets due to unspecific binding of DNA to the capturing beads during library preparation. This is normal and expected since capturing is not a perfect process. It is correct to align against the whole genome but you should simply limit variant calling to exons (or basically the on-target regions) and ignore everything else.

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


title: Finding out-group in unrooted tree

text:
# Question
Title: Finding out-group in unrooted tree
Im trying to find the out-group in the alignment my data files file a and file b both the case i'm getting NA when i run this `pars.a$tip.label[pars.a$edge[which.max(pars.a$edge.length),2]]`

I tried NJ instead of upgma i do still get NA not sure what is the issue ,but at the same time when i used other alignment file with 50 sequence works fine.

```
a<-read.alignment("abhi_seq/final_abhi_tempa1.fas", format="fasta")
a.phydat<-as.phyDat(a)
dist.a.phydat<-dist.dna(as.DNAbin(a.phydat))
upgma.a<-upgma(dist.a.phydat)
parsimony(upgma.a,a.phydat)
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
pars.a$tip.label[pars.a$edge[which.max(pars.a$edge.length),2]]

```

# Answer
> 1 votes
the fact is that the present upgma tree as a midpoint rooting here, corresponding to node 46.

```
getRoot(upgma.a)

```

so for this kind of tree, I recommend you to use the midpoint() rooting function.

```
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
plot(midpoint(pars.a))

```

---
Tags: r, phylogenetics
---


# Question
Title: Combine structures in the same session - PyMOL
I'm a student currently looking at antibody responses against a viral target protein of interest.

I have my own, annotated PyMOL session of my protein and I also have .pdbs of crystallised antibody structures that bind to this protein and rougly know whereabouts it binds. I would like to know if it is possible to combine these two structures into a single session to observe or at least overlay the two structures together. Many thanks!

# Answer
> 1 votes
1. Open the PyMOL session.
2. Use the `load` command to load your PDB file(s): `load file.pdb`.
3. Use the `align` or `super` commands to overlay structures, e.g. `align prot1, prot2`.

---
Tags: structural-biology, pymol
---


# Question
Title: Sophisticated algorithms in biological systems
I am looking for examples of (known) algorithms/data structures in biological systems. My current understanding is that these don't exist - or haven't been found - beyond a relatively low level of complexity, which would seem like strong evidence to me that evolution isn't able to come up with something as complicated (and multi-step) as, e.g., the Schönhage-Strassen algorithm.

I heard of Bloom filters for recognizing smells in fruit flies (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6304992/).

# Answer
There are many examples of biological systems doing computation, according to my understanding:

1. Regulatory networks encode decision trees of significant complexity: wikipedia.
2. Alan Turing did significant work on organismal morphogenesis showing how simple chemical patterns of inducers can lead to complex outcomes, "Turing patterns".
3. All of organismal development, cell division, etc. are tightly regulated procedures of gene regulation and morphogenesis.
4. The inspiration for neural networks is brains.
> 2 votes

---
Tags: microbiology, systems-biology
---


# Question
Title: output Array[Pair[X,Y]] to file inside a WDL task
I have a table containing different metadata information about a set of samples as an Array\[Object\] and I need to extract a couple of columns from this table and output those two columns into a file that will be read and used by the task. I wanted to know whether I can do this within a task rather than within a workflow, that is, without using the scatter construct which is not allowed inside a task.

If I zip two arrays I obtain an Array\[Pair\[String,String\]\] variable.

Now if I want to output to file an Array\[Array\[String\]\] I can use write\_tsv(), if I want to output to file a Map\[String,String\] I can use write\_map().

But WDL does not let me coerce Array\[Pair\[String,String\]\] to either Array\[Array\[String\]\] nor Map\[String,String\] (though I can coerce Map\[String,String\] to Array\[Pair\[String,String\]\]), so neither option is available.

Is there a solution to this that does not involve the scatter construct?

# Answer
`write_json()` is supposed to provide the general way to serialize any data structure for processing/manipulation in the task command script (perhaps leveraging jq). Unfortunately, last I heard, Cromwell's implementation doesn't work on some useful cases right now. (Possibly out of date; I haven't tried lately)
> 0 votes

---
Tags: wdl
---


# Question
Title: Comparing multiple treatments to multiple other treatments in edgeR for simple effects in a complex experimental design
I am working with a RNA-seq data set in maize that has a relatively complex design. There are two levels of treatment A (nitrogen fertilizer level in the field, high or low), two levels of treatment B (nitrogen nutrients in *in vitro* cultures, high and low) and two levels of treatment C (two time points of sampling), all with 3 reps.

```
> library(edgeR)
> load("KC_Raw.RData")
> y <- DGEList(counts = KCraw.data[,2:25])
> keep <- rowSums(cpm(y) > 10) >= 3
> targets <- data.frame(rownames=colnames(KCraw.data)[2:25] ,
+                       Time=rep(c(rep("2DIC",12),rep("5DIC",12))) ,
+                       FieldN=rep(c(rep("FH",6), rep("FL",6)),2) ,
+                       CultureN=rep(c(rep("CL",3),rep("CH",3)),4))
> Group <- factor(paste(targets$FieldN,targets$Time,targets$CultureN,sep="."))
> targets <- cbind(targets,Group=Group)
> targets
   rownames Time FieldN CultureN      Group
1   KC1_H2L 2DIC     FH       CL FH.2DIC.CL
2   KC2_H2L 2DIC     FH       CL FH.2DIC.CL
3   KC3_H2L 2DIC     FH       CL FH.2DIC.CL
4   KC4_H2H 2DIC     FH       CH FH.2DIC.CH
5   KC5_H2H 2DIC     FH       CH FH.2DIC.CH
6   KC6_H2H 2DIC     FH       CH FH.2DIC.CH
7   KC7_L2L 2DIC     FL       CL FL.2DIC.CL
8   KC8_L2L 2DIC     FL       CL FL.2DIC.CL
9   KC9_L2L 2DIC     FL       CL FL.2DIC.CL
10 KC10_L2H 2DIC     FL       CH FL.2DIC.CH
11 KC11_L2H 2DIC     FL       CH FL.2DIC.CH
12 KC12_L2H 2DIC     FL       CH FL.2DIC.CH
13 KC13_H5L 5DIC     FH       CL FH.5DIC.CL
14 KC14_H5L 5DIC     FH       CL FH.5DIC.CL
15 KC15_H5L 5DIC     FH       CL FH.5DIC.CL
16 KC16_H5H 5DIC     FH       CH FH.5DIC.CH
17 KC17_H5H 5DIC     FH       CH FH.5DIC.CH
18 KC18_H5H 5DIC     FH       CH FH.5DIC.CH
19 KC19_L5L 5DIC     FL       CL FL.5DIC.CL
20 KC20_L5L 5DIC     FL       CL FL.5DIC.CL
21 KC21_L5L 5DIC     FL       CL FL.5DIC.CL
22 KC22_L5H 5DIC     FL       CH FL.5DIC.CH
23 KC23_L5H 5DIC     FL       CH FL.5DIC.CH
24 KC24_L5H 5DIC     FL       CH FL.5DIC.CH

```

I have used edgeR in R to calculate differential expression for contrasts involving 3 reps at one treatment combination to 3 reps at another treatment combination, for example

```
> y <- DGEList(counts = KCraw.data[keep,2:25], group = Group)
> y <- calcNormFactors(y)
>
> TMM <- KCraw.data[keep,2:25]
> for (i in 1:24) {
+   TMM[,i] <- TMM[,i] / (y$samples$lib.size[i] * y$samples$norm.factors[i]) * 1e6
+ }
>
> y <- DGEList(counts = TMM,group = Group)
>
> design <- model.matrix(~0+Group)
> colnames(design) <- levels(Group)
> y <- calcNormFactors(y,method = "TMM")
> y <- estimateDisp(y,design)
> fitQL <- glmQLFit(y,design)
> fit <- glmFit(y,design)
> myKC.contrasts <- makeContrasts(
+   H2H.H2L = FH.2DIC.CH - FH.2DIC.CL,
+   L2H.L2L = FL.2DIC.CH - FL.2DIC.CL,
+   H2H.L2H = FH.2DIC.CH - FL.2DIC.CH,
+   H2L.L2L = FH.2DIC.CL - FL.2DIC.CL,
+   H5H.H5L = FH.5DIC.CH - FH.5DIC.CL,
+   L5H.L5L = FL.5DIC.CH - FL.5DIC.CL,
+   H5H.L5H = FH.5DIC.CH - FL.5DIC.CH,
+   H5L.L5L = FH.5DIC.CL - FL.5DIC.CL,
+   H2H.L2L = FH.2DIC.CH - FL.2DIC.CL,
+   H5H.L5L = FH.5DIC.CH - FL.5DIC.CL,
+   H5L.H2L = FH.5DIC.CL - FH.2DIC.CL,
+   H5H.H2H = FH.5DIC.CH - FH.2DIC.CH,
+   L5L.L2L = FL.5DIC.CL - FL.2DIC.CL,
+   L5H.L2H = FL.5DIC.CH - FL.2DIC.CH,
+   levels=design)
> design
   FH.2DIC.CH FH.2DIC.CL FH.5DIC.CH FH.5DIC.CL FL.2DIC.CH FL.2DIC.CL FL.5DIC.CH FL.5DIC.CL
1           0          1          0          0          0          0          0          0
2           0          1          0          0          0          0          0          0
3           0          1          0          0          0          0          0          0
4           1          0          0          0          0          0          0          0
5           1          0          0          0          0          0          0          0
6           1          0          0          0          0          0          0          0
7           0          0          0          0          0          1          0          0
8           0          0          0          0          0          1          0          0
9           0          0          0          0          0          1          0          0
10          0          0          0          0          1          0          0          0
11          0          0          0          0          1          0          0          0
12          0          0          0          0          1          0          0          0
13          0          0          0          1          0          0          0          0
14          0          0          0          1          0          0          0          0
15          0          0          0          1          0          0          0          0
16          0          0          1          0          0          0          0          0
17          0          0          1          0          0          0          0          0
18          0          0          1          0          0          0          0          0
19          0          0          0          0          0          0          0          1
20          0          0          0          0          0          0          0          1
21          0          0          0          0          0          0          0          1
22          0          0          0          0          0          0          1          0
23          0          0          0          0          0          0          1          0
24          0          0          0          0          0          0          1          0
attr(,"assign")
[1] 1 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Group
[1] "contr.treatment"

> myKC.contrasts
            Contrasts
Levels       H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
  FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
  FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
  FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
  FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
  FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
  FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
  FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
  FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
            Contrasts
Levels       L5H.L2H
  FH.2DIC.CH       0
  FH.2DIC.CL       0
  FH.5DIC.CH       0
  FH.5DIC.CL       0
  FL.2DIC.CH      -1
  FL.2DIC.CL       0
  FL.5DIC.CH       1
  FL.5DIC.CL       0

```

After analyzing these contrasts, I wanted to estimate some sort of simple effect, such as the culture media nitrogen level. To do this, I ran the following code.

```
> myKC.contrasts <- cbind(myKC.contrasts,
+                         Development = c(1,1,-1,-1,1,1,-1,-1),
+                         FieldN = c(1,1,1,1,-1,-1,-1,-1),
+                         CultureN = c(1,-1,1,-1,1,-1,1,-1)
+ )
> myKC.contrasts
           H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
           L5H.L2H Development FieldN CultureN
FH.2DIC.CH       0           1      1        1
FH.2DIC.CL       0           1      1       -1
FH.5DIC.CH       0          -1      1        1
FH.5DIC.CL       0          -1      1       -1
FL.2DIC.CH      -1           1     -1        1
FL.2DIC.CL       0           1     -1       -1
FL.5DIC.CH       1          -1     -1        1
FL.5DIC.CL       0          -1     -1       -1

```

Once I rerun the analysis for the CultureN contrast and look at the result for a particular gene, I see that it's estimated log2FC is equal to the sum of every simple contrast.

```
> lrt <- glmQLFTest(fitQL, contrast=myKC.contrasts[,"CultureN"])
> topTags(lrt,n=nrow(y$counts))["GRMZM2G445575",]
Coefficient:  1*FH.2DIC.CH -1*FH.2DIC.CL 1*FH.5DIC.CH -1*FH.5DIC.CL 1*FL.2DIC.CH -1*FL.2DIC.CL 1*FL.5DIC.CH -1*FL.5DIC.CL
                 logFC   logCPM        F       PValue          FDR
GRMZM2G445575 -6.63617 5.417106 151.5261 3.691525e-11 2.825777e-08
# FC is a data frame of the logFC of each constrast in columns for each gene in rows
> sum(FC["GRMZM2G445575",c("H2H.H2L","L2H.L2L","H5H.H5L","L5H.L5L")])
[1] -6.636197

```

My first question is if this analysis is a valid way of summarizing the simple effects of each treatment. I would like to be able to also include the effects of the H2H.L2L and H5H.L5L contrast in the FieldN and CultureN comparison, but I am not sure how to do this, or if this would be valid because each of these contrasts includes treatments that have different levels of two treatment factors.

# Answer
I think the problem is in the design. There is no room for error or variation from a common ground.

I'm not sure of the solution, but I think you'll need to delete a column to give the design a degree of freedom, and probably add the intercept. You'll need to adjust the contrasts accordingly.

However, I would recommend to ask this at support.bioconductor.org. There are more experts on linear modelling and contrasts there than me. (If you ask there make it easier to copy and paste your code)
> 0 votes

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


# Question
Title: Finding orthologues using BLAST on the NCBI database
I'm an informatics student who has essentially zero knowledge of biology. I BLASTed my gene and have 1000s of results with very low E values. Where do I go from here if I want to find orthologues?

# Answer
> 7 votes
BLAST is often a sensible way to find gene/protein homologs, but resolving those as orthologs vs paralogs etc. is a non-trivial task, especially when considering large numbers of species.

Since you are new to bioinformatics I recommend you search an ortholog database that already exists.

Here are some databases:

All of these should allow you to search either with text or a sequence to identify the orthologous groups that your input belongs to at different taxonomic levels. Sequence searches either follow pairwise alignment (e.g. BLAST, DIAMOND, RAPSEARCH) or profile HMM type (e.g. HMMER3) approaches.

A more complete list of ortholog databases is available from "Quest for orthologs".

# Answer
> 0 votes
Two solutions based around the same premise:

1. Repeat the blast using "blastp". I suspect you have used blastn (nucleotides). To translate the gene Expasy is okay and find the common "translate", here.
2. Blasted via Blastx without any modification to the gene.

Blast is capable of doing just about anything if you are aware of all its parameters.

The reasoning is that the gene you have will be under-represented in the database. By switching to amino acids there is a much greater chance of finding orthologous hits.

# Answer
> 0 votes
Without jumping deep into the rabbit hole of orthologology, you can just use biomart

An instruction is here

---
Tags: blast, orthologues
---


# Question
Title: Scatter tasks around predefined batches inside a WDL workflow
Suppose I have two arrays Array\[File\] X and Array\[String\] Y both of length N, with Y being a sort of R-like factor that defines batch membership for X, for a total of M batches with M \< N. A valid example would be:

```
Array[File] X = ["file1", "file2", "file3", "file4", "file5"]
Array[String] Y  =  ["A",     "A",     "B",     "A",     "B"]

```

In which case N=5 and M=2. I want to scatter M tasks so that each task takes only the files corresponding to the same batch. In the previous example, one task would take as input:

```
Array[File] X =["file1", "file2", "file4"]
String Y = "A"

```

And the other task would take as input:

```
Array[File] X =["file3", "file5"]
String Y = "B"

```

Is there an easy way to do this within a WDL workflow in the general case? I understand that there are less elegant ways to do this as its own task but, as this seems like it could be a general purpose scatter/gather problem, an easy solution might be of general interest to the community.

# Answer
Here's an alternative to using a custom task, albeit still rather unappetizing:

```
workflow groupby {
  input {
    File? nil
  }

  Array[File] X = ["file1", "file2", "file3", "file4", "file5"]
  Array[String] Y  =  ["A",     "A",     "B",     "A",     "B"]

  scatter (y in ["A", "B"]) {
    scatter (p in zip(X,Y)) {
      File? x_or_nil = if (p.right == y) then p.left else nil
    }
    Pair[String,Array[File]] group = (y, select_all(x_or_nil))
  }

  output {
    Array[Pair[String,Array[File]]] ans = group
  }
}

```

(It'll be slightly cleaner with the upcoming `None` literal instead of dummy `nil`)

FWIW, here's a lengthy discussion about the idea of allowing developers to define custom library functions which just run a small task that manipulates data structures (hinting that it needn't necessarily offload to a worker node, with the huge latency that typically entails).
> 1 votes

# Answer
It turns out the development version has a function `collect_by_key` that is particularly helpful for scattering around predefined batches:

```
Array[String] X  =  ["A",     "A",     "B",     "A",     "B"]
Array[File] Y = ["file1", "file2", "file3", "file4", "file5"]
scatter (pair in as_pairs(collect_by_key(zip(X,Y)))) {
  call test {input: X = pair.left, Y = pair.right}
}

```

This will call task `t` first with input:

```
String X = "A"
Array[File] Y =["file1", "file2", "file4"]

```

And then with input:

```
String X = "B"
Array[File] Y =["file3", "file5"]

```

This is a readable solution to this problem allowed by `collect_by_key`. This function is not very well advertised (I had completely missed it) but kudos to who introduced it in the language. The only drawback is that `collect_by_key` is not available in version 1.0.
> 0 votes

---
Tags: wdl
---


# Question
Title: Detecting differentially expressed genes with foldchange >= 2 and FDR < 0.05
I'm using edgeR for differential analysis. Using glmTreat function I'm detecting differentially expressed genes between Tumor and Normal. I have set the arguments like below:

```
tr <- glmTreat(fit, contrast=contrast.matrix, lfc=2)
tab2 <- topTags(tr,n=Inf, adjust.method = "BH")
keep <- tab2$table$FDR <= 0.05

```

This gave me differentially expressed genes with logFC \> 2 at FDR \< 0.05 and genes with logFC \< 2 at FDR \< 0.05. \[There are few genes with logFC 1.9, 1.89 also\]

I have also tried with lfc=log2(2).

```
tr <- glmTreat(fit, contrast=contrast.matrix, lfc=log2(2))
tab2 <- topTags(tr,n=Inf, adjust.method = "BH")
keep <- tab2$table$FDR <= 0.05

```

To my surprise I see DE genes with logFC \>=1, FDR \< 0.05 and also logFC \<=1 at FDR \< 0.05.

Here is an example of DE genes I also found when I used lfc = log2(2)

```
GeneSymbols logFC   unshrunk.logFC  logCPM  PValue  FDR
DPYD-IT1    1.14493675  1.560389373 0.67468395  0.002972841 0.018704644
LINC00945   1.144288473 1.525191847 0.693845996 0.00651175  0.038064312
ITPKB-AS1   1.14046752  1.484991704 0.703887505 0.006620218 0.038585862
AC105206.1  1.122495266 1.4871905   0.699167758 0.006545299 0.038194633
LINC02066   1.120660287 1.558458881 0.66016102  0.004344743 0.026424298
AC090985.1  1.10235045  1.58968763  0.64343419  0.002745089 0.017426607
KCNJ6-AS1   1.100399688 1.717028637 0.613324799 0.002224039 0.01438451
LINC01494   1.090270138 1.490899364 0.676937984 0.006770608 0.039348027
AC099541.1  1.080957092 1.53817727  0.651356989 0.008449298 0.048356352
AL133396.2  1.012984207 1.512530231 0.63130643  0.008090259 0.046411928
AC011363.1  1.01207753  1.499708945 0.635740622 0.007226791 0.041797256

```

I'm very confused with these foldchange cutoffs.

Which one should I use to select differentially expressed genes based on fold change \>=2 and FDR \< 0.05?

# Answer
> 3 votes
Fold-change \>= 2 is the same as logFC (`log2(fold-change)`) \>= 1, so your example is doing exactly what you want. logFC is generally easier to think in and work with than fold-change, since (aside from computational reasons) then increases and decreases in expression differ only in sign and are, thus, easier to compare in magnitude (e.g., it's easier to tell that logFCs of 1.5 and -1.5 are the same magnitude change than it is that 3 and 0.333 are the same magnitude change in the fold-change).

For prioritizing the results, typically you would then rank these filtered results by fold-change and relevance (either given background knowledge or for genes with known function).

# Answer
> 3 votes
The edgeR authors recommend that you use a relatively low logFC threshold for `glmTreat` such as `lfc=log2(1.2)`. A `lfc` value as high as `lfc=log2(2)` is seldom required and only would only be appropriate for datasets with really large amounts of DE. The first threshold you used of `lfc=2` is equivalent to `lfc=log2(4)` and is too high for any dataset.

There's a few subtleties here:

1. First, there's the difference between logFC and FC as explained by Devon.
2. The `lfc` in `glmTreat` is not a fold-change cutoff! It is instead a testing threshold. All the significant genes will have actual logFCs well above the threshold.
3. The lfc threshold applies to the raw fold-changes (`unshrunk.logFC` in the output), not to the shrunk versions (`logFC` in the table).
4. The intent of `glmTreat` is to give a better gene ranking than would be provided by FDR or fold-change, so please don't re-sort the results by fold-change -- that would tend to defeat the purpose of the method.

I understand the documentation for `glmTreat` is a bit light, but you can get an idea of how to use it from our workflows.

I understand the simple fold-change cutoff of FC \> 2 with FDR \< 0.05 is often used in the literature, but the edgeR authors do not recommend it and hence we don't provide an automatic way to get it from the edgeR pipeline. For RNA-seq, fold-change cutoffs tend to prioritize genes with low expression and low counts. The `glmTreat` approach is more subtle but we think gives better results, both in terms of biological relevance and statistical correctness.

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


# Question
Title: Do I need to remove water molecules before protein-ligand docking?
I am going to screen small molecules against a protein target. Do I need to remove all the solvent molecules before performing VLS?

# Answer
> 3 votes
Yes, if you are using implicit solvent you ought to remove all of them.

Exception. Unless the crystallographic waters are trapped essential structural components of your active site and your ligand is not going to replace it, in which case keep it because its remove changes the energy of the active site —say you have a catalytic triad, with a water that is pinned between the acid (e.g. aspartate) and the base histidine, its removal would change the tautomerisation of latter and result in the nucleophile being protonated (and not nucleophilic). Having said that, if you don't have time to look at your structure (or cannot because of automation) and are okay with some artefacts from incorrect protonation strip them.

If you are not using implicit solvent, then it depends on the tool used —most require them stripped and re-added. If you do keep them, make sure that the crystallographic water (`HOH`) has the same residue code as the added water (`TP3`, `WAT` etc.).

---
Tags: proteins, docking
---


# Question
Title: Plotting metadata from GenomicRanges object
I have a GRanges object with metadata that I want to plot. Many of the online tutorials which manipulate these GRanges for plotting show how to plot the main columns (chromosome,start/end,strand) but nothing to plot the metadata as a function of the position in each row.

How precisely should I do this? Ideally, I would like to plot metadata from multiple GRanges objects on the same plot to compare them.

EDIT (following comment):

```
GRanges object with 6 ranges and 140 metadata columns:
      seqnames        ranges strand |
         <Rle>     <IRanges>  <Rle> |
  [1]     chr1 100000-199999      * |
  [2]     chr1 400000-499999      * |
  [3]     chr1 600000-699999      * |
  [4]     chr1 700000-799999      * |
  [5]     chr1 800000-899999      * |
  [6]     chr1 900000-999999      * |
                  eigen   domain         arm
              <numeric> <factor> <character>
  [1] -1.62475825405065     open          1p
  [2] -1.62475825405065     open          1p
  [3] -1.62475825405065     open          1p
  [4] -1.62475825405065     open          1p
  [5] -1.63354915136132     open          1p
  [6]  -1.6322733358426     open          1p
      exon.fraction       gc     short      long
          <numeric> <matrix> <numeric> <numeric>
  [1]       0.02968  0.45815         0         0
  [2]       0.04305  0.40192         0         0
  [3]       0.00959  0.42091         0         0
  [4]       0.04593  0.43337         0         1
  [5]       0.02425  0.44822         2         1
  [6]       0.16977  0.62123        27        56
         nfrags             ratio   short.corrected
      <numeric>         <numeric>         <numeric>
  [1]         0               NaN              <NA>
  [2]         0               NaN              <NA>
  [3]         0               NaN              <NA>
  [4]         1                 0 -3.50822623240367
  [5]         3                 2  2.81796238367926
  [6]        83 0.482142857142857  39.8824095079959
...

```

My question regards plotting a column like `eigen` as a function of the position (given by `seqnames` and `ranges`)

# Answer
> 1 votes
You can access the data.frame just using `$`. All of them can be accessed using `values()` . Something like this:

```
library(GenomicRanges)
gr = GRanges(seqnames=1,IRanges(start=1:10,width=1))
values(gr) = data.frame(eigen=rnorm(10),ratio=runif(10))

plot(start(gr),gr$eigen)

```

or

```
plot(start(gr),values(gr)$eigen)

```

---
Tags: r, bioconductor
---


# Question
Title: How to make work programs from the $PATH?
I am trying to analyze my RNAseq reads for defective genomes and I use this program (http://www.di-tector.cyame.eu/) that is suggested by Beauclair et al (https://pubmed.ncbi.nlm.nih.gov/30012569/) for that.

After downloading a .py script and running it with the following command:

```
sudo python3 DI-tector_06.py /mnt/e/nastya/SLX066-02/sequence.fasta /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz

```

It says at the beginning that

```
=================================
Program: DI-tector
Version: 0.6 - Last modified 25/05/2018
=================================
Requirement: (must be in your $PATH)
-bwa
-samtools

Optional: (must be in your $PATH)
-bedtools

```

I have tried to install bwa using:

```
sudo apt-get update -y
sudo apt-get install -y bwa

```

and samtools using

```
cd /usr/local/bin
sudo wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
sudo tar -vxjf samtools-1.9.tar.bz2
cd samtools-1.9
sudo make

```

I did not get exactly what does that mean: **(must be in your $PATH)** but I have searched here and there and discovered that it might mean `/usr/local/bin`

The program starts and, after some calculations it gives me the following error:

```
Input file:               /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz
Host reference:           None
Virus reference:          /mnt/e/nastya/SLX066-02/sequence.fasta
Remove segment:         < 15 nt
Remove reads with MAPQ: < 25
Allow InDel. length:    > 1 nt

=================================
Reading gziped input file...
=================================
Number of reads in input file: 18,035,544
=================================

=================================
Step 1/5 : Alignment against Viral reference...Filtering...
=================================
[E::bwa_idx_load_from_disk] fail to locate the index files
Traceback (most recent call last):
  File "DI-tector_06.py", line 183, in <module>
    Popen(["samtools", "view", "-Sb", "-F", "4"], stdin=out_file_onVirus, stdout=out_file_Virus).communicate()[0]
  File "/usr/lib/python3.6/subprocess.py", line 709, in __init__
    restore_signals, start_new_session)
  File "/usr/lib/python3.6/subprocess.py", line 1344, in _execute_child
    raise child_exception_type(errno_num, err_msg, err_filename)
PermissionError: [Errno 13] Permission denied: 'samtools'

```

I have no idea how samtools must be accessed, I tried `sudo` command with python3 but it did not work.

# Answer
`echo $PATH` will show the directories in your search path. Either copy the samtools executable there or add its location to the $PATH variable. You would be well served with a unix command line course or tutorial.
> 3 votes

# Answer
you can add a directory to `$PATH` by executing following command in the terminal:

`export PATH="$HOME/Documents/Software/bbmap:$PATH"`

In this example, all executable commands in the folder `bbmap` can be executed from any other directory in the terminal.

So you can install `samtools` in an appropriate directory and then execute `export PATH="$HOME/your/directory:$PATH"`. Then check if it succeeded by running `echo $PATH`. This latter command will return a list of directories which are assigned to `$PATH` and your new directory (where `samtools` is located) should be included.
> 2 votes

# Answer
You can run a piece of software by specifying exactly where it is, like `/home/usr/software/samtools_1.1.0/samtools`

Or if that location is in your $PATH, you can run it with just `samtools`. The software you downloaded needs to run those programs, but there's no way for you to give it the specific path for them; it's going to try and run them with just `samtools`. So you need to make sure that will work, by putting their locations in your PATH.

It is not clear to me that running every command with sudo is wise. If you don't understand what $PATH is, I'm pretty sure you should not be running everything with sudo.
> 1 votes

# Answer
The clean way to do this is to make your own `bin` within your user directory

```
cd ~
mkdir bin
vi .profile # or vi .profile_bash or for OS X using `csh` or `tcsh` vi .cshrc

```

For bash (Linux, or OS X) within vi type "i" then type,

```
PATH=$PATH:~/bin

```

or

```
PATH=~/bin:$PATH

```

... followed by `esc` key and "wp" `enter` Then enter for `csh` or `tcsh` e.g. for OS X

```
set path = ($path /Users/username/bin /Users/username/anaconda3/bin)

```

Same instructions for using `vi` as for `bash`. Dump bedtools into any directory (not the desktop) and connect the programs into the bin using `ln`

```
cd ~/progs/bedtools
ln -si $PWD/prog ~/bin

```

If you are using `csh` or `tcsh` you need to type `rehash` but `bash` does this automatically.

The programs will then be cleaning connect to your $PATH via `ln` without massive memory usage and visible anywhere within your home space. Its good practice for any Linux / OS X (time gone by UNIX) program or script and whilst you don't need it for this operation as you add more stuff its better.
> 1 votes

---
Tags: python, samtools, linux
---


# Question
Title: How can I get unmatched reads for defective genomes analysis using bwa and samtools?
I am trying to follow the Beauclair et al. paper (free version here https://rnajournal.cshlp.org/content/24/10/1285.long) for identifying defective genomes using their DI-tector program.

According to figure 2

Before using the program itself, it is highly recommended to filter my raw read using **bwa** and then **samtools**

Unfortunately I have no idea where to start - I have browsed official bwa and samtools sites, I have installed them using the following commands:

```
sudo apt-get update -y
sudo apt-get install -y bwa

```

I have indexed my reference genome with the following command (as I read it is required before the analysis)

```
bwa index sequence.fasta

```

Then, I have installed samtools

```
cd ..
wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
tar -vxjf samtools-1.9.tar.bz2
cd samtools-1.9
make
sudo apt install samtools

```

I wonder what will be the further steps to get these "UNMATCHED READS at the end of step i", since Di-tector program takes virus\_reference genome (I have it) and these unmatched reads.

First of all I want to align my .fastq files vs `GCF_000001405.39_GRCh38.p13_genomic.fna` \- what will be the output file? If I got it right - I need this output file to filter with samtools, am I correct? If I am not mistaken, the output after samtools will be the .sam file that I will have to put into a DI-tector - correct?

# Answer
If you want to keep track of the unmapped reads, I recommend you to use Bowtie2 instead of bwa. Start by building an index:

```
bowtie2-build genome.fa genome.fa

```

Then map your reads, using arguments to save unmapped reads in a specific file like this:

```
bowtie2 -x genome.fa -1 reads_1.fq -2 reads_2.fq -S bowtie2_mapping.sam --un-conc-gz unmappedreads_%.fq.gz

```

It will write the alignment as well as unmappedreads\_1.fq.gz and unmappedreads\_2.fq.gz that you may use in DI-tector.
> 1 votes

---
Tags: samtools, bwa
---


# Question
Title: Encoding a 15bp DNA sequence into a shorter than 15 digits number
I have a series of 15bp DNA sequences from single-cell barcodes and I would like to encode them into a shorter than 15 digits number so that each unique 15bp DNA sequence will create a unique numeric value. What's the best way of achieving this?

# Answer
What you are asking is called a "hash" (wiki).

In python, get your oligo as a string and use the standard library method `hash(oligo)`, this is a number that gets used for a lot of operations. It is however always a 18+1 digit number in decimal, so a different solution is needed.

There are four nucleobases, so binary conversion from scratch is easy and best option and each solution will be unique and small. Each base takes up 2 bits and you concatenate these and convert them to hex or decimal —the length should be 10 decimal digits. One issue is when you have a oligo of a different length which start with adenosine or not.

One issue with this is the order. I am not sure where it comes from but the bases are nearly often expressed in the order ATGC. So...

```
mapping = {'A': '00', 'T': '01', 'G': '10', 'C': '11'}

```

It is satisfying how the first digit is the W/S, while the second is R/Y.

Now simply covert and add them `compressed = functools.reduce(lambda acc,base : acc+mapping[base], oligo, '')` and covert the string to decimal like so `int(compressed, 2)`.
> 3 votes

---
Tags: phylogenetics
---


title: Error pasing the following line in pdb

text:
# Question
Title: Error pasing the following line in pdb
I converted an .sdf file containing around 50,000 anti-viral agents to individual .pdb files using OpenBabel. However, when I tried opening the .pdb files with Autodock I kept getting the error below. What do you think is causing the error.

Thank you for your time and patience.

**Error Message:**

> “ Error pasing the following line in pdb: HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000 ."

# Answer
> 1 votes
A PDB file has strict spacing rules as documented here. The coordinates should be:

```
31 - 38       Real(8.3)      x             Orthogonal coordinates for X.
39 - 46       Real(8.3)      y             Orthogonal coordinates for Y.
47 - 54       Real(8.3)      z             Orthogonal coordinates for Z.

```

So trying this on your line:

```
'HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000'[30:39]

```

I get `9.051 0.`, which is not a number.

The whole line is formatted incorrectly however. So as you say you are using Autodock why not pass a `mol2` to `prepare_ligand4.py` instead which is an acceptable format?

Alternatively, you could convert your multi-entry sdf to pdb files via RDKit in Python:

```
suppl = Chem.SDMolSupplier(sdfile)
for mol in suppl:
      # mol = Chem.AddHs(mol) # has hydrogens?
      # AllChem.EmbedMolecule(mol) # is 3D?
      AllChem.ComputeGasteigerCharges(m) # sdfs do not have partial charges
      Chem.MolToPDBFile(m, f"{mol.GetProp('_Name')}.pdb")

```

---
Tags: pdb, docking, errors
---


# Question
Title: Removing stop codons from alignment file for dn/ds calculation
This is one of the question which is very similar to my problem or almost same to my problem

I ran with dummy data nucleotide alignment and the amino file it works perfectly fine,i can run the dn/ds from phangorn without any issue.

But when i ran with my data nucleotide alignment and the amino file I do get this warning

```
In kaks(carps) : sequence lengths are not a multiple of 3

```

and error.

The code

Any help or suggestion would be highly appreciated.

# Answer
I would align your sequences with the software available here: https://github.com/veg/hyphy-analyses/tree/master/codon-msa
> 2 votes

---
Tags: phylogenetics, perl
---


title: How to anonymise a bam file

text:
# Question
Title: How to anonymise a bam file
I have a bam file of, for instance, RNA-Seq, which contains patient-identifiable data in the form of Single Nucleotide Polymorphisms (SNPs) throughout.

I would like a technique to take this aligned bam file and replace all the non-reference bases with the reference bases at the equivalent positions (hard-clipping any soft-clipped bases while doing so), in order to have an anonymised bam file without the patient-identifiable information.

What is the best way to do this? I want the resulting anonymised bam file to have the reads in the same positions, for downstream tools that care about the alignment positions.

# Answer
You can simply replace the sequence and quality scores with a single `*` in each field. The CIGAR information will still be there, so it will be apparent where there were InDels and possibly mismatches (assuming your aligning uses `=` and `X` CIGAR operations), but since the actual sequence of these variants cannot be known this is quite anonymized.
> 1 votes

---
Tags: sam, ngs
---


# Question
Title: What does instantaneous rate matrix mean?
PhyML is a tool used for maximum likelihood estimation analysis. I was running a PhyML analysis to check the rate of nucleotide change in DNA sequences and in the PhyML output, there are rate parameters and an instantaneous rate matrix. Could you please clarify the difference between the two?

What I've understood so far is that the rate parameters indicate, for ex: that the A -C transition rate is 1.95 and so on. If so what does the instantaneous rate matrix tell us?

# Answer
> 1 votes
**Reversible** At its simplist the GTR rate, is the General Time Reversible model and most importantly infers the matrix is "reversible", thus as many A mutations will go to C as C mutations go to A and the number of mutations is relative to the G-T mutation rate. The G-T mutation rate is one of the least common mutations (transversion) so all other values are above 1. If you plotted this as a matrix the top right-hand matrix is identical to the bottom left hand matrix because its reversible, i.e. A \<-\> C . This is cool because you don't need to tell the algorithm which sequence is the root.

**Directional** The "instantaneous" matrix is a directional matrix, thus A -\> C is not the same as C -\> A. In order to obtain an accurate directional matrix you need to assign a root, its like lets start from the deepest node in the tree (the one between the root/outgroup and all ingroup sequences) and now - if that is correct - we can remove the restrictions of reversibility and impose directionality on the mutation process. This is needed if you are performing molecular clock type calculations. However, the key is assigning the root/outgroup - if it assigns it e.g. by mid-point rooting that is not necessarily a good idea.

**Advice.** My advice is if you are building a basic tree and haven't assigned a root/outgroup, obviously using GTR which is a good move, then you don't to worry about the intananeous matrix because it might not be meaningful. If you are running a molecular clock/dating/rate estimate type model and a very clear what the root/outgroup is and that is specified in the input, then the instananeous matrix becomes important.

---

The mutation rates of GTR and instananeous are not directly comparable because the GTR rate is with respect to G \<-\> T mutation, whilst instananeous is primarily concerned with the directionality of mutation and is a 12-parameter model, as opposed to GTR which 6-perameters, where one of those parameters is the reference.

Thus G-T is reversible in GTR (the reference) and so e.g. A\<-\>C is reversible with respect to G\<-\>T. In instananeous matrix G -\> T and T -\> G are two different rates as are A-\> C and C -\> A, so one approximation based on an assumption with respect to a second approximation based on the same assumption can't directly equate to a mutation rate which removes that assumption.

---
Tags: phylogenetics, phylogeny, hidden-markov-models, nucleotide-models
---


title: Can't align 2 sequences, (MAFFT killed)

text:
# Question
Title: Can't align 2 sequences, (MAFFT killed)
I have 2 .fa sequenes like:

```
>some description of first...
ACTGACTACGTACGACTACGATCTGACTACACGTAGATAGACTAGTCACTACG
ACTGACTACGTACGACTGACTACATGACTGAC...
> some description of second...
ATCGATAGCTACGATCGATAGTAGATAG....

```

This file is 461mb (1 chromosome of homo sapiens & 1 chromosome of pan troglodytes) and I want to align it. I am using mafft.

```
mafft --globalpair --maxiterate 1000 Homo_sapiens_1_Pan_troglodytes_1.fa > ./mafft_results/Homo_sapiens_1_Pan_troglodytes_1_output.fa

```

And it returns me

```
generating a scoring matrix for nucleotide (dist=200) ... done
All-to-all alignment.
/usr/bin/mafft: line 2719:  6030 Killed                  "$prefix/tbfast" _ -u $unalignlevel $localparam -C $numthreads $seqtype $model -g $pgexp -f $pggop -Q $spfactor -h $pgaof -A $usenaivepairscore $focusarg _ -+ $iterate -W $minimumweight -V "-"$gopdist -s $unalignlevel $legacygapopt $mergearg $outnum $addarg $add2ndhalfarg -C $numthreadstb $rnaopt $weightopt $treeinopt $treeoutopt $distoutopt $seqtype $model -f "-"$gop -Q $spfactor -h $aof $param_fft $localparam $algopt $treealg $scoreoutarg $focusarg < infile > /dev/null 2>> "$progressfile"

```

Why? How could I align it for further comparing?

**EDIT 1**

MAFFT working cool, when I tried to use it on test sequence file

```
>test1
ATCGATCGATGACTACGATCGACTACGTACG
ACTAGCTAGCTACGCATAGTAGATCGA
ATCGACTAGCTCAGATCGACTAGATCGATCGT
>test2
TTATCCACTGATCGACTCAGATCGACTCAGAC
ATCGATCAGCTACGATGACTACGTAC
ACTGACTAGCTAGACTGACTAGCTACGTCA

```

and after MAFFT `mafft test.fa > output.fa` it returned

> test1 ------atcgatcgat-------------gactacgatcgactacgtacgactagctagc tacgcatagtagatcgaatcgactagctcagatcgactagatcgatcgt test2 ttatccactgatcgactcagatcgactcagacatcgatcagctac---------gatgac tacgtac----------actgactagct-agactgactagctacgtca-

So it works.

I can use all 4 at the example. Now about my large file:

```
alex@alex:~/Desktop/test/BOTH/homo_sapiens_pan_troglodytes$ mafft --thread 4 --maxiterate 2 Homo_sapiens_1_Pan_troglodytes_1.fa > ./mafft_results/Homo_sapiens_1_Pan_troglodytes_1_output.fa
nthread = 4
nthreadpair = 4
nthreadtb = 4
ppenalty_ex = 0
reallocating...
done.
reallocating...
done.
reallocating...
    done.
reallocating...
done.
reallocating...
...
stacksize: 8192 kb
reallocating...
done.
reallocating...
done.
reallocating...
done.
...
generating a scoring matrix for nucleotide (dist=200) ... done
Gap Penalty = -1.53, +0.00, +0.00

Making a distance matrix ..
    1 / 2 (thread    1)
done.

Constructing a UPGMA tree (efffree=1) ...
    0 / 2
done.

Progressive alignment 1/1...

len1=230480994, len2=223796828, Switching to the memsave mode
STEP     1 / 1 (thread    0)m/usr/bin/mafft: line 2719:  3192 Killed                  "$prefix/disttbfast" -q $npickup -E $cycledisttbfast -V "-"$gopdist -s $unalignlevel $legacygapopt $mergearg -W $tuplesize $termgapopt $outnum $addarg $add2ndhalfarg -C $numthreads-$numthreadstb $memopt $weightopt $treeinopt $treeoutopt $distoutopt $seqtype $model -g $gexp -f "-"$gop -Q $spfactor -h $aof $param_fft $algopt $treealg $scoreoutarg $anchoropt < infile > pre 2>> "$progressfile"

```

Doesnt matter I do not enter `--thread 4` or not, `--maxiterate 2` or `--maxiterate 1000`, it still shows me same error. It takes around 2 minute. Ram usage when running at the first minute is 6656mb, after 1 min 11028mb and then it out of ram and error..

*PS*

discussion moved here

# Answer
`--thread 4` is essential, `--maxiterate` just try and get it working, not getting it to work optimally. Just make sure there are no other processes running on the machine when --thread 4 is operating, e.g. you are using another machine as your desktop

The gist is as you spread the calculation across your processors (4) and simplify the calculation (lowering the iterations), you are getting further into the calculation, for example now you are getting into gap penalties.

If you are agreed then I would try performing the alignment using just one chromosome of the human genome against the chimp genome and I suspect that would work. Then take it from there.

The objective is to get the calculation working, for example I wouldn't trust a result with max. iterations of 2, but once you have it working you can then start refining the parameters towards optimal results. The i3 processor isn't perfect in context.

Using `--thread 4` is super important and you want to increase maxiterate towards the default if the machine supports it.

I suspect the solution is to perform the calculation on an 8+-core machine and set `--thread 8` (or higher depending on the number of cores on the machine), at least it will get further into the calculation. With 4-cores and an i3 (diagnostics of RAM usage would be nice) in any case it just doesn't seem to like it. The alternative is to perform this across the genome via sampling the chromosomes on your current i3.

---

Just to explain why this works. When you have a RAM intensive job and you spread it across cores the RAM overhead massively reduces for each core that is added. Don't ask me how this works in the hardware mecahanics of RAM distribution per core. The reason I know this is because of emperically parallelising Bayesian calculation using each chain on a different core.
> 1 votes

---
Tags: sequence-alignment
---


# Question
Title: Problem with merge data while trying to convert gene names in R
*This question has also been asked on Biostars and StackOverflow*

I've been trying to code (in R) a way to convert gene accession numbers to gene names (from RNAseq data). I've looked at all the related questions and tried to modify my code such, but for some reason it's still not working. Here is my code, where `charg` is a character vector of the gene accession ID's of the data set `resdata`:

```
charg <- resdata$genes
head(charg)

library(biomaRt)

ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")

theBM <- getBM(attributes='ensembl_gene_id','hgnc_symbol',
      filters = 'external_gene_name',
      values = charg,
      mart = ensembl)

resdata <- merge.data.frame(resdata, theBM, by.x="genes",by.y="ensembl_gene_id")

```

Here's some output (where I'm struggling):

```
> head(charg)
[1] "ENSG00000261150.2"  "ENSG00000164877.18" "ENSG00000120334.15"
[4] "ENSG00000100906.10" "ENSG00000182759.3"  "ENSG00000124145.6"

> dim(theBM)
[1] 0 1

> head(theBM)
[1] ensembl_gene_id
<0 rows> (or 0-length row.names)

> dim(resdata)
[1] 20381    11
> resdata <- merge.data.frame(resdata, theBM, by.x="genes",by.y="ensembl_gene_id")
> dim(resdata) #after merge
[1]  0 11 #isn't correct -- just row names! where'd my genes go?

```

Thank you.

# Answer
> 3 votes
This is the code to get a look-up table to convert between Ensembl ID and HGNC:

```
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
theBM <- getBM(attributes=c('ensembl_gene_id','hgnc_symbol'),
               filters = c('ensembl_gene_id'),
               values = gsub("\\..*", "", charg),
               mart = ensembl)

```

What Devon was posting is correct but misses a `c()` around the attributes values. For further help please provide the content of `resdata` which you should always do when posting a question, since we cannot read minds. `Does not work` by the way is not a proper error description.

Once you have the output do:

```
resdata$genes <- gsub("\\..*", "", resdata$genes)

merge(x = theBM,
      by.x = "ensembl_gene_id",
      y = resdata,
      by.y = "genes")

```

Note that I had to go to that SE crosspost to get the content of `resdata`, this is not how this goes. Please post all relevant data up front in the future otherwise your questions might get downvoted and closed. Please also avoid cross-posting. if you provide proper information you usually get a good answer in time.

Edit: Just realized you also cross-posted this to Biostars even twice. Please stop this. I closed the Biostars posts and gave my two cents on this behaviour over there.

# Answer
> 1 votes
Those aren't external\_gene\_name's, they're ensembl\_gene\_id\_versions:

```
theBM <- getBM(attributes='ensembl_gene_id','hgnc_symbol',
               filters = 'ensembl_gene_id_version',
               values = charg2,
               mart = ensembl)

```

Note that you'll get more hits if you strip the gene ID versions off:

```
charg2 = sapply(strsplit(charg, '.', fixed=T), function(x) x[1])
theBM = getBM(attributes='ensembl_gene_id','hgnc_symbol',
              filters = 'ensembl_gene_id',
              values = charg2,
              mart = ensembl)

```

---
Tags: r, rna-seq, gene, biomart
---


# Question
Title: How to extract unmatched reads using bwa and samtools?
I have a single read (**NOT paired**) that I need to pass through the workflow described in Beauclair et al. paper (free version here https://rnajournal.cshlp.org/content/24/10/1285.long) for identifying defective genomes using their DI-tector program.

Here in Materials and methods the procedure is described as the following:

*The first step of the workflow consists of an alignment of the reads against the host genome (Fig. 2i). This step aims at discarding reads that map to the host genome and may partially map to the viral genome after segmentation, and at reducing the working file size. For example, MV and rMV-ΔV data sets were generated from total RNA samples of infected cells and contained mostly reads mapping the human genome (≈99%). This step uses a combination of `bwa mem` and `samtools` view with the parameters `–bS –f4`. An additional step consists of an alignment of the reads against the viral genome of interest, in order to exclude reads perfectly mapping the viral genome. Therefore, only unmapped reads are further analyzed. Of note, clipped reads (i.e., CIGAR motif contain S or H) are also conserved. Some of these reads may map to viral genome recombination junctions that are present in DI genomes.*

I have already been suggested to use `bowtie2` instead of `bwa` but, firstly, the output is not clear for me and secondly, I would like to test the official protocol.

Since the article suggests to use bwa and samtools on this very first step, that is what I have done so far:

0. (Optional) not sure if this is important but as suggested by someone I have transformed .fna --\> .fa

```
cp GCF_000001405.39_GRCh38.p13_genomic.fna GCF_000001405.39_GRCh38.p13_genomic.fa

```
1. Indexed human genome with bwa

```
bwa index GCF_000001405.39_GRCh38.p13_genomic.fa

```
2. aligned my only read with indexed genome from step 1

```
bwa mem GCF_000001405.39_GRCh38.p13_genomic.fa /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz > whole.sam

```
3. Transform `.sam` --\> `.bam`

```
samtools view -S -b whole.sam -o whole.bam

```
4. Separated unmapped reads (as it is recommended in Materials and Methods using `-f4`)

```
samtools view -f4 whole.bam > sample.unmapped.sam

```
5. Converted unmapped reads into `.fastq` format (since this is the format used by the software later)

```
samtools fastq sample.unmapped.sam > unmatched.fastq

```

I wonder if these steps are correct?

# Answer
Yes, correct but overly many steps. There is no need to convert between bam and sam. `samtools` can read from `stdin` and handles both sam and bam and `samtools fastq` can interpret flags, therefore one can shorten this to:

```
bwa mem (...options) | samtools view -o out.bam
samtools fastq -f 4 out.bam > unmatched.fastq

```
> 5 votes

---
Tags: sequence-alignment, samtools, bwa, bowtie2
---


title: .SDF file to .PDB file

text:
# Question
Title: .SDF file to .PDB file
I am currently working on identifying effective antiviral compounds against the corona virus. I have downloaded the Covid-19 AntiViral Candidate Compounds Dataset from the www.cas.org webpage, the link to which I have provided below. Me and a few other doctors are using Docking software to predict the docking capacity of these anti-viral compounds. Docking software requires the anti-viral compounds to be in the .pdb file format only then is it recognized by Automated docking software. **Is there any way I can convert the .sdf file containing the anti-viral compounds to .pdb file or is there a way of obtaining the .pdb files for list of anti-viral compounds**.

*Thank you for your time and patience, Anas Hussain.*

**Anti-Viral Candidate Compounds Dataset** https://www.cas.org/covid-19-antiviral-compounds-dataset

# Answer
> 1 votes
Use OpenBabel.

Something like `babel -isdf file_in.sdf -opdb file_out.pdb`.

# Answer
> 1 votes
> Bizarrely, Obabel seems to have written a PDB that was incorrectly spaced for the OP as seen in Error pasing the following line in pdb

If you can use python, `rdkit` is a great package for cheminformatics —best installed via conda though.

It allows you do all sorts of things, such as match substructures, embed conformers, edit compounds, compute advanced parameters and many forms of shape & colour finger-/foot-printings.

Given a multi-compound sdf you can read each individual molecule and write them out as PDBs.

```
from rdkit import Chem
with i, mol in enumerate(Chem.SDMolSupplier('in.sdf')):
    Chem.MolToPDBFile(f'out_{i}.pdb')

```

We can be fancier with the output name

```
from rdkit import Chem
import re
with mol in Chem.SDMolSupplier('in.sdf'):
    # The headerline of a mol file is retrieved oddly.
    name = mol.GetProp('_Name')
    # make filename safe
    slugify = lambda name: re.sub(r'[\W_.-]+', '-', name)
    Chem.MolToPDBFile(f'{slugify(name)}.pdb')

```

To specify what name (*e.g.* `CA`) each atom gets and what is the residue name, chain and residue position —information absent in the sdf, it gets a bit odd.

```
  names = [' CA ', ' CB ', ' CG ', ' CD '] # say
  for i, atom in enumerate(mol.GetAtoms()):
        info = Chem.AtomPDBResidueInfo(names[i], # or whatever.
                                       residueName='LIG')

```

For more see this blog post.

---
Tags: pdb, docking
---


title: Generating 3D coordinates error

text:
# Question
Title: Generating 3D coordinates error
I have been trying to convert a multi entry sdf file to pdb using openbabel. Everytime I check the Generate 3D coordinates box I end up with .ENT files instead of .PDB files. What am I doing wrong.

Thank you.

# Answer
An ent file is a pdb file. Just change the extension. PDBe for example provides their PDB files as such.
> 2 votes

---
Tags: docking, 3d-structure
---


# Question
Title: How do I set a neural network to loop multiple times and average the resulting values?
I have a script in R/RStudio which creates random datasets of binomial variables, feeds them through a neural network, and calculates their likelihood ratio statistic and deviance. I'd like the script to loop with the `seed` incremented by 1, and collect the value from each run into a data frame, from which the values get averaged. For example, if the script returns `5, 10, and 15` on the first, second, and third iterations of the loop, I'd like the script to create a data frame with those values and then calculate their mean.

Here is my script, annotated:

```
library(neuralnet)

x=1021

#change nrow to number of times code is looped
result <- data.frame(matrix(nrow = 1000, ncol = 2))
colnames(result) <- c("LogL", "D")

repeat{

set.seed(x)

#2 inputs, max number of terms 3
x1=rbinom(2000,1,0.5)
x2=rbinom(2000,1,0.5)
y=rbinom(2000,1,0.5)

data=data.frame(x1,x2,y)

head(data)

#hidden units = 2
nnet=neuralnet(y~x1+x2,data,hidden=2,err.fct = "ce",linear.output = F,likelihood = T)

predictions=compute(nnet,data[,1:2])

predictions$net.result

#likelihood ratio stat
argument=y*log(predictions$net.result)+(1-y)*log(1-predictions$net.result)

head(argument)

#sum of likelihood ratio stat
LogL=sum(argument)

LogL

D=-2*LogL

D

result[x-1022, 1] <- LogL
result[x-1022, 2] <- D

result

#results <- data.frame(LogL, D)

x=x+1
if (x==2021) {
  break
}

print(x)

}

result

```

My repeat loop works, but a thousand iterations, for example, takes multiple minutes. Is there a more efficient way to code such a loop?

# Answer
> 1 votes
The concern with your approach is called 'leakage' because you are parameterising the same data set that you are training. This can easily lead to overtraining is in ANN (artificial neural networks) and is not cool. Overtraining is where the ANN very tightly wraps around the training data set that it performs poorly on any other data set. In other words it is overtightened.

The way they get around it is to split the data set usually 60:20:20 into a training split, parameterisation (validation) split and test (testing) split. If the parameterisation split gets mixed with the training set - leakage happens and overtraining results. Normally the data is trained, and the trained model parameters are assessed/optimised in the parameter split. The trained data and optimised parameters are then tested in the test split.

As long as you are aware of the approach you are using, that I guess that is ok. If you separate your training and parameterisation (validation) data that would be cool.

---
Tags: r, statistics, machine-learning, loop
---


# Question
Title: RNA-Seq data transformation prior to sample correlation analysis
If I’m starting with Deseq2 normalized counts what are some preprocessing steps that I should apply to these data before estimating sample correlation using the cor function in R? For example, would it make sense to quantile normalize the data?

# Answer
You don't absolutely need to perform any further steps before calculating the correlation. However you would likely benefit from the following sort of workflow:

* Transform your data with vst() or robust log (see the DESeq2 vignette) so genes are on a more or less similar scale.
* Select the top 300-500 most variable genes
* Use this subset with to compute sample-wise correlation.

This also what follows what you would find in the DESeq2 user guide, which I strongly encourage you to read.
> 4 votes

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


# Question
Title: How to design DESeq2 LRT model with individuals nested in 2 levels?
We have a complicated experimental design that we would like to perform LRT analysis for. Our main goal is to discover significant genes for the "Injection:Social" interaction term across the entire dataset by removing it from the LRT reduced model, and as a bonus we are also interested in discovering significant genes for that interaction term for each respective brain region.

```
Sample  Injection   Social  Region  Individual  ind.n
HY06    L   ISO HY  S06 S1
NST6    L   ISO NS  S06 S1
TN06    L   ISO TN  S06 S1
HY08    L   ISO HY  S08 S2
NST8    L   ISO NS  S08 S2
TN08    L   ISO TN  S08 S2
HY30    L   KF  HY  S30 S1
NST30   L   KF  NS  S30 S1
TN30    L   KF  TN  S30 S1
HY32    L   KF  HY  S32 S2
NST32   L   KF  NS  S32 S2
TN32    L   KF  TN  S32 S2
HY64    L   KFC HY  S64 S1
NST64   L   KFC NS  S64 S1
TN64    L   KFC TN  S64 S1
HY65    L   KFC HY  S65 S2
NST65   L   KFC NS  S65 S2
TN65    L   KFC TN  S65 S2
HY19    L   NF  HY  S19 S1
NST19   L   NF  NS  S19 S1
TN19    L   NF  TN  S19 S1
HY24    L   NF  HY  S24 S2
NST24   L   NF  NS  S24 S2
TN24    L   NF  TN  S24 S2
HY05    S   ISO HY  S05 S1
NST5    S   ISO NS  S05 S1
TN05    S   ISO TN  S05 S1
HY12    S   ISO HY  S12 S2
NST12   S   ISO NS  S12 S2
TN12    S   ISO TN  S12 S2
HY31    S   KF  HY  S31 S1
NST31   S   KF  NS  S31 S1
TN31    S   KF  TN  S31 S1
HY34    S   KF  HY  S34 S2
NST34   S   KF  NS  S34 S2
TN34    S   KF  TN  S34 S2
HY62    S   KFC HY  S62 S1
NST62   S   KFC NS  S62 S1
TN62    S   KFC TN  S62 S1
HY63    S   KFC HY  S63 S2
NST63   S   KFC NS  S63 S2
TN63    S   KFC TN  S63 S2
HY04    S   NF  HY  S04 S1
NST4    S   NF  NS  S04 S1
TN04    S   NF  TN  S04 S1
HY20    S   NF  HY  S20 S2
NST20   S   NF  NS  S20 S2
TN20    S   NF  TN  S20 S2

```

My first attempt was building simple full (m1) and reduced (m2) models that gets directly at our question of interest but doesn't control for nested individuals.

```
m1 <- model.matrix(~ Region + Social * Injection, colData_filt)
m2 <- model.matrix(~ Region + Social + Injection, colData_filt)

```

We want to control for individual/batch effects, which is nested within both "Injection" and "Social" but not region, as we have three brain regions per individual. I followed the example in the DESeq2 manual for creating a term (ind.n) distiguishing individuals nested within groups, but now I'm not sure how to create the full and reduced model given that I have one more level than the example.

I've tried a really elaborate full model (m1) with the interaction term of interest (Injection:Social) removed for the reduced model (m2), but I'm not sure this is correct based on our design.

```
m1 <- model.matrix(~ Injection + Injection:ind.n + Injection:Social + Injection:Region + Social + Social:ind.n + Social:Region + Region, colData_filt)
m2 <- model.matrix(~ Injection + Injection:ind.n + Injection:Region + Social + Social:ind.n + Social:Region + Region, colData_filt)

```

I'm assuming this is wrong, but even if this was by some miracle the correct formulation, would there be a way to extract genes that explain the "Injection:Social" interaction term for separate brain regions?

As a work-around, I subsetted the data by region and ran three separate LRT analyses for each subset and compared the results. While this simplified the model to look like the first example above, I worry that we lose some power by ignoring the fact we have multiple brain region samples from single individuals across the dataset.

Any guidance is much appreciated. Thanks in advance

# Answer
> 4 votes
From what I can gather, you want to account for the effect of individual, nested within region. That is, you want to see after accounting for these, is there a consistent effect for Injection:Social across all conditions.

So you set up the model like this:

```
m1 <- model.matrix(~ ind.n*Region + Injection + Social + Injection:Social,data=..)

```

The last term should be Injection:Region and you can just use the waldTest (default) in DESeq2 for this term.

What does the terms do? `ind.n*Region` is the equivalent of `ind.n + Region + ind.n:Region` , and with this you effectively get an effect for every region in every individual.

Why don't we need the `Injection:ind.n` or `Social:ind.n` or `Social:Region`. These terms indicate the effect of Injection or Social can vary by individuals or regions. Most likely introducing too many parameters when you are interested in a common effect. Also you do not have the replicates or samples to distinguish this effect from region or other effects.

Since you provided an example we can run DESeq2 and you can see how the results look like:

```
mat = counts(makeExampleDESeqDataSet(n=1000,m=48))
dds = DESeqDataSetFromMatrix(mat,df,~ ind.n*Region + Injection + Social + Injection:Social)

dds = DESeq(dds)

resultsNames(dds)
 [1] "Intercept"            "ind.n_S2_vs_S1"       "Region_NS_vs_HY"
 [4] "Region_TN_vs_HY"      "Injection_S_vs_L"     "Social_KF_vs_ISO"
 [7] "Social_KFC_vs_ISO"    "Social_NF_vs_ISO"     "ind.nS2.RegionNS"
[10] "ind.nS2.RegionTN"     "InjectionS.SocialKF"  "InjectionS.SocialKFC"
[13] "InjectionS.SocialNF"

```

The terms you need are "InjectionS.SocialKF","InjectionS.SocialKFC", "InjectionS.SocialNF", and you can look at each of them:

```
head(results(dds,name="InjectionS.SocialNF"))
log2 fold change (MLE): InjectionS.SocialNF
Wald test p-value: InjectionS.SocialNF
DataFrame with 6 rows and 6 columns
              baseMean     log2FoldChange             lfcSE               stat
             <numeric>          <numeric>         <numeric>          <numeric>
gene1  9.9811166787259   1.25304112986447 0.819806376919295   1.52845984752303
gene2 30.3449455820337 0.0329442893152027 0.705199688255367 0.0467162562092241
gene3 3.83223545055379    1.0281136369045  1.64095596190233  0.626533350543196
gene4  11.232305747171  0.595738624408923   0.8243883031544  0.722643227868976
gene5 6.70950627004097  0.756449993378065   1.0631622863378  0.711509430967263
gene6 26.1431134888287 -0.854784518963918 0.625714541243558  -1.36609342219393
                 pvalue              padj
              <numeric>         <numeric>
gene1 0.126398405431826 0.978671002658464
gene2 0.962739373909937 0.999897888026606
gene3 0.530965168963838 0.978671002658464
gene4 0.469899103018657 0.978671002658464
gene5 0.476768608734069 0.978671002658464
gene6 0.171909642630577 0.978671002658464

```

As mentioned, you can do a LRT if you want to test all the Injection:Social term interactions terms at one go, that is, the null hypothesis is that all of them are zero:

```
dds = nbinomLRT(dds,reduced=~ ind.n*Region + Injection + Social)
results(dds)

```

Usually the individual terms are more intuitively than testing all of them are zero, but you might have a special need for this.

---
Tags: deseq2
---


# Question
Title: What are the state-of-the-art cell-type RNA-Seq deconvolution methods?
I would like to find the proportion of each cell-type in bulk RNA-Seq transcriptomics data. I am looking for some guidance on the following:

* What are the state-of-the-art methods?
* What are their assumptions? (e.g. do they assume that cell-type expression signatures are available in advance?)

# Answer
I see a lot of people using xCell. There should also be papers systematically compare many different methods, which might interest you.
> 1 votes

---
Tags: rna-seq, gene-expression, machine-learning
---


# Question
Title: Do you know a program/script to count how many sequences have mismatches on each position of my primers?
I'm working with primer analysis and design, and so far I've performing it in Excel but: a) I don't like Excel, b) I do it manually, c) I don't like Excel, d) it gets tedious when working with thousands of sequences and tens of primers. My analysis consists of establishing how many sequences have mismatches on each position of my primers, and also count how many mismatches each sequence has with the corresponding primers (see attached image). This last analysis I've been able to (partially) solve it by using ncbi-blast+ tool. However, **I will find very helpful a program/script where you input your sequences and primers fasta files and retrieve the number of mismatches at every primer position. Does anyone know something that could help me achieve all this?** If it also has an output tab-separated file as the one attached that would be perfect for easy visualization of results. Suggestion/orientations are also grateful! I just want to know if there's something to help me out, otherwise I'll try writing a script to automatize my analysis. Thank you!

# Answer
> 1 votes
I'm amazed at your tenacity in getting Excel to do what you want; that looks like many hours of work to count mismatches.

The hammer I use at the moment for this sort of stuff is LAST, which will do *local* alignments of sequences to references. If the number of mismatches are not too large (as would be expected for primer sequences), a local alignment will extend to the end of one of the sequences, and I expect that should work fine for your purposes.

I'm going to demonstrate two options, one that can be done just using Excel and web services, and one that needs command-line access. Note that these are only two ways to fetch mismatch counts; you might want to do more than that (it's not clear from your question), and there are *many* other ways. Also, this may not work for your particular circumstances, but I think it should be quick enough that it's worth a try:

# Option 1 - Web-based LAST

In this, I use the Web version of LAST to generate alignments for the primers \[as "reference"\] mapped to the sequences \[as "query"\].

1. To help create the mismatch (or in this case more correctly, "edit distance") statistic, in the *Advanced Parameters* section, set the `a: gap existence cost` to 1, and the `b: gap extension cost` to 1.
2. Click on `Align!`
3. Once the Alignment Results appear, click on the `Alignments in tabular format` link down the bottom, and save the file to your computer as `last.tab.csv` (the `.csv` at the end is the most important bit, to make sure that Excel recognises it as something it can open).
4. Open that file in Excel, set the delimiter to `<tab>`
5. The first column should be the alignment score, which in this case is the *match* score. This can be converted into a *mismatch* score by subtracting it from the primer sequence length (column 6).

# Option 2 - Local LAST Installation

**Installation**

Install LAST and samtools.

**Index Generation**

The first step is to generate a reference index. LAST doesn't care too much about what's query and what's reference, but because I want all the genome sequences to appear on each line, the genome sequences will be the *query*, and the primer sequences will be the *target*. I use Martin Frith's lastdb mantra to do this:

```
$ lastdb -uNEAR -R01 primers.fa primers.fa

```

**Mapping**

The next step is to do the actual alignment with LAST. I set the match and mismatch to be both 1 (i.e. the default for the Web version), set the gap existence and extension to be 1, change the output format to `tab`, and save the result to `last.tab.csv`:

```
$ lastal -r 1 -q 1 -a 1 -b 1 primers.fa sequences.fa -f tab > last.tab.csv

```

**Aggregating data**

Finally, it's a matter of aggregating the count information. I use `grep -v` to exclude lines that begin with a `#`, use `awk` to subtract field 1 from field 6, then `sort` the results and use `uniq -c` to count the mismatches:

```
$ grep -v '^#' last.tab.csv | awk '{print ($6 - $1)}' | sort | uniq -c
2118 0
   5 1
   5 2

```

---
Tags: sequence-alignment, blast, sequence-analysis, primer
---


# Question
Title: bcl2fastq2 v2.20.0 converstion software on Ubuntu 18.04
bcl2fastq is a program that converts Illumina run output from .bcl to fastq. I'm having difficulties installing the software I downloaded from Illumina on my Ubuntu 18.04. I've followed their installation manual, but the software is not building when I enter the command "make". I receive the error message... "make: \*** No targets specified and no makefile found. Stop. Has anyone had similar experience? Is there another software that will install easier on Ubuntu?

Software: https://support.illumina.com/downloads/bcl2fastq-conversion-software-v2-20.html

Manual: https://support.illumina.com/content/dam/illumina-support/documents/documentation/software\_documentation/bcl2fastq/bcl2fastq2-v2-20-software-guide-15051736-03.pdf

Steps followed:

```
export TMP=/tmp
export SOURCE=${TMP}/bcl2fastq
export BUILD=${TMP}/bcl2fastq2-v2.20.0-build
export INSTALL_DIR=/usr/local/bcl2fastq2-v2.20.0

cd ${TMP}
tar -xvzf bcl2fastq2-v2.20.0.499-Source.tar.gz

mkdir ${BUILD}
cd ${BUILD}
chmod ugo+x ${SOURCE}/src/configure
chmod ugo+x ${SOURCE}/src/cmake/bootstrap/installCmake.sh
${SOURCE}/src/configure --prefix=${INSTALL_DIR}

cd ${BUILD}
make

```

The final command is supposed to be `make install`, but the `configure` step fails with:

```
[ . . . ]
-- Configuring incomplete, errors occurred!
Couldn't configure the project:
/tmp/bcl2fastq2-v2.20.0-build/bootstrap_cmake/bin/cmake -H'/tmp/bcl2fastq/src' -B'/tmp/bcl2fastq2-v2.20.0-build' -G'Unix Makefiles' -DBCL2FASTQ_NAME_SHORT:STRING="bcl2fastq" -DBCL2FASTQ_NAME_LONG:STRING="BCL to FASTQ file converter" -DBCL2FASTQ_COPYRIGHT:STRING="Copyright (c) 2007-2017 Illumina, Inc." -DBCL2FASTQ_VERSION_MAJOR:STRING="2" -DBCL2FASTQ_VERSION_MINOR:STRING="20" -DBCL2FASTQ_VERSION_PATCH:STRING="0" -DBCL2FASTQ_VERSION_BUILD:STRING="422" -DBCL2FASTQ_VERSION:STRING="2.20.0.422" -DBCL2FASTQ_SOURCE_DIR:STRING="/tmp/bcl2fastq/src" -DBCL2FASTQ_PREFIX:PATH="/usr/local/bcl2fastq2-v2.20.0" -DBCL2FASTQ_EXEC_PREFIX:PATH="" -DCMAKE_INSTALL_PREFIX:PATH="/usr/local/bcl2fastq2-v2.20.0" -DBCL2FASTQ_BINDIR:PATH="" -DBCL2FASTQ_LIBDIR:PATH="" -DBCL2FASTQ_LIBEXECDIR:PATH="" -DBCL2FASTQ_INCLUDEDIR:PATH="" -DBCL2FASTQ_DATADIR:PATH="" -DBCL2FASTQ_DOCDIR:PATH="" -DBCL2FASTQ_MANDIR:PATH="" -DCMAKE_BUILD_TYPE:STRING="RelWithDebInfo" -DCMAKE_PARALLEL:STRING="1"
Moving CMakeCache.txt to CMakeCache.txt.removed

```

# Answer
> 4 votes
Install it with conda.

```
conda install -c dranew bcl2fastq

```

# Answer
> 0 votes
Where are you running the software? In general and especially for servers and high performance clusters I would recommend to use a module system like Lmod in conjunction with EasyBuild, which contains configs for bcl2fastq2 that take care of the installation. A module system allows you to use different versions and to keep your environment clean among other benefits.

---
Tags: software-installation
---