# Question
Title: How to find restriction enzyme frequency in whole genome and count the distance between them?
I am having possibly a simple question which is proving difficult.

My objective is to take the reference FASTA file and flat file database with 1500 restriction enzymes with its cutting sites (mostly pattern I would say).

1. I first want to see that how many times each enzyme cut the FASTA file and;
2. second to check the number and size of fragment that each enzyme would give me, this output pregerably would be a CSV file may be.

**Background** I tried to do pattern search using R, but its just collapsing my system. I am looking for bash command may be.

# Answer
Biopython has a module dedicated to working with restriction enzymes. Without any code examples we cannot help you much further.
> 2 votes

# Answer
The `diffHic` package might have solved this issue in R.

```
# Load the genomic data
genome <- biostrings::readDNAStringSet("your_file.fa")

# Run the in silico restriction digest
cuts <- cutGenome(genome, "CCGCGG", overhang=2L) # SacII

```

the `cuts` GRanges object contains the boundaries of each restriction fragment for given genome.
> 2 votes

---
Tags: python, r, linux, fasta
---


title: Using SED command in Nextflow script

text:
# Question
Title: Using SED command in Nextflow script
I was wondering if someone can help me figure out this error message. I’m at the end of my Nextflow pipeline and I want to change the header in a FASTA file and it works when I use this command:

`sed -i "s/^>/>${assembly_id}_${gene_id}_/g" "${assembly_id}/${gene_id}/${assembly_id}_${gene_id}.fasta"`

but it prints out much more than I actually need (i.e., “goodID\_extra\_extra\_extra”). That’s just part of the file name (it includes what assembler was used) but I only want it to add “goodID\_” to the header.

I tried using this command:

`sed -i "s/^>/>${assembly_id%%_*}_${gene_id}_/g" "${assembly_id}/${gene_id}/${assembly_id}_${gene_id}.fasta"`

but I get this error message:

```
N E X T F L O W  ~  version 21.10.6
Launching `geneTargets.nf` [amazing_elion] - revision: b875027a52
Script compilation error
- file : /Users/rimo/Desktop/AWS/nextflow_work/PROJECTS/gene_targets/geneTargets.nf
- cause: Unexpected input: '{' @ line 94, column 32.
   process publish_gene_sequences {
                                  ^

1 error

```

Can someone explain why this is happening? If this is something I can’t do without breaking the script that’s okay the header will still do it’s purpose I just wanted it to be clearer/less noise.

EDIT: this is the process block

```
process publish_gene_sequences {
publishDir params.publish_dir
input:
    tuple val(assembly_id), val(gene_id), path("gene.fasta")
output:
    path("${gene_id}/${assembly_id}_${gene_id}.fasta")
script:
    """
    mkdir -p ${gene_id}
    cp gene.fasta "${gene_id}/${assembly_id}_${gene_id}.fasta"
    sed -i "s/^>/>${assembly_id}_${gene_id}_/g" "${gene_id}/${assembly_id}_${gene_id}.fasta"
    """
 }

```

# Answer
> 5 votes
Inside your script's command string, *assembly\_id* is actually a `java.lang.String` so trying to manipulate it with a Bash operation (i.e. using `${assembly_id%%_*}`) won't produce the expected result. Here's one way to return a substring up to the first underscore; note that we can use single quotes around our sed expression:

```
process publish_gene_sequences {

    publishDir params.publish_dir

    input:
    tuple val(assembly_id), val(gene_id), path("gene.fasta")

    output:
    path("${gene_id}/${assembly_id}_${gene_id}.fasta")

    script:
    def assembly = assembly_id.substring(0, assembly_id.indexOf('_'))

    """
    mkdir "${gene_id}"
    cp gene.fasta "${gene_id}/${assembly_id}_${gene_id}.fasta"
    sed -i 's/^>/>${assembly}_${gene_id}_/g' "${gene_id}/${assembly_id}_${gene_id}.fasta"
    """
}

```

Note also that in the above example, copying the file can be avoided by applying the expression directly and redirecting stdout:

```
script:
def assembly = assembly_id.substring(0, assembly_id.indexOf('_'))

"""
mkdir "${gene_id}"
sed \\
    -e 's/^>/>${assembly}_${gene_id}_/g' \\
    gene.fasta \\
    > "${gene_id}/${assembly_id}_${gene_id}.fasta"
"""

```

---
Tags: fasta, nextflow, errors, command-line, sed
---


# Question
Title: How to identify matching and unique rows across different columns in R data-frame?
I am working on a large data-frame that has 4 columns and each column has variable rows. I am trying to find unique and identical pathways between the 4 columns (each column represents a particular day of treatment with a drug). Here below is a small example.

```
Pathways_1Day <- c("blood","kidney","testis","No","bone","liver","intestine","lungs","ABC","pancreas","Yes","spleen")

Pathways_2Day <- c("blood","kidney","testis","eyes","bone","cells","intestine","cervix","ABC","pancreas","None")

Pathways_3Day <- c("blood","kidney","vessels","lymph","t-cells","liver","intestine","lungs","ABC","epidermis","None")

df<-data.frame(Pathways_1Day,Pathways_2Day,Pathways_3Day)

```

I want to get a summary of the no of pathways that are common between different timepoints (1, 2 and 3 days).

**Important: The no of pathways is not the same for each day.**

I have tried this:

```
All_pathwayNames <- df%>%group_by_all%>%count

```

But the desired output is not what I am trying to get.

There can be different ways to address that. It will be great if I can get matching rows infront of each other across all columns.

Regards

# Answer
When I create the dataframe using the commands you have provided, I get an error, and the data frame can't be created:

```
Error in data.frame(Pathways_1Day, Pathways_2Day, Pathways_3Day) :
  arguments imply differing number of rows: 12, 11

```

What I *think* you were trying to do is something more like this:

```
library(tidyverse)
tibble() %>%
  bind_rows(tibble(day=1, pathway = Pathways_1Day)) %>%
  bind_rows(tibble(day=2, pathway = Pathways_2Day)) %>%
  bind_rows(tibble(day=3, pathway = Pathways_3Day)) -> combined.tbl

```

Producing a long-format data structure that looks like this:

```
> combined.tbl
# A tibble: 34 x 2
     day pathway
   <dbl> <chr>
 1     1 blood
 2     1 kidney
 3     1 testis
 4     1 No
 5     1 bone
 6     1 liver
 7     1 intestine
 8     1 lungs
 9     1 ABC
10     1 pancreas
# ... with 24 more rows

```

When laid out in this way, you can do something like grouping by pathway, and then do a `pivot_wider` to show the sets:

```
> combined.tbl %>%
     group_by(pathway) %>%
     arrange(pathway) %>%
     pivot_wider(names_from=day, values_from=1)
# A tibble: 20 x 4
# Groups:   pathway [20]
   pathway     `1`   `2`   `3`
   <chr>     <dbl> <dbl> <dbl>
 1 ABC           1     2     3
 2 No            1    NA    NA
 3 None         NA     2     3
 4 Yes           1    NA    NA
 5 blood         1     2     3
 6 bone          1     2    NA
 7 cells        NA     2    NA
 8 cervix       NA     2    NA
 9 epidermis    NA    NA     3
10 eyes         NA     2    NA
11 intestine     1     2     3
12 kidney        1     2     3
13 liver         1    NA     3
14 lungs         1    NA     3
15 lymph        NA    NA     3
16 pancreas      1     2    NA
17 spleen        1    NA    NA
18 t-cells      NA    NA     3
19 testis        1     2    NA
20 vessels      NA    NA     3

```

And insert a filter to select only pathways that appear in all three days (i.e. `n() == 3`):

```
> combined.tbl %>%
     group_by(pathway) %>%
     filter(n() == 3) %>%
     arrange(pathway) %>%
     pivot_wider(names_from=day, values_from=1)
# A tibble: 4 x 4
# Groups:   pathway [4]
  pathway     `1`   `2`   `3`
  <chr>     <dbl> <dbl> <dbl>
1 ABC           1     2     3
2 blood         1     2     3
3 intestine     1     2     3
4 kidney        1     2     3

```

Or at least two days (`n() >= 2`):

```
> combined.tbl %>%
     group_by(pathway) %>%
     filter(n() >= 2) %>%
     arrange(pathway) %>%
     pivot_wider(names_from=day, values_from=1)
# A tibble: 10 x 4
# Groups:   pathway [10]
   pathway     `1`   `2`   `3`
   <chr>     <dbl> <dbl> <dbl>
 1 ABC           1     2     3
 2 None         NA     2     3
 3 blood         1     2     3
 4 bone          1     2    NA
 5 intestine     1     2     3
 6 kidney        1     2     3
 7 liver         1    NA     3
 8 lungs         1    NA     3
 9 pancreas      1     2    NA
10 testis        1     2    NA

```
> 1 votes

# Answer
You're misusing data frames. Split the df into 4 vectors and use a Venn diagram generator to denerate a Venn diagram with the 4 vectors.

Adapting the first example from here

```
#untested
venn.diagram(
  x = list(df[,1], df[,2], df[,3], df[,4]),
  category.names = colnames(df),
  filename = 'my_venn_diagram.png',
  output=TRUE
)

```
> 0 votes

---
Tags: r, dataframe
---


# Question
Title: How to count number of motifs in Prosite?
For example I checked motifs in Transferrin P02787 https://prosite.expasy.org/cgi-bin/prosite/ScanView.cgi?scanfile=816050634835.scan.gz

I have hits by profiles: \[2 hits (by 1 profile) on 1 sequence\] and hits by patterns: \[5 hits (by 3 distinct patterns) on 1 sequence\]

Can I say that I detect **4 different motifs**? One motif is searched by profile searching (PS51408 TRANSFERRIN\_LIKE\_4 Transferrin-like domain profile two times) and three motifs are searched by pattern searching (PS00205 TRANSFERRIN\_LIKE\_1 Transferrin-like domain signature 1, PS00206 TRANSFERRIN\_LIKE\_2 Transferrin-like domain signature 2 two times, PS00207 TRANSFERRIN\_LIKE\_3 Transferrin-like domain signature 3 two times).

All the motifs found are equal to 7? Is it right? Can I talk about the different motifs found by the profile and pattern, or are they the same in some way?

# Answer
As per PROSITE database manual:

> Unlike patterns, profiles are usually not confined to small regions with high sequence similarity. Rather they attempt to characterize a protein family or domain over its entire length. This can lead to specific problems not arising with PROSITE patterns. With a profile covering conserved as well as divergent sequence regions, there is a chance to obtain a significant similarity score even with a partially incorrect alignment. This possibility is taken into account by our quality evaluation procedures. **In order to be acceptable, a profile must not only assign high similarity scores to true motif occurrences and low scores to false matches.** In addition, it should correctly align those residues having analogous functions or structural properties according to experimental data.
>
> **Profiles are supposed to be more sensitive and more robust than** **patterns because they provide discriminatory weights not only for the** **residues already found at a given position of a motif but also for** **those not yet found.** The weights for those not yet found are extrapolated from the observed amino acid compositions using empiric knowledge about amino acid substitutability.

Given the less biased and somewhat more balanced result provided by a profile I would say that the sequence contains two motifs that correspond to a PS51408 TRANSFERRIN\_LIKE\_4 Transferrin-like domain profile and provide the bitscore and e-value for the match.

Best,

av
> 2 votes

---
Tags: motifs, prosite
---


title: DESeq2 and EdgeR

text:
# Question
Title: DESeq2 and EdgeR
I am new to using both DESeq2 and EdgeR in Bioconductor used for transforming my RNA expression data.

However, I am struggling to understand their specific purpose, differences between them and underlying theory. The information is required for a thesis with an imminent deadline

# Answer
I think your question has been answered very well here: https://www.biostars.org/p/284775/

If you want to go deeper perhaps you can look at this paper by Li et al. https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02648-4 where they argue that for large scale datasets (multiple conditions) both of these methods are not optimal and instead suggests using Wilcoxon rank-sum test to decrease chances of FDRs.

Best,

av
> 2 votes

---
Tags: deseq2, edger, homework
---


# Question
Title: Running Bionano´s HybridScaffolding from the command line [Error: align0.stdout]
as the title says I´m trying to run Bionanos HybridScaffolding from the command line rather than using their own software interface. This is commercial software to primarily to identify genomewide cancer analysis but also has options for genome assembly.

I installed Perl and all the necessary packages to run the HybridScaffolding Perl script. According to their documentation I started the script with the correct parameters. After 2 minutes it stops giving me the following error:

```
Beginning initial NGS CMAP to BioNano CMAP alignment...
Running command: /home/user/workflow/scripts/bionano/tools/pipeline/1.0/RefAligner/1.0/RefAligner
-ref /home/user/results/bionano/fa2cmap/pac_bio_assembly_CTTAAG_0kb_0labels.cmap
-i /home/user/data/bionano/bCalAnn1_Saphyr_DLE1.cmap.gz
-o align0 -stdout -stderr -maxmem 128 -maxthreads 64
-maxvirtmem 0 -maxmemIncrease 8 -RAmem 3 1 0.9 -M 1 3 -ScaleDelta 0.03 2 -ScaleDeltaBPP
-hashScaleDelta 2 -res 2.6 -resSD 0.7 -FP 0.2 -FN 0.02 -sf 0.10 -sd 0.0 -sr 0.01 -se
0.2 -extend 1 -outlier 0.0001 -endoutlier 0.001 -PVendoutlier -deltaX 12 -deltaY 12
-xmapchim 12 -hashgen 5 11 2.4 1.4 0.05 3.0 1 1 1 -hash -hashdelta 27 -hashoffset 1
-hashrange 0 -hashGrouped 6 11 -hashMultiMatch 50 2 -hashGC 300 -hashT2 1 -hashkeys 1
-HSDrange 1.0 -insertThreads 4 -nosplit 2 -biaswt 0 -T 1e-10 -S -1000 -indel -PVres 2
-rres 0.9 -outlierBC -xmapUnique 12 -AlignRes 2. -outlierExtend 12 24 -resEstimate -f
-mres 0.9 -thetaScale 0.5 -KL 1.5 -KF 1.5 -Kmax 4 -MinFN 0.01 -MinFP 0.16 -MinSF 0.05
-MaxSF 0.15 -MinSD 0.0 -MaxSD 0.0 -MinSR 0.007 -MaxSR 0.015 -MaxSE 0.5 -MinRes 2.0
-MinResSD 0.7 -hashSD 0 1 -hashSR 0 1

ERROR: Cannot open align0.stdout: No such file or directory

```

For context: The HybridScaffolding script contains several scripts which are started automatically over time. One of these is the RefAligner script, which throws the error in this case.

I have a folder called align0 which lies inside the directory that the script is running, therefore I guess this is not the casual "The path is missing" issue..

Inside the Perl Script 'align0.stdout' appears two times:

```
# error check
if (errCheckRefAligner("align0.stdout", "END of output", "Initial alignment complete in $span.", "ERROR: Initial alignment cannot be completed.", 0) != 1)  {
    my $errMsg = "ERROR: Initial alignment cannot be completed.";
    printStatus(\*STATUS_FILE, "align0", $errMsg, 1);
} else  {
    printStatus(\*STATUS_FILE, "align0", "Done", 0);
} # RefAligner check
errCheckRefAligner("align0.stdout", "END of output", "Initial alignment complete in $span.", "ERROR: Initial alignment cannot be completed.", 1);

```

Sadly I don´t have any Perl knowledge at all. Does anyone have experience with Bionanos Software and/or knowledge of Perl and could explain to me what is happening inside this part of the script?

# Answer
I was able to spot the problem. I thought the Server I was working on had at least 256 RAM. Turns out it only had 64. I think HybridScaffolding needs at least 128 RAM. So when I ran the whole HybridScaffolding script on a stronger Server (1000 RAM) it worked just fine without any problems :)

Edit: For the record, I also had to unzip the cmap.gz file because the RefAligner could not read it zipped.
> 2 votes

---
Tags: perl
---


# Question
Title: Removing particular line along with Fasta header in fasta file
I have a fasta file of 3'UTR sequences of the genes. Some of the gene do not have a 3'UTR sequence and reported with `sequence unavailable` but there is FASTA header.

I need to delete or remove all the "sequence unavailable" lines along with their header, i.e. all blank fasta entries.

Please suggest a command with `awk` or `sed` command.

# Answer
> 4 votes
You can use `awk` here, by setting the input (`RS`) and output (`ORS`) record separator, this is what defines a "line" in awk, to `>`. Then, tell it to print records that do not match `Sequence unavailable`. For instance, using this input file:

```
$ cat file
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>baz
Sequence unavailable
>bar1
AGAGATAGAGAGATATA
>asld
Sequence unavailable
>asd
ACTTC

```

You can do:

```
$ awk  'BEGIN{RS=ORS=">"}!/Sequence unavailable/' file
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>bar1
AGAGATAGAGAGATATA
>asd
ACTTC
>

```

That will add a trailing `>`, so you need one more step to remove it:

```
$ awk  'BEGIN{RS=ORS=">"}!/Sequence unavailable/' file | sed 's/>$//'
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>bar1
AGAGATAGAGAGATATA
>asd
ACTTC

```

---

Alternatively, using the `FastaToTbl` and `TblToFasta` scripts I have posted here, you can do:

```
$ FastaToTbl file | grep -v 'Sequence unavailable' | TblToFasta
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>bar1
AGAGATAGAGAGATATA
>asd
ACTTC

```

Both approaches have the advantage of not needing to read the entire file to memory so can be used on files of any size.

# Answer
> 3 votes
Try this, if you have GNU `grep` (the default on Linux systems):

```
grep --no-group-separator -i -B 1 '^[ACGTN]\+$' input.fasta > filtered.fasta

```

Explanation: find all lines that are composed of only A/C/G/T/N (case insensitive search), and show that line plus the preceding line. Don't put any breaks between matches that are not consecutive.

This should work on both multi-line and single line fasta files.

# Answer
> 2 votes
If there are no line breaks within the sequence data (the format I use) then the following will work.

```
perl -0777 -nle 'while (/^(>.+\n^[AGCT-]+)\n/gm){print $1,"\n"}' myfile.fa

```

---

Test

```
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
>baz
Sequence unavailable

```

Output

```
>foo
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA

```

---

**However** If there are line breaks within the sequence data, please use the following ...

```
perl -0777 -nle 'while (/^(>.+\n[AGCT-]{10}.+\n[^>]+)/gm){print $1}' myfile.fa

```

Then pipe to a file,

```
perl -0777 -nle 'while (/^(>.+\n[AGCT-]{10}.+\n[^>]+)/gm){print $1}' myfile.fa > myfile_minus_blanks.fa

```

This works due to a multi-line `m`

Test ...

```
>foo
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
AGAGATAGAGAGATATA
>baz
Sequence unavailable

```

Out ...

```
>foo
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
AGAGTAGATAGAAGATA
>bar
AGAGATAGAGAGATATA
AGAGATAGAGAGATATA

```

**Note** If you've mixed bases, degenerate nucleotides,

`AGCT-` should be replaced with `AGCTRYWSMKNBDHV` (omitting `-` \[not needed\])

# Answer
> 2 votes
One way using :

```
sed -e '/^>/ { N; }; /Sequence unavailable/ { d; }' file.fa

```

For each header line, append the next line of input into the pattern space. Then, if the pattern space contains "Sequence unavailable" delete the pattern space and start a new cycle.

---

Or using :

```
awk '/^>/ { f=$0; next } /^Sequence unavailable/ { next } f { print f; f="" }1' file.fa

```

For each header line, use a variable to store the current record. We can use *next* here to force awk to immediately stop processing the current record and move on to the next record, which is just the next line, since, by default, awk's *record separator* is a newline character. Then, stop processing the current record and move on to the next record if that line begins with "Sequence unavailable". Finally, print the header line if it exists and unset the variable that we used to store it if it does. Once we've printed the header, we can then just print the current record. We can do this by appending *1* to our command string. The *1* on the end just returns true, and by doing so, it performs the default action, which is to print the current record.

---
Tags: phylogenetics, phylogeny, fasta, awk, grep
---


# Question
Title: Non-parametric, background-based tests on proteomics data from Proteome Discoverer v2.5?
I am new to proteomics analysis, so apologies if my question is silly! For context, I am working with proteome samples from the postmortem human brain for a case-control study.

Our lab generally receives an output file of Proteome Discoverer***** (v2.5) from the proteomics center that we work with. From what I understand, PD perform t-tests on the normalized data, but it is a background-based t-test, so there is a weighting factor is applied that takes into account the distribution of ratios for all other proteins and/or the protein abundances themselves? As such, I am able to see an "Abundance Ratio Adj. P-Value".

My problem is that my data **does not meet the assumptions for a parametric test** and it has a small sample size, so the t-test is not valid (as far as I can tell, unless there is some sort of log transformation that PD does to the data before the t-test). Therefore, I would like to run a Mann-Whitney U Test. Can this be done in the PD software?

If not, **how can I do this background adjustment manually such that I can run a Mann-Whitney U Test but still take into protein abundance like in PD?** I do most of my stats in R in case that is relevant.

Thank you so much in advance for your help!

Kelly

---

Package distributed by ThermoFisher

# Answer
Small samples size and parametric tests are not mutually exclusive. Specialized software has been developed to leverage sharing information across measurements to boost power, that is called Empirical Bayes. You can use the limma package for this sort of analysis, contrasting groups. Input is expected to be normalized and on log2-scale. It's an R package hosted at Bioconductor. Don't try a MWU, with small sample size that is massively underpowered.
> 1 votes

---
Tags: proteins, statistics, differential-expression, proteomics, data-preprocessing
---


# Question
Title: Only the first file Pair is processed from Channel.fromFilePairs()
This is my first nextflow pipeline and i am struggling with the .fromFilePairs concept of Nextflow.

I want to process paired-end read files and align them against a reference genome. I use Hisat2 so the index files come in form of 8 files with the suffix .ht2. Hisat needs the file path and the filename without ".ht2". All worked well with a single-end version of the process. The reference also has to be passed through a channel because i use a hisat docker-container and otherwise the files are not mounted.

The idea is to pass read-pairs to a preprocessing step and then pass them (as tuples, like a Channel.fromPairs()) into the hisat process and then further downstream.

My problem is all file-pairs are processed during the first step but only the first pair is passed further downstream. All downstream steps finish with the one pair and the rest simply stays in the output-folder of the preprocessing step.

I tried to simply skip the first step to see if the "connection" between the first and second process is not working but the Channel which works for the first step does not for the hisat step. So i think there is a problem with the combination of multiple read-pairs as input but only one reference file?!

```
process hisat2_paired {
  tag "$sample_id"

  container "${params.container_hisat}"
  publishDir "${params.out_dir}/Nextflow_output/hisat/"

  input:
  tuple val(sample_id), path(input_files)
  path( index_files )

  output:
  path('*.sam') , emit: hisat2_sam_out_ch
  path('*.txt')

  script:
  """
  INDEX=`find -L ${index_files}/ -name "*.1.ht2" | sed 's/.1.ht2//'`
  hisat2 --rna-strandness RF \
          -x \\\$INDEX \
          -1 ${input_files[0]} \
      -2 ${input_files[1]} \
          -S ${sample_id}.sam &> ${sample_id}.hisat_summary.txt
  """

```

```
workflow {

Channel
    .fromFilePairs( "${params.reads_folder}/*_{R1,R2}*f*q.gz")
    .ifEmpty { exit 1, "no fastq files found at given path" }
    .view()
    .set{ read_pair_ch }

Channel
    .fromPath( params.hisat2_index )
    .ifEmpty { exit 1, "no hisat2 index files found - path was empty" }
    //.view()
    .set { hisat2_index_ch }

hisat2_paired( read_pair_ch, hisat2_index_ch )

}

```

I started the workflow with 3 test file-pairs and the Output on the CLI states that only one process is started although if i look at the Channel.view() all 3 file pairs are printed.

```
executor >  local (1)
[2f/30a4d6] process > hisat2_paired (MC3T3-Ko1_S2_L001) [  0%] 0 of 1
[MC3T3-Ko1_S2_L001, [/mnt/c/Path/to/MC3T3-Ko1_S2_L001_R1_001.fastq.gz, /mnt/c/Path/to/MC3T3-Ko1_S2_L001_R2_001.fastq.gz]]
[MC3T3-Ko2_S4_L001, [/mnt/c/Path/to/MC3T3-Ko2_S4_L001_R1_001.fastq.gz, /mnt/c/Path/to/MC3T3-Ko2_S4_L001_R2_001.fastq.gz]]
[MC3T3-Ko3_S6_L001, [/mnt/c/Path/to/MC3T3-Ko3_S6_L001_R1_001.fastq.gz, /mnt/c/Path/to/MC3T3-Ko3_S6_L001_R2_001.fastq.gz]]

```

I was already thinking about memory and resources. Since I prototype on my local PC i set the queue size in the config file to 2. I know that my hisat process needs about 7Gb RAM (for Mus.musculus) per process so there should be enough resources. Also the problem was that the local executor is not even trying to start a second process.

# Answer
> So i think there is a problem with the combination of multiple read-pairs as input but only one reference file?!

Correct. The problem is that the *hisat2\_index\_ch* is a queue channel, but what you want is a value channel so that the (file) values can be read an unlimited number of times without being consumed. Assuming *params.hisat2\_index* is intended to supply a glob pattern, you could simply call the collect operator to produce a list of files as the sole emission:

```
Channel
    .fromPath( params.hisat2_index )
    .collect()
    .set { hisat2_index_ch }

```

You would then need to update your process definition accordingly. Alternatively, have your workflow create the HISAT2 index files for you. This simplifies things considerably, by requiring only a single reference FASTA file as input. Note that value channels are implicitly created as output for a process whose inputs are all value channels. This means in the example below, *hisat2\_index.out* is also a value channel:

```
params.publish_dir = './results'

params.container_hisat = 'quay.io/biocontainers/hisat2:2.2.1--h87f3376_4'

params.reads = '*_R{1,2}.fastq.gz'
params.ref_fasta = 'chr1.fasta'

```

```
process hisat2_index {

    container params.container_hisat
    cpus 8

    input:
    path fasta

    output:
    tuple val(fasta.baseName), path("${fasta.baseName}.*.ht2")

    """
    hisat2-build \\
        -p ${task.cpus} \\
        "${fasta}" \\
        "${fasta.baseName}"
    """
}

```

```
process hisat2_paired {

    tag { sample_id }

    container params.container_hisat
    publishDir "${params.publish_dir}/hisat"

    input:
    tuple val(sample_id), path(input_files)
    tuple val(idx), path('hisat2-idx/*')

    output:
    path('*.sam'), emit: hisat2_sam
    path('*.txt')

    script:
    def (fq1, fq2) = input_files

    """
    hisat2 \\
      --rna-strandness RF \\
      -x "hisat2-idx/${idx}" \\
      -1 ${fq1} \\
      -2 ${fq2} \\
      -S ${sample_id}.sam \\
      &> ${sample_id}.hisat_summary.txt
    """
}

```

```
workflow {

    read_pairs = Channel.fromFilePairs( params.reads )
    ref_fasta = file( params.ref_fasta )

    hisat2_index( ref_fasta )

    hisat2_paired( read_pairs, hisat2_index.out )
}

```
> 1 votes

---
Tags: nextflow, rna-alignment
---


title: Can I discard missing alleles?

text:
# Question
Title: Can I discard missing alleles?
I am converting a biallelic VCF into a SQL table and one of my tables will be something like:

```
SAMPLE | custom variant identifier | Genotype
-------|---------------------------|---------
0001   | 3456                      | ./.

```

My question is that for the diploid `./.` or haploid `.` missing alleles, whether I can just drop them to save space. My thinking is that since the dot represents missing data, the `./.` are just placeholders in the VCF but provide no information at all.

Is this correct? Should I retain the rows with `./.`?

Showing the end of one line of the VCF (really long lines because of high N):

```
FORMAT          SampleX             SampleY
GT:AD:DP:GQ:PL  0/0:3,0:3:4:0,4,45  ./.:3,0:3:.:0,0,0

```

I thought to also show the result after `bcftools norm -m` to split to biallelic and queried per sample.

```
sample_id   chr@pos@ref@alt genotype
sample1 chr1@10158@ACCCT@A       0/0
sample2 chr1@10158@ACCCT@A       ./. <- does this provide any info or can I discard?
sample3 chr1@10158@ACCCT@A       0/0
sample4 chr1@10158@ACCCT@A       ./.
sample5 chr1@10158@ACCCT@A       0/0

```

# Answer
> 1 votes
If you can, do not discard missing allele information.

Generally you often benefit from being able to distinguish between a position that you interrogated and genotyped as REF (reference) versus one that you interrogated but could not determine the genotype (e.g due to low coverage in the sequencing or probing issues in a microarray setup) and which usually is encoded as a missing (./. if ploidy = 2) genotype. A missing genotype tell us that the genotype at this position could be REF, but could also be some other base, we simply do not know.

Knowing this is critical and/or useful depending the downstream analysis. For example it is quite important and informative when merging datasets. Since formats like VCF or PLINK do not record non-variant positions, if you discarded missing genotypes prior to merging, how can you tell between genotypes that truly were REF in one dataset (fully invariant position) and some recorded missing genotypes at the same position but which you decided to discard? You will end up assuming REF status for everything not recorded and this will inevitably lead to batch effects in your combined dataset. For a better discussion see this blog post

On the other hand, if your relational database (your particular situation) for instance will not incorporate external databases in order to produce joint datasets, perhaps discarding missing genotypes is all fine.

---
Tags: vcf, genomics, bcftools
---


# Question
Title: Using information criteria weights to create consensus tree
While I am familiar with calculating information criteria for various ml models of a multiple sequence alignment, I am not aware of how to use the derived weights to create a tree that combines multiple models weighted accordingly.

Thus, I have X trees, computed by multiple different models on the same alignment. Each has an associated BIC. BICs can be converted to weights. I wish to average the trees according to the BIC derived weights.

I would prefer to be able to do this with R.

# Answer
> 0 votes
Calculating a consensus is easy, `RAxML` has a tree consensus module, however calculating it according to weights is a rare procedure. `ETE3` (Python) will have a consensus library and I assume in R it will be `ape`

The only way I can see this being done is to use the relative BIC as a duplicate factor of the number of trees of a given model and then perform a global tree consensus on the resulting treefile. Thus if Bayes Information Criterion top value is twice the size as the lower value the number of 'top value BIC trees' is double that of the lower value. That could be coded up to be cleaner. . Alternatively a multiplication factor could be introduced to ensure the number of duplicated trees within a scale of 100 - to give a cleaner more interpretable 'bootstrap' (but its not a bootstrap) score for each branch.

Here the consensus frequency is a direct function of BIC, essentially you just have an alternative measure to the "non-parametric bootstraps".

---

I better finish this question. The solution to the problem is to produce separate trees for each model, provide a method of robustness, e.g. bootstrapping and then produce a single tree with all the bootstraps for each model along the branches. That is the formal approach if you're not going to simply select the model with highest probability, highest BIC in this case. This however is not the question.

The problem with the solution above the line ... is the question. I was sort of being polite when I said 'rare procedure'. Weighting a consensus tree is not done in this context, weighting a locus or specific site - alignment position - can be done with good justification within the tree building method. The problems sort of multiple when it is weighted according to exact difference in BIC because the probability of a given branch being robust - which does use consensus is calculated differently. What you are doing is assigning shifts in the BIC/maximum likelihood not to the global topology but to specific incongruent branches. Thats not how it works. The BIC/maximum likelihood - whatever - is for the whole tree. The difference between models is the whole tree - every single branch change - not just a branch here and a branch there because its the score (probability ML) of a singular hypothesis. So a branch change here or there is a different hypothesis - for which you might not have generated a BIC score. The branch here and a branch there is MCMCMC consensus or bootstrap resampling consensus not precise shifts in BIC/maximum likelihood. You would need to dedicate the entire project to this specific method and even then it might be not accepted ... its not as a passing calculation.

There could be biological situations where weighting a consensus tree is doable in the way described here. If a specialist reviewer saw the calculation they would refuse to permit publication because it violates the foundations of tree building theory. What you may not be aware is the long historic that has led to the current acceptance of modern methods and what you are proposing is dismantling that theory.

---
Tags: phylogenetics, phylogeny
---


# Question
Title: What are the columns in bedtools coverage output hist file?
I am using bedtools to caculate the coverage of my targets of my WES data, to later plot. But to plot, I need to know what to plot and what it is I am seeing. I have unsuccefully tried to find what exactly the columns are in the output of the bedtool coverage, and how to read them. Also, the documentation does not help me a lot, and I find it confusing. Does anyone know, or where to find this a bit more elaboratly explained?

I used the following code:

`bedtools coverage -hist -abam input.bam -b input.bed > output >2 log`

**Correction note**: the code above is actually wrong (thanks to @Alexlok) and should be: `bedtools coverage -hist -a input.bed -b input.bam > output 2> log`

And the output is (first line):

Last 5 lines:

Thank you in advance

# Answer
Found it, it is an extended .bed file format, so not specific to bedtools output: source bed file output. With 4 added columns to the end from bedtools -hist option: bedtools -hist option
> 1 votes

---
Tags: sequence-analysis, bedtools, coverage, exome
---


title: Filter with bcftools

text:
# Question
Title: Filter with bcftools
I want to filter a SNP, specifically `CHROM1:POS:630128`. I tried to use `bcftools` for this.

I came up with the following variant of a bcftools call:

```
bcftools filter -e "CHROM=1&POS=63018" /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz

```

which leads to:

```
Error: cannot use arithmetic operators to compare strings and numbers

```

and the call with unquoted parameters leads to

```
$ bcftools filter -e CHROM=1&POS=63018 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz

sh: 1: /home/oem/Documents/GitHub/Genomics_prac_guide/trio/file/genome_file.vcf.gz: Permission denied
Failed to read from standard input: unknown file type

```
* Q1: Is bfctools filter the correct tool here to use?
* Q2: Regarding this comment the first error might indicate that I have to annotate the CHROM & POS Tags in the header file. Is this correct?

An example of the file:

```
##fileformat=VCFv4.2
##fileDate=20221104
##source=PLINKv1.90
##contig=<ID=1,length=248928329>
##contig=<ID=2,length=242157857>
##contig=<ID=3,length=198157342>
##contig=<ID=4,length=190106194>
##contig=<ID=5,length=181288140>
##contig=<ID=6,length=170595721>
##contig=<ID=7,length=159331484>
##contig=<ID=8,length=145075470>
##contig=<ID=9,length=138211490>
##contig=<ID=10,length=133659511>
##contig=<ID=11,length=135075227>
##contig=<ID=12,length=133261768>
##contig=<ID=13,length=114341522>
##contig=<ID=14,length=106873816>
##contig=<ID=15,length=101922277>
##contig=<ID=16,length=90172890>
##contig=<ID=17,length=83203771>
##contig=<ID=18,length=80252738>
##contig=<ID=19,length=58586567>
##contig=<ID=20,length=64323573>
##contig=<ID=21,length=46679699>
##contig=<ID=22,length=50776369>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0
1   74792   rs13328684  N   .   .   .   PR  GT  ./.
1   630128  rs9283150   A   .   .   .   PR  GT  0/0
1   791532  i713426 N   .   .   .   PR  GT  ./.
1   792461  rs116587930 G   .   .   .   PR  GT  0/0
1   817341  rs3131972   G   .   .   .   PR  GT  0/0
1   818725  rs12184325  C   .   .   .   PR  GT  0/0
1   820888  rs12567639  A   .   .   .   PR  GT  0/0
1   823656  rs114525117 G   .   .   .   PR  GT  0/0
1   841166  rs12124819  G   A   .   .   PR  GT  0/1
1   858952  rs12127425  G   .   .   .   PR  GT  0/0
1   866156  rs79373928  T   .   .   .   PR  GT  0/0
1   880041  rs72888853  N   .   .   .   PR  GT  ./.

```

# Answer
> 3 votes
First of all, bcftools might be a bit of overkill if you want one specific variant. All you need is GNU `grep` (the default on Linux):

```
zgrep -P '^(#|1\t630128\t)' file.vcf.gz

```

Or, with any `grep`:

```
zgrep -E '^(#|1'$'\t''630128'$'\t'')' file.vcf.gz

```

That will print all lines that start with a `#` (the headers) or that start with `1`, followed by a tab, then `630128` and another tab.

Alternatively, you can use `awk`:

```
zcat file.vcf.gz |
    awk -F'\t' '/^#/ || ( $1==1 && $2==630128)'

```

---

All that said, yes, you can also do it with `bcftools` and in most cases that is the right tool for the job so you might as well use it here too to get used to it. However, you need the right command and that's `view`, not `filter`:

```
bcftools view -r 1:630128 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/files/genome_file.vcf.gz

```

Note that the file needs to be indexed, so if you have not done so already, run

```
tabix file.vcf.gz

```

---

Also, the reason your unquoted command failed as it did is because the `&` has special meaning to the shell, it means "run this command in the background". Therefore, when you ran `bcftools filter -e CHROM=1&POS=63018 file.vcf.gz` that actually ran two separate commands:

```
bcftools filter -e CHROM=1&

```

And then:

```
POS=63018 file.vcf.gz

```

Which means "set the variable `POS` to the value `63018` and execute the command `file.vcf.gz`. Since the file isn't executable, you got:

```
sh: 1: file.vcf.gz: Permission denied

```

---
Tags: vcf, filtering, bcftools
---


# Question
Title: Using very closely related strains to increase coverage for short read, de novo assembly
Do you think it's possible to combine short read Illumina libraries (WGS) from multiple closely related eukaryotic microbial strains (e.g. libraries from a re-sequencing study, \>99% ITS1 sequence) to boost read coverage during a de novo genome assembly run (e.g. with SPAdes)? My target genome is much larger than anticipated, and I'm wondering if increasing read coverage by an order of magnitude could improve length accuracy, contiguity, etc.

Maybe differences in repetitive regions or transposable element movement prevent such an approach? Or maybe kmer frequencies are affected enough that de brujin graph construction gets messed up.

There seem to be several pitfalls, but maybe benefits could outweigh the costs here. If anyone else has done this or considered it, I'd be interested in your thoughts. Thanks!

# Answer
Given that the differences between your strains probably involve rearrangements and copy number variations, this approach seems inadvisable if your goal is a highly contiguous genome. You will make the dBG into more of a hairball if such differences exist.

Illumina is quite cheap these days, I would vastly prefer simply sequencing more deeply if your issue in inadequate coverage. Online tools exist for figuring out how to get there.

Unfortunately, you will swiftly plateau in terms of contiguity even for quite simple genomes using Illumina. So you can expect quickly diminishing returns.

Knowing more stats on your assembly would help evaluate this. For example, if you are working with effectively 10X coverage with N50 ~ 500bp, I might be willing to look the other way while you co-assemble such a synthetic metagenome, because you will at least be able to assemble gene length contigs, which are unlikely to vary too dramatically between your strains. However, if the issue is that your assembly has an N50 ~ 10Kbp and you want it to be higher, then you should consider using other approaches to improve your assemblies, such as long reads. In a pinch, you could do something like optical maps or Hi-C, though I believe that is not preferred these days given much cheaper/better long reads. (Full disclosure: I used to work at a Hi-C company and still hold stock.)

The field of metagenomics is based on the assumption that your suggested approach is ok to get a basic idea of e.g. gene complement, and if you are clever you can even deconvolute such mixtures into genomes, under the assumption that different genomes are at most distantly related. But that doesn't mean that it's a good idea to do it in silico by combining isolates.
> 2 votes

---
Tags: assembly, genomics
---


# Question
Title: What does the "filters" argument in "getBM" do?
In the `biomaRt` package, there is a function `getBM` which among a few things is useful for mapping between different gene representations. Currently, I'm using it as follows, to map Ensembl transcript IDs to gene names and Ensembl gene IDs.

```
ensemble2gene <- getBM(attributes=c("ensembl_transcript_id",
                                    "external_gene_name",
                                    "ensembl_gene_id"),
                       values = as.list(transcripts),
                       mart = mart)

```

How is the following code different:

```
ensemble2gene <- getBM(attributes=c("ensembl_transcript_id",
                                    "external_gene_name",
                                    "ensembl_gene_id"),
                       filters = "ensembl_transcript_id",
                       values = as.list(transcripts),
                       mart = mart)

```

It's unclear to me from the description of the `filters` argument here. But just based on results, I get less results when I run the second one.

# Answer
Hi An Ignorant Wanderer,

Filters allow you to narrow down the dataset to query a specific subset of features within your dataset. A description and tutorial of the filtering steps within BioMart can be found here: https://www.ensembl.org/info/data/biomart/how\_to\_use\_biomart.html

If you remove filters, you will query the full dataset (all human genes, in the example above). However, this is not recommended for BioMart as large, genome-wide queries will overload the BioMart servers. For queries of this size, we suggest using the Ensembl REST API or FTP site.
> 0 votes

---
Tags: ensembl, biomart
---


# Question
Title: Why do I have missing values returned from getBM when converting Ensembl transcript IDs to gene names?
When I run the following R code to convert Ensembl transcript IDs to gene names and Ensembl gene IDs:

```
mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
ensemble2gene <- getBM(attributes=c("ensembl_transcript_id",
                                    "external_gene_name",
                                    "ensembl_gene_id"),
                       filters = "ensembl_transcript_id",
                       values = as.list(transcripts),
                       mart = mart)

```

some of rows returned have NA values for both the gene name and Ensembl gene ID. This is a bit odd because for example, one of the transcript IDs with missing gene IDs is ENST00000414213 but when I google it, I find that it corresponds to gene ENSG00000235143 (according to this). Why am I getting missing results in `getBM` then? Is it because maybe the transcript is from an old version of the database? If so, is that the only reason or are there other reasons this could be happening?

# Answer
Hi An Ignorant Wanderer,

Yes, you're right - ENST00000414213 is a deprecated ID from Ensembl release 75.

https://www.ensembl.org/Homo\_sapiens/Transcript/Idhistory?t=ENST00000414213

To convert a list of IDs from GRCh37, you will need to use the GRCh37 dataset. You can do this using the following documentation/example: https://www.ensembl.org/info/data/biomart/biomart\_r\_package.html
> 0 votes

---
Tags: r, ensembl, biomart
---


# Question
Title: Multivariate analysis of cell culture bioprocess
I am seeking statistical support for the following bacterial growth time series.

In the lab, we have cell cultures growing in bioreactors. Every day we measure the concentration in *m*mol<sup>-1</sup> of a set of metabolites and amino acids in the media with a cell culture analyser.

Some metabolites are consumed by the cells and some are produced during the time-course of the culture. Along with this, we also count the number of cells every day to calculate the cell growth dynamics.

We wish to assess the relationship between the consumption or production of any of these compounds and cell growth. For example, if the consumption of a specific metabolite from the media has a significant effect on increasing the cell count and would assist media formulation.

One measurement per day is taken over four consecutive days from multiple bioreactors.

What kind of analysis is would be appropriate?

# Answer
This is a machine learning application. Its clear cut. Whether you want to perform unsupervised learning is your call.

Input data is for the cell growth, metabolite concentration time series would be:

> T1 ... T4 (cell growth) + T1 ... T4 (metabolite 1) etc ...

The training target is:

> media type

I would start with a simple naive Bayes model and see how it goes. If you obtain a high accuracy score you're in business, its worked. You would simply plot the result. The classification may not be linear in which case an alternative model is used.

You can manually test it, by inputing the data of cell growth etc ... and seeing if it correctly predicts the media type.

You then ask the ML algorithm which the metabolite/cell growth the 'media type' is influencing and it will tell you. Setting it up is complex I agree, but thats the gist and was the question.

Just to finish, if you want to limit what a given media is in particular influencing, i.e. not loads of measurements but just a few inputs/variables, e.g. a metabolite at a particular time point and/or cell growth at a given time point, then ML delivers that information. Thats a key part of its strength. Thus you can then trial loads of different media - even stuff not trained - and just measure a small number of select inputs/variables rather than the entire kitchen sink - which is what the OP is doing right now.

---

The overall input of ML into biological problems is a bit worrying in truth and tends to be over-represented in a select number of fields, protein modelling, and zero in everything else.
> 1 votes

---
Tags: statistics, machine-learning, bacteria
---


# Question
Title: eQTL data extraction from database JSON file
I have a long list of JSON data, with repeats of contents similar to followings.

Due to the original JSON file is too long, I will just shared the hyperlinks here. This is a result generated from a database called RegulomeDB.

Direct link to the JSON file

I would like to extract specific data (eQTLs) from `"method": "eQTLs"` and `"value": "xxxx"`, and put them into 2 columns (tab delimited) exactly like below. Note: `"value":"xxxx"` is extracted right after `"method": "eQTLs"`is detected.

```
eQTLs   firstResult, secondResult, thirdResult, ...

```

In this example, the desired output is:

```
eQTLs   EIF3S8, EIF3CL

```

I've tried using a python script but was unsuccessful.

```
import json
with open('file.json') as f:
    f_json = json.load(f)
    print 'f_json[0]['"method": "eQTLs"'] + "\t" + f_json[0]["value"]

```

Thank you for your kind help.

# Answer
Here a `jq` solution if you can use it:

```
jq -r --arg f "eQTLs" '$f + "\t" + ([."@graph"[] | select(.method == $f) | .value] | join(", "))' file.json

```

Explanations:

* `-r`: Raw string result (so result string will be without between `"`)
* `--arg name value`: Your filter on `eQTLs` method value
* Search in `@graph` array elements with `[]`
* Select only array element with `method` attribute equals to `eQTLs` (in `$f` argument)
* The result, 2 in this case, are put in a array `[...]` and join with string `,`
* `$f` string, a tabulation `"\t"` and the join result are concatenated to one string with `+` operators

Result:

```
eQTLs   EIF3S8, EIF3CL

```
> 3 votes

# Answer
It's a little hard to know what your goal is here based on the sample code, but I do find that your JSON snippet does not match any spec. I was able to edit it to add brackets to make the entries elements of a list:

```
[
  {
  "biosample_ontology": {
    "@id": "https://www.encodeproject.org/biosample-types/primary_cell_CL_0000576/",
    "@type": [
      "BiosampleType",
      "Item"
    ],
    "aliases": [],
    "cell_slims": [
      "hematopoietic cell",
      "monocyte",
      "myeloid cell",
      "leukocyte"
    ],
    "classification": "primary cell",
    "dbxrefs": [],
    "developmental_slims": [
      "mesoderm"
    ],
    "name": "primary_cell_CL_0000576",
    "organ_slims": [
      "blood",
      "bodily fluid"
    ],
    "references": [],
    "schema_version": "1",
    "status": "released",
    "synonyms": [],
    "system_slims": [
      "immune system"
    ],
    "term_id": "CL:0000576",
    "term_name": "monocyte",
    "uuid": "20685264-1187-4c5c-b4ce-dc85ea53fc11"
  },
  "chrom": "chr11",
  "dataset": "https://www.encodeproject.org/annotations/ENCSR464XOM/",
  "dataset_rel": "/annotations/ENCSR464XOM/",
  "documents": [],
  "end": 59936979,
  "file": "ENCFF814NRI",
  "method": "eQTLs",
  "start": 59936978,
  "strand": null,
  "targets": [],
  "value": "MS4A4A"
},

    {
  "biosample_ontology": {
    "@id": "https://www.encodeproject.org/biosample-types/primary_cell_CL_0000576/",
    "@type": [
      "BiosampleType",
      "Item"
    ],
    "aliases": [],
    "cell_slims": [
      "hematopoietic cell",
      "monocyte",
      "myeloid cell",
      "leukocyte"
    ],
    "classification": "primary cell",
    "dbxrefs": [],
    "developmental_slims": [
      "mesoderm"
    ],
    "name": "primary_cell_CL_0000576",
    "organ_slims": [
      "blood",
      "bodily fluid"
    ],
    "references": [],
    "schema_version": "1",
    "status": "released",
    "synonyms": [],
    "system_slims": [
      "immune system"
    ],
    "term_id": "CL:0000576",
    "term_name": "monocyte",
    "uuid": "20685264-1187-4c5c-b4ce-dc85ea53fc11"
  },
  "chrom": "chr11",
  "dataset": "https://www.encodeproject.org/annotations/ENCSR464XOM/",
  "dataset_rel": "/annotations/ENCSR464XOM/",
  "documents": [],
  "end": 59936979,
  "file": "ENCFF814NRI",
  "method": "eQTLs",
  "start": 59936979,
  "strand": null,
  "targets": [],
  "value": "MS4A4B"
}
]

```

When I did this, I was able to slightly modify your code to print out each element's `method` key, value pair:

```
import json
with open('/Users/max.press/Downloads/file.json') as f:
    f_json = json.load(f)
    for i in range(len(f_json)):
        print(f_json[i]["method"] + "\t" + f_json[i]["value"])

```

```
% python ~/Downloads/readfile.py
# prints this:
eQTLs   MS4A4A
eQTLs   MS4A4B

```
> 0 votes

# Answer
A (very) simple Python solution:

The script file `readfile.py`:

```
#! /usr/bin/python3

import sys
import json

filename=sys.argv[1]
filter=sys.argv[2]

results=[]
with open(filename) as json_file:
    json_content = json.load(json_file)
    for json_graph_element in json_content["@graph"]:
        if json_graph_element["method"] == filter:
            results.append(json_graph_element["value"])
print("%s\t%s" % (filter, ", ".join(results)))

```

Make the script executable:

```
chmod +x readfile.py

```

Use the script:

```
./readfile.py /path/to/the/file.json  eQTLs

```

Result:

```
eQTLs   EIF3S8, EIF3CL

```
> 0 votes

---
Tags: python, bash, text-processing, data-retrieval
---


# Question
Title: Analysis of 4665 proteins in String database
**Background** I have a data set of 4665 proteins from *Phytophthora cactorum* and want to analyze them in the String database as other database does not have *Phytophthora cactorum* for GO or KEGG analysis. The String database states I can upload a maximum of 2000 proteins at once.

Thus I subset my data into 3 sets and upload. However, it gives/risks statistically biased analysis.

**Question** Is there any possibility I can analyze all of my 4665 protein sequences at once?

Looking forward for the help.

Best regards,

Bikal

# Answer
> 1 votes
I didn't even know STRING has the kind of analysis you are requesting, I don't see it on their webpage. In any case, the standard tool for GO or KEGG analyses, at least back when I was doing this sort of thing, is DAVID.

I have confirmed that it does have Phytophthora cactorum so all you need to do is upload your file there and use their tools.

# Answer
> 0 votes
Just to confirm the OPs question going to https://string-db.org/cgi/organisms and typing "Phytophthora cactorum", then clicking 'random network' results in..

https://string-db.org/cgi/network?taskId=bbh7BVTwU0tH&sessionId=bA71KKVVotqM

I've encountered String for metabolic pathways, so I don't have direct experience with the OPs query.

---

All databases are all limited by skews in underlying wet-lab work, e.g. there are some large taxonomic groups of bacteria that are complete blanks, i.e. no-one has ever worked on them. This however doesn't concern the OP because they are focussing on a single plant species (I think thats what it is, but could be wrong).

The limit on uploading is understandable - these are not our computers and someone else pays for the OPs analysis.

What the OP suggests is the '3 uploads' comprise a sampling bias and this too is very, very understandable regardless of what analysis is performed.

**Example** If you've loads of `COG` 'V' type genes (disease resistance), for example, in one of the uploads, but not evenly mixed this will certainly skew the results.

**Solution** The solution in this scenario is to randomly sample the genes amongst the 3 uploads. This is easy to do in Python via `random` and `itertools` libraries/packages . It would be easy to write a custom script. (separate question).

Normally in a randomisation procedure there needs to be several, e.g. 10, iterations. So essential the OP would need to upload and perform the analysis 30 time comprising 3 groups of 10 replicates.

**Example 2** Equally if the OP is honing a particularly metabolic pathway and the genes got split between the three uploads then visibility of the pathway would be partial.

**Solution** The OP could cluster genes into one of the three uploads by membership to a given metabolic pathway. This is more difficult. This will depend on what the upstream annotation software/algorithm was because often the pathway membership is described in that annotation file. `string` is a strong database for pathway membership.

**Example 3** The skew is taxonomic membership, I don't this this counts here because it is a single species.

**Conclusion**

1. Ultimately, the preferred solution is to perform the analysis locally this way the OP will not need to break up their file. With KEGG the databases can be downloaded.
2. The OP could contact `string` essentially asking if they could supply an account that lifts the upload limit. `string` would likely want a fee if they will do this at all. I honestly have no idea.
3. There are other servers as @terdon recommends, I think BRITE is another one which might have the plant (which is what I think it is) under investigation.
4. The OP could go through solutions 1 and 2 above. Solution 1 is likely to be 'the solution', but there needs to be a supporting script written (separate question).

**APIs** I don't know about `string` but if there was an API this would enable a multiple replicate analysis to be coded without all the manual point and click that it would otherwise involve.

---

My opinion of `Go` is that it is sparsely populated, just saying.

# Answer
> 0 votes
Hej Bikal,

I am not sure why you would like to inspect all 4,000+ proteins at once. I would expect that you have a look on a subset of "interesting" genes, eg from a differential expression experiment.

Yet, since you've asked, yet it is possible to get the data for all proteins at once.

You can download the network for entire genomes and analyze with what ever tool you seem fit.

Alternatively, may I recommend you to have a look at the Cytoscape plugin stringApp. It should help you with the download, analysis (including pathway enrichment), and visualizations. Unlike the String network in your web browser, Cytoscape can easily handle large networks.

The tool's authors have a nice series of video tutorials of how to use it. Hopefully that helps in your work.

---
Tags: sequence-analysis, kegg, stringdb
---


# Question
Title: Do you know any virtual screening libraries with small, soluble compounds?
I had a quick question on virtual libraries. I'm hoping to perform a virtual screen using a library that contains very small, soluble compounds - think glycerol and smaller (\<100MW). As far as I can tell, the majority of virtual libraries contain mostly larger, drug-like compounds. Is anyone aware of any libraries that might contain ideally 10,000+ or so of these tiny soluble compounds?

If not, any idea what the process might be (and how difficult it would be) to construct a custom virtual library enriched for these compounds?

Thank you!!

# Answer
## Lipinski and logP

Most datasets will be filtered by the Lipinski rule, plus in the some cases, such as Enamine Real whether they can be bought. One criterion in the Lipinski rule filter is the logP, or octanol-water partition coefficient. This is super important, because a compound that cannot cross the membrane is likely of limited uses —as a consequence prodrugs are a commonly taken approach.

However, albeit obvious, not all datasets in Zinc are filtered by the Lipinski rule, they just need more browsing...

## Natural?

However, you *may* have a valid reason to eschew logP. Nature tries its best to have intermediate metabolites that have hydrophilic handles attaches, eg. CoA, phosphates and sugars. Therefore, if you want a natural metabolite this is a completely different question. Many of these are not very purchasable, but a PDB-derived dataset of natural ligands will likely be purchasable or easily extractable. Alternatively Zinc has several natural catalogues...

## Too small

Glycerol is a 3 carbon molecule and has a MW of 90 Da. Many of the old drugs are small, but not that small. For example aspirin is 180 Da, which is twice the size of glycerol. Lipinski rule says not to go over 500 Da and 250-400 Da is a lovely range. Say you have a perfectly bound glycerol, with a perfect ligand efficiency of -1.5 kcal/mol/atom, that would mean it would bind at –9 kcal/mol in the best case scenario, i.e. it is too small to bind anything with nanomolar affinity and will actually bind a lot of things. Glycerol after all is commonly cocrystallised in protein...

The GDB datasets have enumerated compounds of a given size, but don't go that low.

## Docking

There will be two further complications that will hinder the accuracy.

* protonated and deprotonated atoms
* polarisable groups

The former requires correct conformers, the latter requires a force field that has a Drude particle, e.g. in openMM.

## Caveat against sorting by binding score

Also, it should be said that if you sort by binding score (kcal/mol or kJ/mol) you will get bigger compounds (hence the whole point of ligand efficiency) and compounds with high logP. The former is because they will have more contact points to make, the latter is because given that ∆∆G is ∆G\_bound - ∆G\_unbound, an unbound high-logP compound is highly unhappy in *aqueous* solution.
> 1 votes

# Answer
ZINC20 is a free database of commercially-available compounds for virtual screening. It contains 1.3M compounds with molecular weight less than 200. https://zinc20.docking.org/

Cartblanche22 is an interface to ZINC22. Here you can find more compounds than ZINC20. https://cartblanche22.docking.org/cartblanche
> 1 votes

---
Tags: docking, computational-biochemistry, cheminformatics
---


# Question
Title: finding matching motifs in a sequence and their position
I am trying to find some matching motifs on a sequence, as well as the position that the motif is located in and then output that into a fasta file. The code below shows that the motif `[L*L*L]` is present in the sequence, when I run it returns as "YES" but I do not know where it is positioned.

```
peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

if re.search(r"L*L*L", peptide1):
    print("YES")
else:
    print("NO")

```

The code that I wrote to find the position is below, but when I run it says invalid syntax. Could you please assist as I am still new in Bioinformatics, and I have no clue whether in the right track or not.

```
for position in range(len(peptide1)):
    if peptide1[position:].startswith(re.search(r"L*L*L")):
        print(position+1)

```

Lastly once the position of these motifs has been identified, how do I highlight them such that they are visible in the output fasta file.

I am searching for `LxLxL` where `x` can be any amino acid. I used the `*` to indicate to python that it should use any amino acid, but I think I am incorrect. Could you please assist on how I should actually be writing that motif?

# Answer
> 2 votes
`re.search` method do not return a simply boolean but a useful object contains many information:

```
#! /usr/bin/env python3

import re

peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

re_search_result = re.search(r'L*L*L', peptide1)
if re_search_result:
    print("YES")
    motif_start_position=re_search_result.start()
    motif_end_position=re_search_result.end()
    motif_length = re_search_result.end() - re_search_result.start()
    print("Motif start at position: %d" % motif_start_position)
    print("Motif end at position  : %d" % motif_end_position)
    print("SO length is ......... : %d" % motif_length)
    print("The motif is : <%s>" % peptide1[motif_start_position:motif_end_position])
else:
    print("NO")

```

Result is:

```
YES
Motif start at position: 20
Motif end at position  : 488
SO length is ......... : 468
The motif is : <LTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNL>

```

The motif `L*L*L` is a quite strange: `L` zero or many times followed by `L` zero or many times followed by followed by a `L`. It's an equivalent of `L*L`, equivalent of `L+`.

Do you search for `L.*L.*L`, `L[^L]*L[^L]*L` or simply `LLL` ?

Notes:

* When you post code, use block delimited by two \`\`\` (no by two simple \`)
* In Python, code indentation is very important, always use spaces (4)

---

If you want to loop over all motifs (like strings), use `findall` method:

```
#! /usr/bin/env python3

import re

peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

re_findall_result = re.findall(r'L[^L]*L[^L]*L', peptide1)
for element in re_findall_result:
    print("Element: <%s>" % element)

```

Result will be:

```
Element: <LTGALL>
Element: <LYRDAAVESPQHEVSNESGSTLKEHDYFGL>
Element: <LNLNESATTL>
Element: <LGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNIL>
Element: <LDLKVHL>
Element: <LSFDNMQGPLNPFFRKSL>
Element: <LVKIYMDGVPFGRKVDLKTNDSYDKL>
Element: <LEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFL>
Element: <LIYEDHEGDSMLVGDVPWEL>

```

---

If you want all find objects like with `search` method, use `finditer` method:

```
#! /usr/bin/env python3

import re

peptide1="MKFSNEVVHKSMNITEDCSALTGALLKYSTDKSNMNFETLYRDAAVESPQHEVSNESGSTLKEHDYFGLSEVSSSNSSSGKQPEKCCREELNLNESATTLQLGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNILEEFRAMKAQAHMTKSPKPVHTMQHNMHASFSGAQMAFGGAKNNGVKRVFSEAVGGNHIAASGVGVGVREGNDDVSRCEEMNGTEQLDLKVHLPKGMGMARMAPVSGGQNGSAWRNLSFDNMQGPLNPFFRKSLVSKMPVPDGGDSSANASNDCANRKGMVASPSVQPPPAQNQTVGWPPVKNFNKMNTPAPPASTPARACPSVQRKGASTSSSGNLVKIYMDGVPFGRKVDLKTNDSYDKLYSMLEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFLEGSEYVLIYEDHEGDSMLVGDVPWELFVNAVKRLRIMKGSEQVNLAPKNADPTKVQVAVG"

re_findall_result = re.finditer(r'L[^L]*L[^L]*L', peptide1)
for element in re_findall_result:
    print("peptide1[%s:%s]=<%s>" % (element.start(), element.end(), peptide1[element.start():element.end()]))

```

Result will be:

```
peptide1[20:26]=<LTGALL>
peptide1[39:69]=<LYRDAAVESPQHEVSNESGSTLKEHDYFGL>
peptide1[90:100]=<LNLNESATTL>
peptide1[101:163]=<LGPPAAVKPSGHADGADAHDEGAGPENPAKRPAHHMQQESLADGRKAAAEMGSFKIQRKNIL>
peptide1[248:255]=<LDLKVHL>
peptide1[278:296]=<LSFDNMQGPLNPFFRKSL>
peptide1[378:404]=<LVKIYMDGVPFGRKVDLKTNDSYDKL>
peptide1[407:443]=<LEDMFQQYISGQYCGGRSSSSGESHWVASSRKLNFL>
peptide1[449:469]=<LIYEDHEGDSMLVGDVPWEL>

```

---
Tags: fasta, python, motifs, sequence
---


# Question
Title: Is there any book or reference that explains evolutionary dynamics and fixation probability?
I am a master's student in computational and mathematical biology, with a bachelor's degree in pure mathematics. I am really struggling with my computational biology course, we are supposed to study problems related to evolutionary dynamics like: 1)the rock-paper-scissors game which explains the evolution of male mating strategies in side-blotched lizards 2)fixation probability of a new mutation. Could someone link me to some books or references that might be helpful in solving such problems?

# Answer
* Maynard Smith, J. (1989) Evolutionary Genetics. Oxford: Oxford University Press. ISBN 0-19-850231-1
* Maynard Smith, J. and Harper, D. (2003) Animal Signals. Oxford University Press. ISBN 0-19-852685-7\[16\]

Then perform a PubMed search for Adam Eyre-Walker.
> 2 votes

# Answer
I'm fond of Joe Felsenstein's book on evolutionary genetics as a relatively recent source that presents a lot of these problems at a very simple level for the newer student, and then adds complexity gradually. Page 260 seems relevant.
> 2 votes

---
Tags: phylogenetics, phylogeny, genetics, population-genetics
---


# Question
Title: Large amount of REF prefix differences after lifting
I lifted some unmerged trio files from hg19 to hg38 using liftOverPlink. It looks like that several houndred allels on each chromosome of one file shifted/switched to non-standard alleles. I recorded some of the differences. And will add them to the post. I will also add a sample of the file to this post. My goal is to merge the lifted trio files to one singe file.

My command to merge is

```
bcftools  merge -l merge.txt --force-samples -Ob -o file.bcf

The REF prefixes differ: C vs N (1,1)
Failed to merge alleles at 1:3940550 in file_1.vcf.gz

```

First, I tried to remove those alleles in question but I realized that these are just to many to do this by hand.

I used R to edit the file

```
file <- read_delim(
  gzfile(
    str_glue("{some_path}file.vcf.gz")),
  delim="|",  quote='', escape_backslash = T, col_names = F) %>%
#A list of all the REF prefix differences: number equals position in chromsome 1
  #630128 N vs A, 791532, N vs A, 841166 G vs A, 937572 C vs T, 1110226 N vs A, 1486611 N vs G, 1554290 C vs G,
  #1654663 N vs T, 1671724 C vs T, 1755523 C vs T, 1774697 C vs T, 1916092 C vs T, 1949922 C vs T, 1959595 C vs n
  #1959717 C vs T, 1983048 A vs G, 1991690 A vs G, 1999315 C vs T 2096621 C vs T 2107385 A vs G, 2124676 A vs G
  #2274142 A vs C, 2284952 A vs G, 2293812 C vs T, 2303949 N vs C, 2321888 T vs N, 2351827 C vs T,2364905 N vs G
  #2364905 T vs C, 2420707 T vs C, 2431388 T vs C, 2441501 T vs C, 2461209 C vs G, 2614790 A vs G, 2620461 G vs A
  #2622185 C vs T, 2792599 A vs C, 2861494 C vs T 2864702 G vs N, 2884008 T vs C , 2894837 A vs C, 2925918 T vs C,
  #2944366 T vs C, 3097126 G vs A,, 3107701 C vs A, 3109422  C vs A, 3112423 G vs A, , 3119756 C vs T, 3146311 G vs T
  #3282557 G vs A,, 3294386 A vs G, 3294894 A vs C, 3303337 A vs G, 3320610 T vs C, 3335540 T vs C, 3344952 T vs C
  #3345118 T vs C 3355977 G vs A, 3401871 G vs A, 3443779 G vs T,, 3534441 T vs N, 3591297 A vs G,
  #3595370  A vs G, 3601923 T vs C, 3624792 G vs N, 3689917 T vs C, 3689917 G vs A, 3695328 G vs A, 3696503 G vs A
  #3762340 T vs C, 3775163 A vs C, 3870048 T vs Cl 3924314 G vs A3935847 T vs C, 3940550 G vs A
   filter(!grepl('.+630128.+|.+791532.+|.+841166.+|.+937572.+|.+1110226.+|.+1486611.+|.+1554290.+|.+1654663.+|.+1671724.+|.+1755523.+|.+1774697.+|.+1916092.+|.+1949922.+|.+1959595.+|.+1959717.+|.+1983048.+|.+1991690.+|.+1999315.+|.+2096621.+|.+2107385.+|.+2124676.+|.+2274142.+|.+2284952.+|.+2293812.+|.+2303949.+|.+2321888.+|.+2324778.+|.+2351678.+|.+2351827.+|.+2364905.+|.+2364905.+|.+2420707.+|.+2431388.+|.+2441501.+|.+2461209.+|.+2614790.+|.+2620461.+|.+2622185.+|.+2792599.+|.+2861494.+|.+2864702.+|.+2884008.+|.+2894837.+|.+2925918.+|.+2944366.+|.+3097126.+|.+3107701.+|.+3109422.+|.+3112423.+|.+3119756.+|.+3146311.+|.+3282557.+|.+3294386.+|.+3294894.+|.+3303337.+|.+3320610.+|.+3335540.+|.+3344952.+|.+3345118.+|.+3355977.+|.+3401871.+|.+3443779.+|.+3534441.+|.+3591297.+|.+3595370.+|.+3601923.+|.+3624792.+|.+3689917.+|.+3689917.+|.+3695328.+|.+3696503.+|.+3762340.+|.+3775163.+|.+3870048.+|.+3924314.+|.+3935847.+|.+3940550.+',X1)) %>%
  write_delim(file = str_glue("{some_path}some_file_v1.vcf"),
                                 k                        col_names = F, quote = 'none', escape = 'none')

#zip the file
call = str_glue("bcftools view {some_path}some_file_v1.vcf -o {some_path}some_file_v1.vcf.gz")
system(call)
#index the new file
call = str_glue("bcftools index -f {f}.gz")
system(call)

```

**Q1: How can one fix the large amount of REF prefix differences in the chromosome files?**

**Q2: If one uses the fixref plugin of bcftools which is the correct reference file to use after lifting was performed?**

**EDIT:** JRodrigoF recommended to use the `bcftools norm --check-ref s --fasta-ref GRCh38_full_analysis_set_plus_decoy_hla.fa`  to fix the issue. It seems this produces the following error on my end:  `[E::faidx_adjust_position] The sequence "1" was not found faidx_fetch_seq failed at 1:69869`

To make sure this is not a problem created during liftover I will post my liftover procedure here:

```
library(tidyverse)

#set working directory to data directory
trio_wd <- str_glue(here::here(),'/trio/SomeName/')

#create file list for raw data
file_list <- str_c(trio_wd,dir(trio_wd)) %>% str_extract('.+txt') %>%
  {.[!is.na(.)]}

    #generate map files from 23andme text file
    for(f in file_list){
      call = str_glue('plink --23file {f} --recode --out {str_remove(f,".txt")}
                      ')
      system(call)
    }

    #create file list for map data
    file_list <- str_c(trio_wd,dir(trio_wd)) %>% str_extract('.+map') %>%
      {.[!is.na(.)]}

    for (f in file_list){
    call <- str_glue('python3 reference/liftOverPlink.py --map {f} --out {str_replace(f,".map","_lifted")} --chain reference/hg19ToHg38.over.chain.gz  ')
    system(call)

    #create valid mapping file for hg19
    call <- str_glue('python3 reference/rmBadLifts.py -m {str_replace(f,".map","_lifted")}.map --out {str_replace(f,".map","_good_lifted")}.map -l {str_replace(f,".map","_bad_lifted.dat")}')
    system(call)

    call <- str_glue('cut -f 2 {str_replace(f,".map","_bad_lifted.dat")} > {str_replace(f,".map","_exclude_lifted.dat")}')
    system(call)
    call <- str_glue('cut -f 4 {str_replace(f,".map","_lifted")}.bed.unlifted | sed "/^#/d" >> {str_replace(f,".map","_exclude_lifted.dat")} ')
    system(call)
    # Note: this will clobber the lifted MAP file generated by `liftOverPlink`:
    call <- str_glue('plink --file {str_remove(f,".map")} --recode --out {str_replace(f,".map","_lifted")} --exclude {str_replace(f,".map","_exclude_lifted.dat")} ')
    system(call)
    call <- str_glue('plink --ped {str_replace(f,".map","_lifted.ped")} --map {str_replace(f,".map","_good_lifted")}.map --recode --out {str_replace(f,".map","_final")}')
    system(call)
    }

```

**EDIT II:** Sample of the file including header

```
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20221108
##source=PLINKv1.90
##contig=<ID=1,length=248928329>
##contig=<ID=2,length=242157857>
##contig=<ID=3,length=198157342>
##contig=<ID=4,length=190106194>
##contig=<ID=5,length=181288140>
##contig=<ID=6,length=170595721>
##contig=<ID=7,length=159331484>
##contig=<ID=8,length=145075470>
##contig=<ID=9,length=138211490>
##contig=<ID=10,length=133659511>
##contig=<ID=11,length=135075227>
##contig=<ID=12,length=133261768>
##contig=<ID=13,length=114341522>
##contig=<ID=14,length=106873816>
##contig=<ID=15,length=101922277>
##contig=<ID=16,length=90172890>
##contig=<ID=17,length=83203771>
##contig=<ID=18,length=80252738>
##contig=<ID=19,length=58586567>
##contig=<ID=20,length=64323573>
##contig=<ID=21,length=46679699>
##contig=<ID=22,length=50776369>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_viewVersion=1.13+htslib-1.13+ds
##bcftools_viewCommand=view -r 1 -o /home/oem/Documents/GitHub/Genomics_prac_guide/trio/some_name/some_name//genome_some_name_v5_Full_20191220113716_final_1.vcf.gz /home/oem/Documents/GitHub/Genomics_prac_guide/trio/some_name/some_name_v5_Full_20191220113716_final.vcf.gz; Date=Tue Nov  8 13:02:01 2022
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0
1   74792   rs13328684  N   .   .   .   PR  GT  ./.
1   792461  rs116587930 G   .   .   .   PR  GT  0/0
1   817341  rs3131972   G   .   .   .   PR  GT  0/0
1   818725  rs12184325  C   .   .   .   PR  GT  0/0
1   820888  rs12567639  A   .   .   .   PR  GT  0/0
1   823656  rs114525117 G   .   .   .   PR  GT  0/0
1   858952  rs12127425  G   .   .   .   PR  GT  0/0
1   866156  rs79373928  T   .   .   .   PR  GT  0/0
1   880041  rs72888853  N   .   .   .   PR  GT  ./.
1   889018  rs7538305   A   .   .   .   PR  GT  0/0
1   894801  rs28444699  A   .   .   .   PR  GT  0/0
1   895351  i713449 N   .   .   .   PR  GT  ./.
1   899450  rs116452738 G   .   .   .   PR  GT  0/0
1   899712  rs72631887  T   .   .   .   PR  GT  0/0
1   903285  rs28678693  T   C   .   .   PR  GT  0/1
1   905373  rs4970382   T   C   .   .   PR  GT  0/1
1   911428  rs4475691   C   .   .   .   PR  GT  0/0
1   916010  rs72631889  G   .   .   .   PR  GT  0/0

```

**EDIT III:** Output of `grep '^>' reference/GRCh38_full_analysis_set_plus_decoy_hla.fa`

```
>chr1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38
>chr2  AC:CM000664.2  gi:568336022  LN:242193529  rl:Chromosome  M5:f98db672eb0993dcfdabafe2a882905c  AS:GRCh38
>chr3  AC:CM000665.2  gi:568336021  LN:198295559  rl:Chromosome  M5:76635a41ea913a405ded820447d067b0  AS:GRCh38
>chr4  AC:CM000666.2  gi:568336020  LN:190214555  rl:Chromosome  M5:3210fecf1eb92d5489da4346b3fddc6e  AS:GRCh38
>chr5  AC:CM000667.2  gi:568336019  LN:181538259  rl:Chromosome  M5:a811b3dc9fe66af729dc0dddf7fa4f13  AS:GRCh38  hm:47309185-49591369

```

**EDIT IV:** Edited the fasta-file and the corresponding index file as mentioned in the comments. Still, the error from EDIT I sustains.

First five chromosomes of the modified FASTA-file:

```
grep '^>' reference/GRCh38_full_analysis_set_plus_decoy_hla_modified.fa
>1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38
>2  AC:CM000664.2  gi:568336022  LN:242193529  rl:Chromosome  M5:f98db672eb0993dcfdabafe2a882905c  AS:GRCh38
>3  AC:CM000665.2  gi:568336021  LN:198295559  rl:Chromosome  M5:76635a41ea913a405ded820447d067b0  AS:GRCh38
>4  AC:CM000666.2  gi:568336020  LN:190214555  rl:Chromosome  M5:3210fecf1eb92d5489da4346b3fddc6e  AS:GRCh38
>5  AC:CM000667.2  gi:568336019  LN:181538259  rl:Chromosome  M5:a811b3dc9fe66af729dc0dddf7fa4f13  AS:GRCh38  hm:47309185-49591369

```

Addjusted the `bcftools`-call of EDIT I afterwards so the modified FASTA-file will be used

```
call = str_glue('bcftools  norm --check-ref s {f} --fasta-ref reference/GRCh38_full_analysis_set_plus_decoy_hla_modified.fa -o {out} ')
[E::faidx_adjust_position] The sequence "1" was not found
faidx_fetch_seq failed at 1:69869

```

View region in question

```
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20221108
##source=PLINKv1.90
##contig=<ID=1,length=248928329>
##contig=<ID=2,length=242157857>
##contig=<ID=3,length=198157342>
##contig=<ID=4,length=190106194>
##contig=<ID=5,length=181288140>
##contig=<ID=6,length=170595721>
##contig=<ID=7,length=159331484>
##contig=<ID=8,length=145075470>
##contig=<ID=9,length=138211490>
##contig=<ID=10,length=133659511>
##contig=<ID=11,length=135075227>
##contig=<ID=12,length=133261768>
##contig=<ID=13,length=114341522>
##contig=<ID=14,length=106873816>
##contig=<ID=15,length=101922277>
##contig=<ID=16,length=90172890>
##contig=<ID=17,length=83203771>
##contig=<ID=18,length=80252738>
##contig=<ID=19,length=58586567>
##contig=<ID=20,length=64323573>
##contig=<ID=21,length=46679699>
##contig=<ID=22,length=50776369>
##INFO=<ID=PR,Number=0,Type=Flag,Description="Provisional reference allele, may not be based on real reference genome">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##bcftools_viewVersion=1.13+htslib-1.13+ds
##bcftools_viewCommand=view -r 1 -o /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name//genome_Some_name_v5_Full_20191220113716_final_1.vcf.gz /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/genome_Some_name_v5_Full_20191220113716_final.vcf.gz; Date=Tue Nov  8 13:02:01 2022
##bcftools_viewCommand=view -o /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name/genome_Some_name_v5_Full_20191220113716_final_1_1.vcf.gz /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name/genome_Some_name_v5_Full_20191220113716_final_1_1.vcf; Date=Tue Nov  8 18:16:55 2022
##bcftools_viewCommand=view -r 1:69869 /home/oem/Documents/GitHub/Genomics_prac_guide/trio/Some_name/Some_name/genome_Some_name_v5_Full_20191220113716_final_1_1.vcf.gz; Date=Thu Nov 10 10:56:54 2022
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  FAM001_ID001
1   69869   rs548049170 T   .   .   .   PR  GT  0/0

```

# Answer
It is important that after lifting over a VCF/BCF file from one coordinate system to another, the header of the resulting file contains the chromosome/contig names as present in the coordinate system of the new reference genome. it works as a double check that the assumptions of the coordinate system (particular reference) are consistent with the actual data (VCF lines). PLINK files for example do not have a header section and this doesn't not apply.

Looks like BCFTools (merge) complains whenever if finds positions that have one base ('C') as reference in one of the input files, and another base ('N') in another file. When working with VCF/BCF files, **bcftools norm** can be used to fix 'N' reference position into their actual base (if that is the case) and/or fix bases listed as reference (e.g. coming from a PLINK major-allele based BED file), but which actually correspond to an alternate allele based on either GRCh37/GRCh38, ... .

```
bcftools norm --check-ref s --fasta-ref <reference_file> <my_VCF_file>

```
> 1 votes

# Answer
Thanks for @JRodrigoF for pointing towards `bcftools norm --check-ref s --fasta-ref`. This will indeed fix the initial REF prefix diff issue of my question.

If anyone else stumbles over this question while facing a similar error to :`[E::faidx_adjust_position] The sequence "1" was not found faidx_fetch_seq failed at 1:69869`

Please try to re-index your FASTA-file using samtools `samtools faidx reference_file.fa` Most of the times this seems to solve this kind of error!
> 0 votes

---
Tags: bcftools, liftover, merge
---


# Question
Title: Is current AI able to generate RNA sequences of viruses?
After seeing a video of a combination of CRISPR and AI, and an article of someone who made two babies immune against HIV and still healthy, I wondered about something. If such complicated thing, to edit someone's genome, so that the person isn't likely to get other diseases or complications, while having a selected function, is (almost) possible, will it be possible currently to let AI generate RNA of virus particles like from a Generative Adversarial Network?

> We show that these tools capture important structures of the data and, when applied to designing probes for protein binding microarrays (PBMs), allow us to generate new sequences whose properties are estimated to be superior to those found in the training data.

https://arxiv.org/abs/1712.06148 In this research article, it is stated that the researchers were successful in training the model so that it would generate new sequences, while having a specific functioning protein. Note that this was from the time when ProGAN wasn't even invented yet, so they used a Wasserstein-GAN.

> While we found these techniques to generally be well-behaved, we did notice some undesirable cases during our explorations.

They concluded that even one mistake in nucleotides could mess the protein up. But now in 2022, we have StyleGAN-Ada (https://github.com/NVlabs/stylegan2-ada), which is far better in generating stuff, we've already (almost) overcome that stuff? So if you basically train such model on RNA of virus particles (most viruses have really short genome compared to an eukaryote cell, so a supercomputer is not required for training such thing, and there are databases online: Publicly available genome sequence database for viruses?), such as every mutation of SARS-CoV-2, HIV, ebola and so on, it would basically be possible to let it generate something that is as contageous as SARS-CoV-2, whereas as deadly as HIV? Keep in mind that one single generated sequence that functions, would already be devastating. It doesn't need to be 100% accurate, because after all we can just experiment with it and find out whether it is working or not as supposed to be.

It doesn't have to be a GAN:

> We present three approaches: creating synthetic DNA sequences using a generative adversarial network (GAN); a DNAbased variant of the activation maximization (“deep dream”) design method; and a joint procedure which combines these two approaches together.

The RNA-sequence can be brought into a cell, and the cell will automatically create the components of the virus, which will become the virus ultimately.

Anyone can recall that this is possible? I have done my searches on Google, but I haven't found any article particularly describing this concept.

I hope this is on-topic, I don't know how the forums work. Thanks in advance.

# Answer
> 1 votes
Reinforced learning is an interesting technology, but they are prone to extreme results - which are not necessarily at all biological feasible. Google can use it to outperform, or very possibly outcompete, the world's best chess player, because all of the rules are known. In biology we simply do not know all the rules.

To reiterate, the problem with reinforced learning in biology is not only do we not know all the rules, in fact in reality we know very few of the rules. HIV - why does it cause persistent infection? SARS-CoV-2 why doesn't the vaccine confer complete protection? Moreover, the GAN originator needs to know all of the biological *known* rules and apply them within the GAN and there are a lot of rules. Does the originator know the conformational change of env or spike and on which amino acid that are invariant? I doubt it. That information would be pivotal to generate AI sequence for a biologically feasible virus.

Even in the original Atari video game experiment ***** whilst GANS outcompeted humans by \>1000%, equally they under-performed humans considerable for around \>20% games tested. Thus even if all the rules are known GANs can mess-up a large proportion of the time.

These are simply hypotheses in an emerging technology within biology that require biological verifiable results and those proofs will be difficult, i.e. time consuming. This is not the case within video games, or board games.

In the first instance, GAN output can be investigated for biological validity against strongly established bioinformatics analysis - the stuff we use to assess whether there are errors in the data. However, there needs to be a lot incentive to perform this research at this point in the technological cycle. When the authors have something concrete, we might then reinvestigate the fundamental issues of their predictions.

*****, the old Atari video games were subject to GANs and the computer then compared its performance in the video game against an average human. The results were highly variable, moreover humans are highly variable - which to my knowledge was ever investigated - so it was never at all clear the "the average human" represented in these tests.

---

A GAN isn't 'trained against \[biological entity\]' in the technical way supervised learning is trained. A GAN 'trains' against itself via a 'win/lose' quantification. Establishing a valid 'win/lose' is a pivotal part of the GAN. Thus in chess thats easy ... you win or lose the game, in biological that is not trivial.

**Note** To address the point below, yes GANs devise their own "rules" to maximise win/loss WITHIN the rules already in place. However, they needs 'ground rules' otherwise how could it play a game of chess if it doesn't follow the rules of chess? How could it play an Atari video game if it does not know what a "joystick" is? The problem is without knowing the ground rules (what is biologically permitted), what's the validity of the subsequent rule set, i.e. predicted sequence in this case?

So a virus that can't attach to the host cellular receptor, isn't a virus because it will never be able to enter the cell.

**Further note**

The OP pointed out in the notes

> "The idea of style transfer has largely only been explored in image-based tasks, which we attribute in part to the specific nature of loss functions used for style transfer."

The author is talking about a variant of transfer learning. What the author is reminding us is that overall transferring-learning is the technique which powered the success of computer vision. They are saying they have replicated that in their GANs.

*Please note any further points should be separate questions*

---
Tags: rna, machine-learning, virus
---


# Question
Title: Converting Pairwise single- linkage clustering distance data to "newick" format
I have pairwise distance data for single linkage cluster and I would like to convert it to newick format.

I was unable to find a conversion algorithm.

I require newick format for visualization and annotation of a phylogeny.

Here is the output.

```
INFO  : ======================================
INFO  : Pairwise single linkage clustering
INFO  : ======================================
INFO  : Hierarchical Tree
INFO  : ======================================
INFO  : Node     Item 1   Item 2      Distance
INFO  :     0 :       53       52        0.001  5zxf  5z02
INFO  :    -1 :       64       60        0.002  7bj6  7bir
INFO  :    -2 :       62       -1        0.002  7biv
INFO  :    -3 :       34       26        0.002  4wt2  4ode
INFO  :    -4 :       65       -2        0.003  7bmg
INFO  :    -5 :       25       13        0.003  4occ  4erf
INFO  :    -6 :       29       -3        0.009  4ogt
INFO  :    -7 :       39       38        0.010  4zyi  4zyf
INFO  :    -8 :       50       48        0.011  5oc8  5ln2
INFO  :    -9 :       54       -8        0.012  6ggn
INFO  :   -10 :       42       41        0.015  5hmi  5hmh
INFO  :   -11 :       -4       61        0.020            7bit
INFO  :   -12 :       70       40        0.028  7qdq  5c5a
INFO  :   -13 :       47       68        0.032  5laz  7na3
INFO  :   -14 :       -9       45        0.035            5law
INFO  :   -15 :       49       59        0.038  5oai  6q9o
INFO  :   -16 :       32       -5        0.039  4qo4
INFO  :   -17 :      -15        8        0.044            3uu1
INFO  :   -18 :       23      -16        0.047  4oas
INFO  :   -19 :      -14       -7        0.053
INFO  :   -20 :       24      -18        0.057  4oba
INFO  :   -21 :      -19        6        0.057            3lbl
INFO  :   -22 :       43      -10        0.058  5hmk
INFO  :   -23 :       44       19        0.059  5lav  4jvr
INFO  :   -24 :      -23      -20        0.061
INFO  :   -25 :      -11       15        0.061            4hg7
INFO  :   -26 :       57      -24        0.062  6q9h
INFO  :   -27 :       27       -6        0.065  4odf
INFO  :   -28 :      -25      -17        0.067
INFO  :   -29 :       14      -28        0.067  4hbm
INFO  :   -30 :      -13      -21        0.067
INFO  :   -31 :      -27      -29        0.067
INFO  :   -32 :      -26      -30        0.068
INFO  :   -33 :       58       21        0.068  6q9l  4mdn
INFO  :   -34 :      -31      -32        0.070
INFO  :   -35 :       17        5        0.071  4jv9  3lbk
INFO  :   -36 :       30      -34        0.074  4ogv
INFO  :   -37 :      -33        2        0.077            1t4e
INFO  :   -38 :       35      -36        0.079  4zfi
INFO  :   -39 :       20      -38        0.086  4jwr
INFO  :   -40 :      -39      -37        0.094
INFO  :   -41 :        7      -40        0.106  3tj2
INFO  :   -42 :       51      -41        0.111  5trf
INFO  :   -43 :       69      -42        0.112  7na4
INFO  :   -44 :       56      -43        0.113  6q96
INFO  :   -45 :        4      -44        0.127  3jzk
INFO  :   -46 :       33      -45        0.129  4qoc
INFO  :   -47 :       22      -46        0.137  4mdq
INFO  :   -48 :      -12      -47        0.151
INFO  :   -49 :       37      -48        0.155  4zyc
INFO  :   -50 :       31      -49        0.156  4oq3
INFO  :   -51 :       12      -50        0.158  4ere
INFO  :   -52 :      -22       46        0.161            5lay
INFO  :   -53 :      -35      -51        0.182
INFO  :   -54 :       55      -53        0.186  6i29
INFO  :   -55 :       10      -54        0.224  3w69
INFO  :   -56 :       18      -55        0.238  4jve
INFO  :   -57 :        9      -56        0.265  3vzv
INFO  :   -58 :       28      -57        0.271  4ogn
INFO  :   -59 :      -58        1        1.000            1rv1
INFO  :   -60 :       66      -59        1.000  7na1
INFO  :   -61 :        0      -60        1.000
INFO  :   -62 :      -52      -61        1.000
INFO  :   -63 :       11      -62        1.000  4dij
INFO  :   -64 :       67      -63        1.000  7na2
INFO  :   -65 :       16      -64        1.000  4jv7
INFO  :   -66 :       36      -65        1.000  4zgk
INFO  :   -67 :        3      -66        1.000  2lzg

```

# Answer
> 2 votes
The problem with the tree format presented is that it does not conform to a known standard. Whilst its possible to write a parser I don't think its a good idea.

This output is hierarchical clustering, UPGMA, that sort of thing. Its not phylogenetics because this is a linkage tree.

Can you access the distance matrix prior to the clustering? If so the information below the line will help.

The output should provide the raw distance matrix prior to hierarchical clustering and will certainly have generated it in the calculation. BTW that particular clustering method isn't cool compared with nj (below) for complex reasons.

The calculation will be:

1. Perform pairwise linkage analysis
2. Perform cluster analysis on the resulting matrix

---

You need neighbor-joining to convert the distance matrix into a newick-format. There is no possibility of a direct conversion.

If you are using `R` its in the `ape` library via `treenk <- nj(M)` with `M` is the matrix. The treenk will be a format that can be read straight into `ggtree`

Briefly, a distance matrix requires a clustering algorithm and neighbor-joining is by far the best clustering method developed in phylogenetics.

The following is a phylogenetics example ... this is different from the example you have 'cause thats about linkage.

```
data(mynucleotides)
tr <- nj(dist.dna(mynucleotides))
plot(tr)

```

I your case the final line will be something like `ggtree(tr)` ... with various options.

---
Tags: phylogenetics, phylogeny, ggtree
---


# Question
Title: Turning GMT file into data frame so that third column is comma seperated
Similar to: How to convert data in gmt format to dataframe?

I have downloaded a GMT file from the GSEA database (https://www.gsea-msigdb.org/gsea/msigdb/download\_file.jsp?filePath=/msigdb/release/2022.1.Hs/c2.all.v2022.1.Hs.symbols.gmt)

I have imported this into R using the GSA.read.gmt tool from GSA as I want to do pathway analysis using GAUSS (https://github.com/diptavo/GAUSS/blob/master/example\_gmt.txt).

GAUSS needs the table to have the genes in a comma separated manner as well as there being headers in the table.

I have tried using the gsub function on the third column but it doesn't work, I feel this is because a GMT file is not the same as a dataframe.

Would anyone be able to help me turn this column into a comma separated one please?

Many thanks

# Answer
I used the unite function from tidyr but first I read the gmt file in as a data table using the read.delim function (with header set to False). Then I used unite:

```
library(tidyr)

gmt_file <- read.delim("~/path_to_file/file.gmt", header = F)
gmt_file <- gmt_file %>% unite(genes,3:ncol(gmt_file), sep=",")
colnames(gmt_file) <- c("pathway", "source", "genes")
write.table....

```
> 2 votes

---
Tags: r, pathway, dataframe
---


# Question
Title: How can I clean a PDB file using VMD or NAMD?
I have 3884 PDB files generated by pd2\_ca2main. PD2 CA2main is server for reconstructing the backbone of Calpha proteins using Gaussian mixture models.

All of these files have duplicate atoms or alternate coordinates.

How can I clean these PDB files using Visual Molecular Dynamics (VMD) or Scaleable Molecular Dynamics (NAMD) *****?

Can anyone tell me what command I should use in console mode?

***** note, this server may be unavailable

# Answer
You can use Chimera for this. The following code cleans your pdb file and write new pdb file.

```
import chimera
from DockPrep import prep
from chimera import runCommand, openModels

models = chimera.openModels.list(modelTypes=[chimera.Molecule])
prep(models)
runCommand("write format pdb 0 molecule.pdb")

```

You can run the code with

```
 chimera --nogui file.pdb code.py

```
> 1 votes

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