# Question
Title: How can I measure the distance between each pair of carbon alphas in overlayed protein structures?
I am trying to make a figure in which the protein structures are colored by the distance between the calphas in the overlayed structures.

I am comparing our cryoEM structure to that of a previous crystal structure and trying to show the major distance shifts in an easy to look at manner. I can get the RMSD using some online servers but it doesn't give me a good visual way of looking at the data.

# Answer
Ok, quite late.... my code stolen from https://pymolwiki.org/index.php/Spectrum as per in my comment:

> Not an expert here, but you could try to create a pdb with only c-alpha where B factor column has value according to bond lenght (proportional to it) and then color by b factor spectrum, I am sure there is some command about this : pymolwiki.org/index.php/Spectrum –

code using a fake peptide created with `cmd.fab('AFGVAL', 'peptide_1')` :

`peptide_1.pdb` :

```
ATOM      1  N   ALA     1      19.613   5.728  -6.281  1.00  0.00           N
ATOM      2  CA  ALA     1      20.289   7.022  -6.281  1.00  0.00           C
ATOM      3  C   ALA     1      21.789   6.848  -6.281  1.00  0.00           C
ATOM      4  O   ALA     1      22.320   5.731  -6.292  1.00  0.00           O
ATOM      5  CB  ALA     1      19.781   7.814  -5.064  1.00  0.00           C
ATOM      6  H   ALA     1      20.159   4.796  -6.281  1.00  0.00           H
ATOM      7  HA  ALA     1      20.021   7.561  -7.209  1.00  0.00           H
ATOM      8 1HB  ALA     1      18.685   7.964  -5.100  1.00  0.00           H
ATOM      9 2HB  ALA     1      20.005   7.300  -4.109  1.00  0.00           H
ATOM     10 3HB  ALA     1      20.237   8.819  -5.007  1.00  0.00           H
ATOM     11  N   PHE     2      22.540   7.897  -6.270  1.00  0.00           N
ATOM     12  CA  PHE     2      23.990   7.730  -6.270  1.00  0.00           C
ATOM     13  C   PHE     2      24.690   9.068  -6.255  1.00  0.00           C
ATOM     14  O   PHE     2      24.050  10.156  -6.244  1.00  0.00           O
ATOM     15  CB  PHE     2      24.437   6.942  -7.538  1.00  0.00           C
ATOM     16  CG  PHE     2      23.824   5.546  -7.729  1.00  0.00           C
ATOM     17  CD1 PHE     2      22.610   5.411  -8.413  1.00  0.00           C
ATOM     18  CD2 PHE     2      24.430   4.416  -7.174  1.00  0.00           C
ATOM     19  CE1 PHE     2      22.002   4.165  -8.523  1.00  0.00           C
ATOM     20  CE2 PHE     2      23.826   3.167  -7.294  1.00  0.00           C
ATOM     21  CZ  PHE     2      22.611   3.043  -7.966  1.00  0.00           C
ATOM     22  H   PHE     2      22.130   8.821  -6.260  1.00  0.00           H
ATOM     23  HA  PHE     2      24.287   7.177  -5.358  1.00  0.00           H
ATOM     24 2HB  PHE     2      25.541   6.842  -7.539  1.00  0.00           H
ATOM     25 3HB  PHE     2      24.233   7.549  -8.444  1.00  0.00           H
ATOM     26 1HD  PHE     2      22.116   6.283  -8.820  1.00  0.00           H
ATOM     27 2HD  PHE     2      25.362   4.509  -6.633  1.00  0.00           H
ATOM     28 1HE  PHE     2      21.050   4.077  -9.026  1.00  0.00           H
ATOM     29 2HE  PHE     2      24.295   2.298  -6.857  1.00  0.00           H
ATOM     30  HZ  PHE     2      22.134   2.079  -8.049  1.00  0.00           H
ATOM     31  N   GLY     3      25.980   9.099  -6.253  1.00  0.00           N
ATOM     32  CA  GLY     3      26.656  10.393  -6.239  1.00  0.00           C
ATOM     33  C   GLY     3      28.155  10.220  -6.239  1.00  0.00           C
ATOM     34  O   GLY     3      28.685   9.105  -6.251  1.00  0.00           O
ATOM     35  H   GLY     3      26.513   8.242  -6.261  1.00  0.00           H
ATOM     36 3HA  GLY     3      26.365  10.908  -5.371  0.00  0.00           H
ATOM     37  HA  GLY     3      26.378  10.953  -7.151  1.00  0.00           H
ATOM     38  N   VAL     4      28.906  11.269  -6.227  1.00  0.00           N
ATOM     39  CA  VAL     4      30.356  11.101  -6.228  1.00  0.00           C
ATOM     40  C   VAL     4      31.055  12.440  -6.212  1.00  0.00           C
ATOM     41  O   VAL     4      30.429  13.507  -6.201  1.00  0.00           O
ATOM     42  CB  VAL     4      30.798  10.251  -7.484  1.00  0.00           C
ATOM     43  CG1 VAL     4      32.326  10.136  -7.737  1.00  0.00           C
ATOM     44  CG2 VAL     4      30.276   8.794  -7.470  1.00  0.00           C
ATOM     45  H   VAL     4      28.497  12.192  -6.217  1.00  0.00           H
ATOM     46  HA  VAL     4      30.654  10.571  -5.305  1.00  0.00           H
ATOM     47  HB  VAL     4      30.366  10.744  -8.381  1.00  0.00           H
ATOM     48 1HG1 VAL     4      32.862   9.673  -6.888  1.00  0.00           H
ATOM     49 2HG1 VAL     4      32.559   9.537  -8.638  1.00  0.00           H
ATOM     50 3HG1 VAL     4      32.799  11.120  -7.921  1.00  0.00           H
ATOM     51 1HG2 VAL     4      29.177   8.752  -7.361  1.00  0.00           H
ATOM     52 2HG2 VAL     4      30.510   8.264  -8.413  1.00  0.00           H
ATOM     53 3HG2 VAL     4      30.710   8.201  -6.643  1.00  0.00           H
ATOM     54  N   ALA     5      32.345  12.471  -6.211  1.00  0.00           N
ATOM     55  CA  ALA     5      33.021  13.765  -6.196  1.00  0.00           C
ATOM     56  C   ALA     5      34.521  13.592  -6.197  1.00  0.00           C
ATOM     57  O   ALA     5      35.052  12.475  -6.219  1.00  0.00           O
ATOM     58  CB  ALA     5      32.511  14.544  -4.971  1.00  0.00           C
ATOM     59  H   ALA     5      32.879  11.614  -6.220  1.00  0.00           H
ATOM     60  HA  ALA     5      32.753  14.314  -7.119  1.00  0.00           H
ATOM     61 1HB  ALA     5      31.415  14.694  -5.007  1.00  0.00           H
ATOM     62 2HB  ALA     5      32.734  14.020  -4.022  1.00  0.00           H
ATOM     63 3HB  ALA     5      32.967  15.549  -4.903  1.00  0.00           H
ATOM     64  N   LEU     6      35.271  14.641  -6.173  1.00  0.00           N
ATOM     65  CA  LEU     6      36.722  14.473  -6.173  1.00  0.00           C
ATOM     66  C   LEU     6      37.420  15.811  -6.143  1.00  0.00           C
ATOM     67  O   LEU     6      36.792  16.881  -6.122  1.00  0.00           O
ATOM     68  CB  LEU     6      37.174  13.638  -7.404  1.00  0.00           C
ATOM     69  CG  LEU     6      36.580  12.214  -7.580  1.00  0.00           C
ATOM     70  CD1 LEU     6      37.231  11.514  -8.782  1.00  0.00           C
ATOM     71  CD2 LEU     6      36.743  11.343  -6.323  1.00  0.00           C
ATOM     72  H   LEU     6      34.862  15.564  -6.154  1.00  0.00           H
ATOM     73  HA  LEU     6      37.013  13.950  -5.245  1.00  0.00           H
ATOM     74 2HB  LEU     6      38.279  13.552  -7.384  1.00  0.00           H
ATOM     75 3HB  LEU     6      36.958  14.216  -8.325  1.00  0.00           H
ATOM     76  HG  LEU     6      35.495  12.316  -7.792  1.00  0.00           H
ATOM     77 1HD1 LEU     6      37.097  12.092  -9.716  1.00  0.00           H
ATOM     78 2HD1 LEU     6      38.318  11.367  -8.642  1.00  0.00           H
ATOM     79 3HD1 LEU     6      36.787  10.518  -8.962  1.00  0.00           H
ATOM     80 1HD2 LEU     6      37.799  11.259  -6.005  1.00  0.00           H
ATOM     81 2HD2 LEU     6      36.175  11.759  -5.470  1.00  0.00           H
ATOM     82 3HD2 LEU     6      36.357  10.318  -6.476  1.00  0.00           H
TER
END

```

`peptide_2_moved.pdb` , same pdb with some just carbon alpha moved by hand :

```
ATOM      1  N   ALA     1      19.613   5.728  -6.281  1.00  0.00           N
ATOM      2  CA  ALA     1      20.289   7.022  -6.281  1.00  0.00           C
ATOM      3  C   ALA     1      21.789   6.848  -6.281  1.00  0.00           C
ATOM      4  O   ALA     1      22.320   5.731  -6.292  1.00  0.00           O
ATOM      5  CB  ALA     1      19.781   7.814  -5.064  1.00  0.00           C
ATOM      6  H   ALA     1      20.159   4.796  -6.281  1.00  0.00           H
ATOM      7  HA  ALA     1      20.021   7.561  -7.209  1.00  0.00           H
ATOM      8 1HB  ALA     1      18.685   7.964  -5.100  1.00  0.00           H
ATOM      9 2HB  ALA     1      20.005   7.300  -4.109  1.00  0.00           H
ATOM     10 3HB  ALA     1      20.237   8.819  -5.007  1.00  0.00           H
ATOM     11  N   PHE     2      22.540   7.897  -6.270  1.00  0.00           N
ATOM     12  CA  PHE     2      24.369   7.337  -6.515  1.00  0.00           C
ATOM     13  C   PHE     2      24.690   9.068  -6.255  1.00  0.00           C
ATOM     14  O   PHE     2      24.050  10.156  -6.244  1.00  0.00           O
ATOM     15  CB  PHE     2      24.437   6.942  -7.538  1.00  0.00           C
ATOM     16  CG  PHE     2      23.824   5.546  -7.729  1.00  0.00           C
ATOM     17  CD1 PHE     2      22.610   5.411  -8.413  1.00  0.00           C
ATOM     18  CD2 PHE     2      24.430   4.416  -7.174  1.00  0.00           C
ATOM     19  CE1 PHE     2      22.002   4.165  -8.523  1.00  0.00           C
ATOM     20  CE2 PHE     2      23.826   3.167  -7.294  1.00  0.00           C
ATOM     21  CZ  PHE     2      22.611   3.043  -7.966  1.00  0.00           C
ATOM     22  H   PHE     2      22.130   8.821  -6.260  1.00  0.00           H
ATOM     23  HA  PHE     2      24.287   7.177  -5.358  1.00  0.00           H
ATOM     24 2HB  PHE     2      25.541   6.842  -7.539  1.00  0.00           H
ATOM     25 3HB  PHE     2      24.233   7.549  -8.444  1.00  0.00           H
ATOM     26 1HD  PHE     2      22.116   6.283  -8.820  1.00  0.00           H
ATOM     27 2HD  PHE     2      25.362   4.509  -6.633  1.00  0.00           H
ATOM     28 1HE  PHE     2      21.050   4.077  -9.026  1.00  0.00           H
ATOM     29 2HE  PHE     2      24.295   2.298  -6.857  1.00  0.00           H
ATOM     30  HZ  PHE     2      22.134   2.079  -8.049  1.00  0.00           H
ATOM     31  N   GLY     3      25.980   9.099  -6.253  1.00  0.00           N
ATOM     32  CA  GLY     3      26.656  10.393  -6.239  1.00  0.00           C
ATOM     33  C   GLY     3      28.155  10.220  -6.239  1.00  0.00           C
ATOM     34  O   GLY     3      28.685   9.105  -6.251  1.00  0.00           O
ATOM     35  H   GLY     3      26.513   8.242  -6.261  1.00  0.00           H
ATOM     36 3HA  GLY     3      26.365  10.908  -5.371  0.00  0.00           H
ATOM     37  HA  GLY     3      26.378  10.953  -7.151  1.00  0.00           H
ATOM     38  N   VAL     4      28.906  11.269  -6.227  1.00  0.00           N
ATOM     39  CA  VAL     4      30.719  10.841  -6.509  1.00  0.00           C
ATOM     40  C   VAL     4      31.055  12.440  -6.212  1.00  0.00           C
ATOM     41  O   VAL     4      30.429  13.507  -6.201  1.00  0.00           O
ATOM     42  CB  VAL     4      30.798  10.251  -7.484  1.00  0.00           C
ATOM     43  CG1 VAL     4      32.326  10.136  -7.737  1.00  0.00           C
ATOM     44  CG2 VAL     4      30.276   8.794  -7.470  1.00  0.00           C
ATOM     45  H   VAL     4      28.497  12.192  -6.217  1.00  0.00           H
ATOM     46  HA  VAL     4      30.654  10.571  -5.305  1.00  0.00           H
ATOM     47  HB  VAL     4      30.366  10.744  -8.381  1.00  0.00           H
ATOM     48 1HG1 VAL     4      32.862   9.673  -6.888  1.00  0.00           H
ATOM     49 2HG1 VAL     4      32.559   9.537  -8.638  1.00  0.00           H
ATOM     50 3HG1 VAL     4      32.799  11.120  -7.921  1.00  0.00           H
ATOM     51 1HG2 VAL     4      29.177   8.752  -7.361  1.00  0.00           H
ATOM     52 2HG2 VAL     4      30.510   8.264  -8.413  1.00  0.00           H
ATOM     53 3HG2 VAL     4      30.710   8.201  -6.643  1.00  0.00           H
ATOM     54  N   ALA     5      32.345  12.471  -6.211  1.00  0.00           N
ATOM     55  CA  ALA     5      32.965  13.942  -6.049  1.00  0.00           C
ATOM     56  C   ALA     5      34.521  13.592  -6.197  1.00  0.00           C
ATOM     57  O   ALA     5      35.052  12.475  -6.219  1.00  0.00           O
ATOM     58  CB  ALA     5      32.511  14.544  -4.971  1.00  0.00           C
ATOM     59  H   ALA     5      32.879  11.614  -6.220  1.00  0.00           H
ATOM     60  HA  ALA     5      32.753  14.314  -7.119  1.00  0.00           H
ATOM     61 1HB  ALA     5      31.415  14.694  -5.007  1.00  0.00           H
ATOM     62 2HB  ALA     5      32.734  14.020  -4.022  1.00  0.00           H
ATOM     63 3HB  ALA     5      32.967  15.549  -4.903  1.00  0.00           H
ATOM     64  N   LEU     6      35.271  14.641  -6.173  1.00  0.00           N
ATOM     65  CA  LEU     6      37.139  14.337  -6.385  1.00  0.00           C
ATOM     66  C   LEU     6      37.420  15.811  -6.143  1.00  0.00           C
ATOM     67  O   LEU     6      36.792  16.881  -6.122  1.00  0.00           O
ATOM     68  CB  LEU     6      37.174  13.638  -7.404  1.00  0.00           C
ATOM     69  CG  LEU     6      36.580  12.214  -7.580  1.00  0.00           C
ATOM     70  CD1 LEU     6      37.231  11.514  -8.782  1.00  0.00           C
ATOM     71  CD2 LEU     6      36.743  11.343  -6.323  1.00  0.00           C
ATOM     72  H   LEU     6      34.862  15.564  -6.154  1.00  0.00           H
ATOM     73  HA  LEU     6      37.013  13.950  -5.245  1.00  0.00           H
ATOM     74 2HB  LEU     6      38.279  13.552  -7.384  1.00  0.00           H
ATOM     75 3HB  LEU     6      36.958  14.216  -8.325  1.00  0.00           H
ATOM     76  HG  LEU     6      35.495  12.316  -7.792  1.00  0.00           H
ATOM     77 1HD1 LEU     6      37.097  12.092  -9.716  1.00  0.00           H
ATOM     78 2HD1 LEU     6      38.318  11.367  -8.642  1.00  0.00           H
ATOM     79 3HD1 LEU     6      36.787  10.518  -8.962  1.00  0.00           H
ATOM     80 1HD2 LEU     6      37.799  11.259  -6.005  1.00  0.00           H
ATOM     81 2HD2 LEU     6      36.175  11.759  -5.470  1.00  0.00           H
ATOM     82 3HD2 LEU     6      36.357  10.318  -6.476  1.00  0.00           H
TER
END

```

code :

```
import pymol

import math

from pymol import (

                    cmd ,
                    stored
                    )

print('########## PYMOL VERSION ##########################################')
print('         ',  cmd.get_version() )
print('###################################################################')

pymol.finish_launching()

cmd.set('internal_gui_width' ,  '500')

cmd.load('peptide_1.pdb' , 'peptide_1')

cmd.load('peptide_2_moved.pdb' , 'peptide_2_moved')

diff_len = lambda x,y : math.sqrt((x[0]-y[0])*(x[0]-y[0]) + (x[1]-y[1])*(x[1]-y[1]) + (x[2]-y[2])*(x[2]-y[2]))

cmd.select('peptide_1_CA' ,'peptide_1 and name CA')

cmd.select('peptide_2_moved_CA' ,'peptide_2_moved and name CA')

stored.peptide_1_CA_coord = []

stored.peptide_2_moved_CA_coord = []

cmd.iterate_state('1' , 'peptide_1_CA' , 'stored.peptide_1_CA_coord.append((x,y,z))')

cmd.iterate_state('1' , 'peptide_2_moved_CA' , 'stored.peptide_2_moved_CA_coord.append((x,y,z))')

print('stored.peptide_1_CA_coord : ' , stored.peptide_1_CA_coord)

print('stored.peptide_2_moved_CA_coord : ' , stored.peptide_2_moved_CA_coord)

stored.newB = list(map(lambda x,y: diff_len(x,y), stored.peptide_2_moved_CA_coord, stored.peptide_1_CA_coord))

print('\nstorednewB to be passed : ', stored.newB , '\n')

cmd.alter('peptide_2_moved_CA', 'b=stored.newB.pop(0)')

cmd.save('peptide_2_moved_CA.pdb' , 'peptide_2_moved_CA')

cmd.spectrum('b', 'rainbow_rev' , 'peptide_2_moved_CA' , '0' , '1')

cmd.color('lightblue' , 'peptide_1' )

cmd.set('cartoon_transparency' , '0.5', 'peptide_1')

cmd.alphatoall('peptide_2_moved')

cmd.save('peptide_2_moved_CA_alphatoall.pdb' , 'peptide_2_moved')

cmd.load('peptide_2_moved_CA_alphatoall.pdb' , 'peptide_2_moved_CA_alphatoall')

cmd.spectrum('b', 'rainbow_rev' , 'peptide_2_moved_CA_alphatoall' , '0' , '1')

```

output:

```
########## PYMOL VERSION ##########################################
          ('2.3.0'
...
stored.peptide_1_CA_coord :  [(20.288999557495117, 7.021999835968018, -6.281000137329102), (23.989999771118164, 7.730000019073486, -6.269999980926514), (26.6560001373291, 10.392999649047852, -6.238999843597412), (30.356000900268555, 11.10099983215332, -6.228000164031982), (33.020999908447266, 13.765000343322754, -6.196000099182129), (36.72200012207031, 14.472999572753906, -6.172999858856201)]
stored.peptide_2_moved_CA_coord :  [(20.288999557495117, 7.021999835968018, -6.281000137329102), (24.368999481201172, 7.336999893188477, -6.514999866485596), (26.6560001373291, 10.392999649047852, -6.238999843597412), (30.7189998626709, 10.840999603271484, -6.508999824523926), (32.96500015258789, 13.942000389099121, -6.048999786376953), (37.138999938964844, 14.336999893188477, -6.385000228881836)]

storednewB to be passed :  [0.0, 0.5984269572074781, 0.0, 0.5275691186189471, 0.23679966390593804, 0.48716415818795433]

```

showing in pymol window:

that is the CAs trace of moved (rainbow colored) vs original (lightblue \[lightgray out of ray to png\]) according to Palette https://pymolwiki.org/index.php/Palette\_Colorbars `rainbow_rev` :

that goes from red = 0 to BLUE = 1 as set out in

`cmd.spectrum('b', 'rainbow_rev' , 'peptide_2_moved_CA' , '0' , '1')`

againg transferred calculated CAs to all atom using `AlphaToAll` :

`cmd.alphatoall('peptide_2_moved')` (check it into code)

To use `spectrum` from min to max of the actual distance differences use in code above:

```
.....

stored.min = min(stored.newB)

stored.max = max(stored.newB)

print('\nmin of storednewB : ' , stored.min)

print('\nmax of storednewB : ' , stored.max)

print('\n')

cmd.alter('peptide_2_moved_CA', 'b=stored.newB.pop(0)')

cmd.save('peptide_2_moved_CA.pdb' , 'peptide_2_moved_CA')

cmd.spectrum('b', 'rainbow_rev' , 'peptide_2_moved_CA' , stored.min , stored.max)
....

```
> 1 votes

---
Tags: python, protein-structure, pymol
---


# Question
Title: Filtering cohorts in PySpark returns Py4JJavaError
I'm working with UK Biobank data. I made a cohort with the Cohort Browser and tried to extract the data with this code `df = participant.retrieve_fields(names=field_names, filter_sql = cohort.sql, coding_values = "replace", engine=dxdata.connect())`.

This is part of the text of the error: `Py4JJavaError: An error occurred while calling 04100.showString. org.apache.spark. SparkException: Exception thrown in awaitResult:`. The full error message is cut off.

However, if I remove the filter - `df = participant.retrieve_fields(names=field_names, coding_values = "replace", engine=dxdata.connect())` – there is no issue.

What kind of error is this and why would filtering for a cohort help?

# Answer
> 1 votes
It looks like the issue comes from one of the fields I was filtering for. When I remove that filter the code works again. Still need to investigate why that happened.

# Answer
> 0 votes
Generally remote biological servers, for example UCSC genome server, will not perform SQL queries for a free account. They'll allow downloads, but not queries. That is probably the same reason here.

The reason is simply these are often huge databases and queries can take up a lot of their computational resources, particularly if everyone is doing them. I don't know BioBank UK, but suggest its common practice.

If thats true the only difference with Biobank and UCSC is the UCSC server refusals are stated very clearly, "this operation ain't permitted for free accounts" (paraphrasing). What might be happening in an API error is not being picked up and relayed to you. Depends how the API is written so it's just reporting a timeout.

---
Tags: python
---


# Question
Title: How to subsample BAM based on read line number?
Consider a sorted, indexed, only mapped reads containing BAM file.

Is there a way to get a sub BAM based on read line numbers?

For example get a sub BAM file which contains read line numbers \[1, 3, 8, 70\] of the original BAM.

This can be done iteratively using a counter but its too slow. Is there a way to do this using samtools view or some tool like it.

Edit:

1. Size of the BAM file is ~7GB. Its human exome.
2. I'm using Julia (https://github.com/BioJulia/XAM.jl/) to iterate over the BAM and it takes me ~300 seconds to iterate over all the records. I feel this is too slow.
3. I'm running a function on each read of the BAM file and it decides which bucket the read has to go into. So this is specific subsampling.

I think `samtools view -` and `awk` route seems viable.

# Answer
> 2 votes
Samtools+awk:

```
samtools view -h aln.bam \
  | awk 'BEGIN{while((getline<"lineno.txt")>0)a[$1]=1}/^@/{++n;print}!/^@/&&a[NR-n]' \
  | samtools view -b --threads=4 > sub.bam

```

I would recommend mawk which is much faster than other awk.

I believe using Julia should be faster as it would not need to convert BAM to SAM and then convert back. Five min to read through a 7GB BAM file sounds about right. That is not slow.

---
Tags: sam, samtools, pysam
---


# Question
Title: Is there a data analysis software (free) or code that I can use to view normalized CDR3 amino acid length distributions?
Is there a data analysis software (free) or code that I can use to view normalized CDR3 amino acid length distributions?

Regarding code, preferably using Python or R.

EDIT: Added a picture of the non-normalized CDR3 amino acid length distribution. The issue is that ES0017\_metastasis has more samples than all other individuals which is why it has the highest peak at every amino acid length from 10-20. Was wondering if normalization was possible? Plot generated via Immunarch package in R.

EDIT2: Now able to access OriginPro and Graphpad Prism via university.

EDIT3:

EDIT4:  Here are the results, what do you think... The first chart is the log transformation (one extra as well, don't know how to delete). The bottom chart is the min-max transformation.

The x-axis refers to amino acid length. Each discrete block is an amino acid length, e.g. 12 amino acids long, etc.

# Answer
### Final update

Okay the answer below is based on the second graph from the top. Looking at the update the correct answer is closer to the first answer presented. **EDIT 3** is normalised data. It appears that the y-axis is probably natural log transformed (loge) or log2 and this transformation has normalised it.

Thus there is not need to normalise the data, the program has already done it in EDIT3 and in the very top graph now presented.

---

This looks impressive, thanks. I get it. There are two way to do this. If you are normalising the frequency (y-axis) is what you are wanting to normalise then either (in this instance its the wrong answer \[see below the ----\]):

* a log10 transformation (even Excel will log this data)
* a max-min transformation (thats in Python, likely in R - I dunno a general package maybe SPSS??, dunno)

You then test via the normal distribution using with Shapiro–Wilk test or Kolmogorov–Smirnov test

I suspect the max-min transformation will work. There are other ways in but it will involve discarding data that is \>3SD from the mean. Try the above and see how it works.

Note, a specialist statistician will request you perform a QQ plot, because the above tests of normalisation are not perfect.

This is not standardisation this is conformation to the normal distribution.

---

This could be about bin size if you have a discrete range for gene length and doesn't need to be very big - normally you need to try a few values, in this case its just 1 amino acid bin ... This can be used as a sliding window. It will pull the data into a normal distribution. That may not be very well explained, but I'm pretty sure that will work.

1 amino acid sounds weird but part the failure to see a normal distribution is a result of your calculation: you are splitting an amino acid into fractions and thats whats causing the irregularities in the data. Thus between values 10 and 20 on the x-axis (presumably size range difference) you should only have 20 bars - but you've got loads \> 100 which isn't possible.

Do that try SW and KS tests and you're sorted.
> 0 votes

---
Tags: rna-seq, sequence-analysis
---


# Question
Title: How to convert Foldseek .m8 alignment files to FASTA format?
I want to use sequence alignments found by Foldseek in PROSS. I know there are other ways to get alignments (also tried with results of HHPred), but I had in mind that Foldseek may find remote homologs that would still be structurally relevant.

Currently Foldseek only gives alignment results in m8 format, however PROSS only accepts FASTA alignment files.

I tried looking into BioPython but it seems it doesn't open m8 files, only m9. I also tried importing the alignment in Clustal X but it can't read it.

I could probably use Foldseek CURL API to get the alignment as json and make a script to write the results in FASTA format, but this seems like a heavy handed solution.

What would be other tools I could try to convert m8 to FASTA?

# Answer
## tl;dr

Use an awk script (below) or "one-liner"

```
awk '{name = "> "; field = 2; while ( $field + 0 != $field ) { name = name " " $field; field = field+1 } ; field = NF; if ( $field + 0 == $field ) { field = field -1 }; while ( $field + 0 != $field ) { field = field - 1 }; field = field - 1; seq = $field; gsub(/-/, "", seq); print name; print seq}' the_m8_file.m8 > the_fasta_file.fasta

```

Caveat emptor: I've tested this on my own Foldseek outputs and it works for my pedestrian purposes of multisequence alignment with COBALT. I haven't tested it with PROSS but expect it should work for you too.

## Long answer:

Using this definition of m8 (which is a modification of the BLAST m8 format):

https://github.com/steineggerlab/foldseek#output-search

And this definition of FASTA:

https://en.wikipedia.org/wiki/FASTA\_format

We see we want to take lines like this (I've added \[...\] here to show what I've cut off for brevity):

```
job.pdb_A   Some matched structure name 51.100  311 138 8   2   305 25  328 1.000   2.790E-35   1272    306 334 VDVSVSIFINKIYGVNTLEQTYKVDGYIVAQWTGKPR--KTPGDKPLIVENTQIERWINNGLWVPALEFINVVGSPDTGNKRLMLFPDGRVIYNARFLGSFSNDMDFRLFPFDRQQFVLELEPFSYNNQQLRF--SDIQVYTENIDN[...] -27.488,-8.287,-35.993,-26.838,-4.500,-36.421,-25.278,-3.958,-32.938,-22.398,-6.436,-33.555,-21.375,-4.861,-36.892,-21.373,-1.326,-35.385,-19.237,-2.585,--18.044,-1.289,8.073,-14.345,-2.072,7.367,-15.099,-5.282,5.328,-16.872,-3.358,2.522,-16.876,-5.018,-0.77 [...] MIKKEILSAIASLVWFCSPALYAQYVVDTSVSINKIYGVNTVDQTYKVDGYLVTSWDEPDPTYRPASGERVYENHHADQLMNDGLWVPAFEFINIIGQRQTANRRVVISADGRVTYNERFQGLFTTTMDFRQFPFDYQTFVLMMEPFSYEQKQL[...]   1918945 Vibrio spartinae

```

and map them to entries like:

```
>some matched structure name
MIKKEILSAIASLVWFCSPALYAQYVVDTSVSINKIYGVNTVDQTYKVDGYLVTSWDEPDPTYRPASGERVYENHHADQLMNDGLWVPAFEFINIIGQRQTANRRVVISADGRVTYNERFQGLFTTTMDFRQFPFDYQTFVLMMEPFSYEQKQL[...]

```

We can do this with an awk script (The awk manual is your friend here https://www.gnu.org/software/gawk/manual/html\_node/Long.html). Put the following in a file somewhere handy and call it using `awk -f m8_to_fasta.awk my_m8_file.m8 > my_fasta_file.fasta` you may want to combine this with head, tail, and/or grep to reduce the number of sequences that end up in your output.

```
#! /bin/awk
{
        name = "> "
        # start on the second field, foldseek will put the job name in $1
        field = 2
        # Check each field until we see numbers
        while ( $field + 0 != $field ) {
           name = name " " $field
           field = field+1
        }
        # name now contains the name string which is the first line of a fasta entry

        field = NF
        # Now find the matched sequence
        # Sometimes an entry ENDS in a number, so we strip that
        if ( $field + 0 == $field ) { field = field -1 }
        while ( $field + 0 != $field ) {
            field = field - 1
        }
        # field is now pointing at some undocumented number. Decrement once more to get the sequence
        field = field - 1
        # field is now pointing at the sequence
        seq = $field

        #Now I/we want to sanitize the sequence to remove any gaps
        gsub(/-/, "", seq)

        print name
        print seq
    }

```

Or we can coerce this into a monstrous one-liner:

```
 awk '{name = "> "; field = 2; while ( $field + 0 != $field ) { name = name " " $field; field = field+1 } ; field = NF; if ( $field + 0 == $field ) { field = field -1 }; while ( $field + 0 != $field ) { field = field - 1 }; field = field - 1; seq = $field; gsub(/-/, "", seq); print name; print seq}' the_m8_file.m8 > the_fasta_file.fasta

```
> 2 votes

---
Tags: sequence-alignment, protein-structure, format-conversion
---


# Question
Title: What is the required input peak file format for IDR (Irreproducible Discovery Rate)
I have managed to install IDR (https://github.com/nboley/idr) correctly (idr --version returns 2.0.3), but I can't get it to run on my peak files (\*.bed) without crashing.

My peak files are generated using HOMER (findPeaks) and they are converted to proper \*.bed format via HOMER/pos2bed.pl. But when I run idr on my \*.bed peak files as described in the IDR guide page like below:

idr --samples peak1.bed peak2.bed

it always returns an error inside the function bed\_load() : signal = float(data\[signal\_index\]) IndexError: list index out of range

I suspect it's related to peak files (likely columns) not being in the format the program idr.py is expecting, and I couldn't find much about the required "input" file format in the idr help pages.

Can someone provide examples of an idr.py friendly peak file, maybe just a couple of rows with the proper/expected column types and numbers? Thank you in advance.

I am including the guesstimate guidance from ChatGPT, which is not necessarily correct, but maybe a pointer in the right direction - it's suggesting the IDR format has to be consistent with narrowPeak formats having 10 columns, and the actual column that IDR uses in it's analysis is column 7 - which it calls it the "signal value". But what exactly is this measure (the peak strength, the peak height etc.)? Any insight would be appreciated? Also, on where this type of information may have been discussed, seeing how critical it is to running the tool successfully.

The IDR readme on usage describes very clearly what the OUTPUT format is, but is not clear what the INPUT format is (column descriptions) - which is what one needs to successfully form the INPUT and run IDR. My typical HOMER/pos2bed.pl converted row is less than 8 columns which at this point, I am certain is the cause of IDR failures. I am attempting to figure out which column IDR actually uses, and which are not essential for it's considerations (they just get duplicated in the output).

Here are the error messages I get:

https://github.com/nboley/idr/issues/41

# Answer
> 2 votes
As mentioned in the IDR readme on usage, there are examples in the IDR github repository here. Here's one of them:

```
chr16   31085416    31086009    .   0   .   362.546677850757    -1  4.12368834166759    250
chr12   109490129   109490582   .   0   .   360.329496957805    -1  4.12368834166759    215
chr13   41345146    41345561    .   0   .   358.624215585062    -1  4.12368834166759    207
chr19   50083722    50084001    .   0   .   266.603941224671    -1  4.5777559727315 186

```

From the looks of it, it seems that the output is a line-filtered version of the input format (i.e. output and input formats are identical).

I have no idea how IDR works (never used it before), but the error you've added to your post suggests there's a range mismatch. My guess based on that error message is that the two input files don't have overlapping specified regions.

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


# Question
Title: obtaining a fragment data file bed from bam
Im trying to follow this pipeline https://github.com/epifluidlab/cragr on cfDNA WGS. I have bam files from paired end reads aligned with bwa mem.

The workflow requires bgzipped and indexed fragment data files as input.

I have been searching how can I obtain BED files from fragments coordinates instead of reads coordinates, and came up with no answer.

Does someone know how can I obtain such bed file?

# Answer
There's an example in the README (https://github.com/epifluidlab/cragr#example-datasets):

> For you information, below are steps from the BAM file to the fragment file:

```
samtools sort -n -o qsorted.bam -@ 4 sample.bam
samtools view -h -f 3 -F 3852 -G 48 --incl-flags 48 qsorted.bam |\
  bamToBed -bedpe -mate1 -i stdin |\
  awk -F'\t' -v OFS="\t" '{if ($1!=$4) next; if ($9=="+") {s=$2;e=$6} else {s=$5;e=$3} if (e>s) print $1,s,e,$7,$8,$9}' |\
  sort -k1,1V -k2,2n |\
  bgzip > sample.frag.bed.gz
tabix -p bed sample.frag.bed.gz

```
> 2 votes

---
Tags: sam, bed, bedtools, wgs, epigenetics
---


title: stats.sh error in BBmap package

text:
# Question
Title: stats.sh error in BBmap package
I am trying to calculate the N50 value from the assembled FASTA file. I used stats.sh from the BBmap package. I executed the following command

```
./stats.sh  in=/Home/Downloads/BBMap_39.01/bbmap/resources/contig.fasta

```

I got the following error.

```
Exception in thread "main" java.lang.RuntimeException: Input file does not appear to be valid: /Home/Downloads/BBMap_39.01/bbmap/resources/contig.fasta
    at jgi.AssemblyStats2.process(AssemblyStats2.java:248)
    at jgi.AssemblyStats2.main(AssemblyStats2.java:40)

```

How to resolve this error?

The sequences are shown below. My file also contains some lengthy sequences.

```
>NODE_1_length_3013_cov_4.540906
TGGGGGATCATGTAACTCGCCTTGATCGTTGGGAACCGGAGCTGAATGAAGCCATACCAA
ACGACGAGCGTGACACCACGATGCCTGTAGCAATGGCAACAACGTTGCGCAAACTATTAA
CTGGCGAACTACTTACTCTAGCTTCCCGGCAACAATTAATAGACTGGATGGAGGCGGATA
AAGTTGCAGGACCACTTCTGCGCTCGGCCCTTCCGGCTGGCTGGTTTATTGCTGATAAAT
CTGGAGCCGGTGAGCGTGGGTCTCGCGGTATCATTGCAGCACTGGGGCCAGATGGTAAGC
CCTCCCGTATCGTAGTTATCTACACGACGGGGAGTCAGGCAACTATGGATGAACGAAATA
GACAGATCGCTGAGATAGGTGCCTCACTGATTAAGCATTGGTAACTGTCAGACCAAGTTT
>NODE_2_length_2292_cov_16.637908
CAAAGGAATAATAATTGACTGATGATTGGTTCATTAATTACTTCATTAGTTAATAGATTG
CAGTCTTTATAAATAAGTTTATAGTTGTATTCCACAGAATTTCTGATCCTGTCTCTGCCC
TACCCAATGTTTTATCATTAAACTTCATTTACATAAGATAGTATGTTCAACAAATATAAA
GATTTGCTGGAATGGAAGGCACAGTAAGATATTTCACAGGATAAAATTAAAGTAGATGTC
AACAGAGTAAATCAAAGTGAAATTTAATAGAGTTAAACAGTAAGTGTTTTAAATACCAGG
ATTAAATTCCATGAAGCTGCAGCTTTTAAAAAAAGGCACAAAGTGGAATAAGAAAAAATG
ACTTTCAGAAAGTTAATGTAAAGAAAGCTCACGGTTTTCATTGGCTGCCTGCTCAGTTTG

```

# Answer
> 1 votes
```
./stats.sh in=notExistingFile.fa
Exception in thread "main" java.lang.RuntimeException: Input file does not appear to be valid: notExistingFile.fa
        at jgi.AssemblyStats2.process(AssemblyStats2.java:248)
        at jgi.AssemblyStats2.main(AssemblyStats2.java:40)

```

The error comes when the file does not exist. My best guess is that it is a typo on your side. Is the path really `/Home` with capital "H"? Most OS like Linux and macOS use lowercase. Anyway, make sure the file exists.

---
Tags: fasta, sequence-analysis, linux
---


# Question
Title: what does 8 stands for in bwa mem command line
What does this command implies in bwa mem

```
bwa mem -t 8 index/referencegenome.fna read1.fastq.gz read2.fastq.gz > output.sam"""

```

What does `-t 8` stands for here?

Thank you

# Answer
> 0 votes
`-t` is an important parameter here and is important to understand. bwa is a CPU calculation and the more CPUs the job (calculation) is spread across the faster the calculation will take to complete and the less chance of the calculation exceeding the RAM capacity of the machine. "Exceeding RAM capacity" (I don't know the technical term), is when the process "grinds to a halt" and never completes.

Basically what can happen is the indexing step instead of taking ~3 hours, could easily be still running after 24 hours if the machine has little RAM and `-t` isn't set sufficiently high.

One CPU is known as "one core".

For a regular desktop if you set `-t X` where X is number of CPUs - 1 in my opinion that is good practice. So if the machine has 8 CPUs then X = 7. If you are working on a cluster, e.g. GCP or AWS, then X is the maximum number of cores (CPUs) available.

# Answer
> 1 votes
When you have questions on software then the first thing is to do your due diligence and read the manual:

https://bio-bwa.sourceforge.net/bwa.shtml

Then simply checking for the "-t" argument will tell you that it is the number of threads used for the alignment process.

---
Tags: bwa, command-line
---


# Question
Title: What's the best way to download data from the SRA? Is it really this slow?
I'm trying to download three WGS datasets from the SRA that are each between 60 and 100GB in size. So far I've tried:

* Fetching the .sra files directly from NCBI's ftp site
* Fetching the .sra files directly using the aspera command line (`ascp`)
* Using the SRA toolkit's `fastqdump` and `samdump` tools

It's excruciatingly slow. I've had three `fastqdump` processes running in parallel now for approximately 18 hours. They're running on a large AWS instance in the US east (Virginia) region, which I figure is about as close to NCBI as I can get. In 18 hours they've downloaded a total of 33GB of data. By my calculation that's ~500kb/s. They do appear to still be running - the fastq files continue to grow and their timestamps continue to update.

At this rate it's going to take me days or weeks just to download the datasets. Surely the SRA must be capable of moving data at higher rates that this? I've also looked, and unfortunately the datasets I'm interested have not been mirrored out to ENA or the Japanese archive, so it looks like I'm stuck working with the SRA.

Is there a better way to fetch this data that wouldn't take multiple days?

# Answer
Proximity to NCBI may not necessarily give you the fastest transfer speed. AWS may be deliberately throttling the Internet connection to limit the likelihood that people will use it for undesirable things. There's a chance that a home network might be faster, but you're likely to get the fastest connection to NCBI by using an academic system that is linked to NCBI via a research network.

Another possibility is using Aspera for downloads. This is unlikely to help if bandwidth is being throttled, but it might help if there's a bit of congestion through the regular methods:

https://www.ncbi.nlm.nih.gov/public/

NCBI also has an online book about best practises for downloading data from their servers.

On a related note, just in case someone sees this and EBI/ENA is an option, there's a great guide for how to do file transfer using Aspera on the EBI web site:

https://www.ebi.ac.uk/ena/browse/read-download#downloading\_files\_aspera

> Your command should look similar to this on Unix:

```
ascp -QT -l 300m -i <aspera connect installation directory>/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:<file or files to download> <download location>

```

In my case, I've just started downloading some files from a MinION sequencing run. The estimated completion time via standard FTP was 12 hours for about 32GB of data; `ascp` has reduced that estimated download time to about an hour. Here's the command I used for downloading:

```
ascp -QT -l 300m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/ERA932/ERA932268/oxfordnanopore_native/20160804_Mock.tar.gz .

```
> 12 votes

# Answer
In general, the best way to download SRA data is: don't download from SRA. However, as ENA has not be sync'd yet, I would recommend to download from SRA ftp and then convert to fastq locally. You can find files in the SRA format here. Downloading and then converting locally is much faster than direct retrieval from NCBI for some mysterious reasons.
> 6 votes

# Answer
By far the fastest method in my experience has been to use the SRAdb library in R. For most entries, you can download fastq files directly. Some older experiments don't have them, but I've still found it much faster to download SRA files via `getSRAfile()` and then to convert them using `fastqdump` than to use `fastqdump` directly.
> 3 votes

# Answer
NCBI is moving slowly to cloud (AWS/Google) for hosting. As such the tradiational FTP urls do not always work. My solution involves using pysradb.

An entire project can be downloaded by:

```
pysradb download -p <SRP_ID>

```

A Google Colab notebook with metadata retrieval commands is here.
> 2 votes

# Answer
Another possibility in 2020 is to use fasterq-dump which is promoted to be faster than previous `fastq-dump`.
> 1 votes

# Answer
I recommend using fastq-dl: downloads from either ENA or SRA, multithreaded, automatic retries.
> 1 votes

---
Tags: data-download
---


# Question
Title: Why doesn't bedtools getfasta pick up the first base?
I am using bedtools getfasta to convert genomic coordinates to sequences. I feed it a bedfile containing say chr1 s.start s.end. When I check the sequence output, it does not include the s.start base. For example if the bedfile entry is chr1 23456 23678, it gets the sequence from 23457-23678. Is this a feature and I should use a tag to make it include the first base? The command I am using is `bedtools getfasta -fi genome.fa -fo gene.fa -bed input.bed` I tried this in bedtools versions 2.27.1 and 2.30.1

# Answer
BED files are 0-based. I link this old but gold biostars thread for more background. Essentially, from normal counting with start 1 you subtract 1 from the start coordinate.

For example in the sequence `ATCGGG` the `TCGGG` substring would have "normal" (aka human-readable) coordinates 2-6, but in BED format that would be:

```
chrSomething    1    6

```

Here is the bedtools "proof":

```
$ cat foo.bed
chrSomething    1       6

$ cat foo.fa
>chrSomething
ATCGGG

$ bedtools getfasta -fi foo.fa -bed foo.bed
>chrSomething:1-6
TCGGG

```
> 3 votes

---
Tags: genome, bed, bedtools
---


# Question
Title: What are the output files of RNA-Seq from facility?
*This question was also asked on Reddit*

I am new in our lab and I am going to do bulk RNA-Seq. What type of files will we get from the company (Genewiz)? Will it be a bunch of Fastq files? or they give something like bam files?

# Answer
This cannot be predicted unless you have interacted with them previously. Most commonly a sequencing facility will send raw data as gzipped fastq files, often together with md5 checksums so the user can check file integrity after transfer/download from them.

We worked with a facility before who sent data in unaligned bam format, which then for typical workflows first had to be converted back to fastq.

After all, you're a paying customer so feel free to shoot them an email and ask whatever question comes up. Ask what exactly they will provide so that you can prepare your workflows accordingly. Most likely it will be fastq.gz.
> 1 votes

---
Tags: rna-seq, ngs, fastq
---


# Question
Title: How evaluate biclustering algorithm over real datasets containing labels?
I have a set of real data containing **labels**. However, in some biclustering algorithms (e.g., Cheng and Church Algorithm), originally authors applied gene expression data like Yeast without having labels. My question is that what are appropriate *evaluation measures* in this scenario?

# Answer
> 0 votes
Briefly evaluation measures, which I presume means "features", are determined using supervised learning using a technique called "feature selection". Unsupervised learning generally is used to identify categories for a supervised learning algorithm. The orthodoxy has changed.

*"Unsupervised learning feature selection"*

1. To identify suitable features using unsupervised learning *a priori* then genome location which must have a major outcome, thereafter it would be a case of deleting a given feature and re-evaluating the clusters. If the clusters disintegrate then that is a key feature (not very cool but would work for datasets with limited numbers of features).
2. Weighted unsupervised learning is popular, this is basically "feature selection", i.e. which features contribute to the clustering. This is basically the formal version of point 1. Personally I have had little success with this approach, i.e. the results were not informative, but it does work. The results are not discriminatory between different conditions, which is what I was seeking.

Regarding labels, labels are not used in unsupervised learning approach - except weighted unsupervised learning - hence their presence or absence usually makes little difference. The only "label" used is the sample type, i.e. the "id" and not the "feature". Without intrinsically understanding the experiment it is difficult to forward an accurate *a priori* assessment.

---
Tags: gene-expression, machine-learning, algorithms
---


title: ncbi fastq dump error in loop

text:
# Question
Title: ncbi fastq dump error in loop
I have a very basic objective which I want to give list of id to a ncbi-srafastq tool kit which will prefetch the id which is in .sra file extension and then it will run fastq-dump on the file to split the files into fastq files, single file if there is single end file or paired end file if its paired end or two set of files.

This is my code

```
#!/bin/bash

input=$1

a=$(basename "$input")
id=${a/.id/}

# Download SRA data
prefetch "$id"

# Extract SRA ID from the downloaded data
#sraID=$(ls -1 | grep -v run.sh)
sraID=$(ls -1 | grep -v run.sh | sed 's/\.sra$//')
# Validate the downloaded SRA data
vdb-validate "$sraID"

# Dump fastq data from the SRA file
fastq-dump "$sraID"  --gzip

# Concatenate fastq files for read1, read2, and unpaired reads (if available)
for s in "$sraID"
do
  echo "$s"
  if [ "$s" != "$id" ]
  then
    [ -f "${s}_1.fastq" ] && cat "${s}_1.fastq" >> "${id}_R1.fastq"
    [ -f "${s}_2.fastq" ] && cat "${s}_2.fastq" >> "${id}_R2.fastq"
    [ -f "${s}.fastq" ] && cat "${s}.fastq" >> "${id}.fastq"
  fi
done

# Compress fastq files using pigz
#for i in "${id}"*.fastq
#do
  #echo "$i"
 # pigz "$i" & sleep 1
#done
#wait

# Remove the original uncompressed fastq file if both read1 and read2 files are present
[ -f "${id}_R1.fastq.gz" ] && [ -f "${id}_R2.fastq.gz" ] && [ -f "${id}.fastq.gz" ] && rm "${id}.fastq.gz"

# Test deleting intermediate folder
echo "Space before tmp folder delete"
echo "Removing tmp folder $sraID "
ls *
rm -rf "$sraID"*
ls *
echo 'done'

```

I run using as small example file as such

```
./run.sh SRR12332899

```

The error I get is this

```
2023-08-05T12:37:32 prefetch.3.0.6: Current preference is set to retrieve SRA Normalized Format files with full base quality scores.
2023-08-05T12:37:33 prefetch.3.0.6: 1) Downloading 'SRR12332899'...
2023-08-05T12:37:33 prefetch.3.0.6: SRA Normalized Format file is being retrieved, if this is different from your preference, it may be due to current file availability.
2023-08-05T12:37:33 prefetch.3.0.6:  Downloading via HTTPS...
2023-08-05T12:40:48 prefetch.3.0.6:  HTTPS download succeed
2023-08-05T12:40:48 prefetch.3.0.6:  'SRR12332899' is valid
2023-08-05T12:40:48 prefetch.3.0.6: 1) 'SRR12332899' was downloaded successfully
2023-08-05T12:40:48 prefetch.3.0.6: 'SRR12332899' has 0 unresolved dependencies
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'SRR12332899.sra' metadata: md5 ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Table 'SEQUENCE' metadata: md5 ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'ALTREAD': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'QUALITY': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'READ': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'SPOT_GROUP': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'sra/SRR12332899.sra' contains only unaligned reads
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'SRR12332899.sra' is consistent
2023-08-05T12:40:49 fastq-dump.3.0.6 err: item not found while constructing within virtual database module - the path 'sra' cannot be opened as database or table
fastq-dump quit with error code 3
sra
Space before tmp folder delete
Removing tmp folder sra
run.sh

sra:
SRR12332899.sra
run.sh
done

```

If I do just

`prefetch SRR12332899` and then `fastq-dump SRR12332899`

I get a SRR12332899.fastq SRR12332899.sra

This is for one file since I want to download list of files using the script which is not working.

Any suggestion would be really helpful

# Answer
> 3 votes
The main problem is that your script is inadvertently treating a set of file and directory names as a single SRA accession number, all in one variable. You should update it to work with each accession separately and make sure they really are accessions.

---

On this line you're taking all the names of things in your working directory (minus the script), each on a separate line, and trimming .sra off the end of each line:

```
sraID=$(ls -1 | grep -v run.sh | sed 's/\.sra$//')

```

In your case that list was apparently just a directory named "sra". (Even though I don't have that directory on my end, you can see it in your output, from the first of two `ls *` at the end, where you have `sra/SRR12332899.sra`.) So later on it tries to `fastq-dump sra`, which gives you that "cannot be opened as database or table" error. (`vdb-validate` can take a directory as input, so it doesn't complain.) If I run your script from scratch in an empty directory with SRR12332899, it actually seems to run OK, aside from removing everything it creates at the end. But if you run it with a second accession you see the problem, since it then tries things like `vdb-validate "several accessions and filenames separated by newlines"`.

I think a proper fix would involve a bunch of changes, but here are some tips, at least:

* use subdirectories to keep things organized so you don't have to juggle multiple types of file in one place
* `set -x` at the start of a script will show you each command as it happens, which can be very helpful for troubleshooting
* be careful with quotes and whitespace in bash. If you want to loop over each word in a variable you can do `for thing in $var; do ...; done` (unquoted) but if you run into filenames with whitespace, it breaks them up. Another pattern is `command | while read thing; do ...; done` which will go line-by-line. (In that case if you were working with a variable storing the command output, you might do something like `echo "$output" | while...`, quoted that time, since you'd want to preserve the newlines.)
* `ls *` will expand each file/directory name as separate arguments (so bash really ends up calling `ls "sra" ...whatever else...`), and when ls sees directory names is lists their contents too. If you just want to show all the files/directories you can leave off the asterisk altogether, or if you wanted to give directory names but not have it list their contents, you can use `ls -d`.
* you don't necessarily have to call `ls` to store the names as a block of text. You can loop over each filename, including those with spaces, directly in bash (e.g. `for thing in *; do echo "path: '$thing'"; done`)
* be careful with the rm command, especially combined with wildcards, and especially with `-rf`. If `sraID` happened to be empty somehow, think of what `rm -rf "$sraID"*` would actually mean. (And yes, scripts can remove themselves!)

---
Tags: fastq, shell, ncbi
---


# Question
Title: Interpretation of clusterProfiler results
I am using clusterProfiler for differential gene expression using this code: https://yulab-smu.top/biomedical-knowledge-mining-book/enrichplot.html

I have a question about the third line of code:

```
 de <- names(geneList)[abs(geneList) > 2]

```

Here all gene names are taken with absolute fold change \>2. As an example, if in a hypothetical extreme case you have 100 downregulated genes (log2 fold change \<-2) and zero upregulated genes (log2 fold change \>2), then you may incorrectly think that identified pathways are activated, but they are actually repressed. How is this accounted for? Thank you.

# Answer
The scenario you described would be reported by the code because:

```
... [abs(geneList) > 2]

```

The `abs` will change any negative value to positive and retain all positive values as positive. Therefore values \<-2 will be reported as well as values \>2. `abs` = absolute and to take the absolute value of a number is exactly the operation just described.

---

From comments. Just remove `abs`

```
negative <- names(geneList)[geneList < -2]
positive <- names(geneList)[geneList > 2]

```
> 1 votes

---
Tags: r, differential-expression, gene-expression, pathway, clusterprofiler
---


title: Reference genome (bwa -mem)

text:
# Question
Title: Reference genome (bwa -mem)
I want to do mapping with bwa mem. I downloaded a reference genome from ncbi database which is (Danaus plexippus). I got one `ncbi_database.zip` file which contain 2 folders and 3 files of doc,assembly\_report\_json file. I read one file named `read.md` which came after downloading but its not helping. What should be done in this case?

I don't have an issue with using bwa-mem; the problem is with the reference genome from NCBI. I downloaded it from https://ncbi.nlm.nih.gov/datasets/genome/GCF\_009731565.1 but after that I checked the downloaded file and I got 5 folders into it and I'm confused which one to use. 2 files have `.fna` extension and the rest are reports that are in the json file format.

# Answer
I downloaded genome from https://www.ncbi.nlm.nih.gov/datasets/genome/GCF\_009731565.1/ I got 5 files in it so if you are confused which one to use in that so GCF.fna having genomic data with fna extension is used as an reference genome if you are also new to this things this might be useful.

Thank you
> 0 votes

# Answer
What you need is the genome sequence in fasta format, typically suffixed `.fa`, for example from https://www.ncbi.nlm.nih.gov/datasets/genome/GCF\_009731565.1/

This file needs to be indexed via `bwa index`.
> 1 votes

---
Tags: read-mapping, reference-genome, bwa, ncbi
---


# Question
Title: Using CSV files continaing scRNA-seq count data from GEO
I'm completely new to scRNA-seq analysis (and to most of Bioinformatics as well) so I apologize if this was asked a million times before.

I am trying to get the single cell expression data of lymph node and bone marrow from this GEO entry

In Samples, there are multiple organs and accessing any of them gets me a CSV file compressed in gz format. The CSV file appears to contain the barcodes and the counts, but I do not know where to go from there. I need to annotate each cell type and I am not sure how. Seurat tutorials indicate that I have to have a few other file types which do not seem to be available in this repository. Does this mean I cannot use it, or is there a way to get this data to count/cell type just using what's available? I imagine that I don't need to perform any alignment since the rows already contain gene names.

Any and all guidance would be very much appreciated. Thank you so much in advance.

# Answer
> 1 votes
Please understand that your question is open-ended and essentially you're asking for a complete guidance through your analysis which is not feasible. I recommend to go through either the Seurat tutorials or the Bioconductor single-cell workflow to learn the basics, such as typical workflows, formats of required input files etc.

With this knowledge you can then start finding out which sort of data are provided in GEO and start getting your analysis going. You're welcome to ask specific questions here, but please be precise and not open-ended.

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


# Question
Title: What is the meaning of "Filtering allele frequency" in the gnomAD database?
Recent gnomAD versions include a "filtering allele frequency" which tells you when a variant can be safely adjudged not to be disease-causing. Unfortunately, I'm having trouble making sense of the definition, which is

> "Technically, this is the highest disease-specific maximum credible population AF for which the observed AC is not compatible with pathogenicity. More practically, If the filter allele frequency of a variant is above the maximum credible population AF for a condition of interest, then that variant should be filtered (ie not considered a candidate causative variant). See http://cardiodb.org/allelefrequencyapp/ and Whiffin et al. 2017 for additional information."

I interpret this statement as "If

> gnomAD\_WG\_AF\_{population} \> gnomAD\_WG\_FAF95\_{population},

then the variant should be filtered."

(Here "FAF95" stands for "filtering allele frequency at 95% confidence".) However, the above inequality seems to hold for every variant in gnomAD, so I'm clearly missing something. On the other hand, perhaps the statement "maximum credible population AF for a condition of interest" means "proportion of people in the population who have the condition," but in that case the use of the term AF doesn't make sense. Does anyone have a better understanding of this FAF95 number than I do?

# Answer
My understanding is from Whiffin et al. 2017.

Imagine a dominant mendelian disease A with complete penetrance, caused only by one allele, which is currently unknown. In your search for that unknown allele, a "credible" candidate variants is one whose population allele frequency is lower than the prevalence of that disease in the population. In that sense,

$$\text{max. credible. AF}\_\text{A} = \text{prevalence}$$

In a more realistic scenario, a similar disease B can be caused by more than one allele. A credible candidate variant would therefore have a lower expected allele frequency, because only some proportion of diseases in the population should be attributable to that variant (the "allelic contribution" of that variant),

$$\text{max. credible AF}\_\text{B} = \text{prevalence} \times \text{max. allelic contrib.}$$

Finally, consider a similar disease C with incomplete penetrance. Intuitively, a credible candidate variant would be allowed to have higher population allele frequency than the disease prevalence, due to the presence of unaffected carriers,

$$\text{max. credible AF}\_\text{C} = \text{prevalence} \times \text{max. allelic contrib.} \times 1/\text{penetrance}$$

Thus far, we have talked about the maximum credible *population* allele frequency. Population reference databases like gnomAD only have a sample of the population. Whiffin and colleagues observe that one can set a value for the maximum credible allele frequency, and subsequently compute by modelling a Poisson distribution, the one-tailed 95% confidence interval for the expected allele count in a sample of a given size.

Conversely, they suggest, you could work backwards from an allele count and find the highest possible value for the maximum credible population allele frequency that gives the allele count as the 95% confidence interval. I believe this is the FAF95 recorded in gnomAD, or at least its definition should be similar.

In other words: find a reasonable estimate for $\text{max. credible AF}\_\text{C}$, and filter FAF95 against that.
> 1 votes

---
Tags: phylogenetics, allele-frequency, variant-filtration
---


# Question
Title: How to find the file names referred to in the R script produced by GEO2R in GEO Accession
I am new to both coding and GEO so this may be a silly question.

I need the raw\_counts data set for the GEO data set I am using but data in the downloads box is not the raw\_counts. How would I be able to find the raw counts?

As I am new to both I decided to use GEO2R to run an analysis and then to see if I could create the same results using code in R as I new that I could look at the R script for GEO2R if needed.

I was struggling with which of the files I needed to download so I looked at the R script used in GEO2R but the files that are being used are completely different to the ones that I have downloaded.

I am using the data set for Series GSE178646 "CRISPR gene editing and inducible pluripotent stem cell neuronal disease modelling for rare disease diagnosis: EMHM1 genetic variant analysis in Kleefstra Syndrome" and am only interested in CF7-12 but the files that I have downloaded are called `"GSM5395222_CF7_kallisto_abundance.tsv.gz"` where as in the R script in GEO2R the code has used the file `"GSE178646_raw_counts_GRCh38.p13_NCBI.tsv.gz"` and I am really not sure where I could access this file.

# Answer
> 2 votes
All the data associated with this project is here:

https://www.ncbi.nlm.nih.gov/gds/?term=GSE178646

The actual download page and project description is here https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE178646

The dataset you seek **GSE178646\_raw\_counts** is here:

https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE178646&format=file

The NPC EHMT1\_SNV CF7 cells for **GSM5395222\_CF7\_kallisto\_abundance.tsv.gz** is on the same page here.

If you now click on the **GEO2 analysis** on this page it will take your here where you select your samples and begin the GEO2 analysis.

---

From the comments,

How to do this?

Its easy simply go to https://www.ncbi.nlm.nih.gov

In the check box search for the option **GEO datasets**, in this example I queried **GSM5395222** in the box. This goes to the search results here, the top hit called "**1**" is document you want here. The very top hit (above the **1**) is actual data in question, minus the project information page.

All the information is within this document and simply use a "search text in webpage" using **GSM5395222**

---
Tags: geo
---


# Question
Title: Bioprocess data visualisation and analysis
I'm looking for some kind of software, ideally FOSS, that will allow us to plot, compare, annotate and add additional points to time series data from bioprocessing runs (bioreactors). Many of the bioreactors come with some of their own quite limited software to do this, but it cannot be installed on other computers.

I have looked around a fair amount but can only expensive-looking biopharma software. Is there any academic software to achieve this? I was considering that writing it myself wouldn't be that difficult with python, but I don't want to waste effort.

(maybe this SE is not the best for bioprocessing, not sure where else to ask.)

# Answer
> 0 votes
As a general answer the technique here is called "Internet of Things" (IoT), where a machine, or an "IoT" device reports to either the cloud or another computer. In bioinformatics this technology has not embraced at present outside the industrial sector.

Google Cloud Platform (GCP) does have an IoT platform where the devices can be registered and data is streamed to the cloud. However, doing this operation across a LAN (local area network) is automatic and can be directed via `cron` notably in automating and collecting the data and analysing it periodically. For a generic IoT stack consider https://www.mdpi.com/2076-3417/11/24/11932

---
Tags: software-recommendation
---


# Question
Title: phytools add.species.to.genus could not match your species to a genus
I have a species tree to which I am trying to add species, which includes a tip "Geissospermum laeve". When I run: `with_polytomes = phytools::add.species.to.genus(tree, "Geissospermum polytome1", where='root')`, the original tree is returned with the error message: 'could not match your species to a genus check spelling, including case'.

Reproducible example:

```
test_tree = ape::rtree(5)
test_tree$tip.label = c("t3 sp1", "t4 sp1", "t5 sp1", "t2 sp1", "t1 sp1")
with_polytomes = phytools::add.species.to.genus(test_tree, "t1 polytome1", where='root')

```

I don't understand why this is happening, if anyone could shed some light on what's causing this error that would be great

# Answer
Two things to try,

1. setting `genus=NULL`, The logic being it doesn't know which genus to add it to. Just as a precaution setting all blank space within a taxa id in the tree to `_` ...

```
add.species.to.genus(test_tree, "t1_polytome1", genus=NULL, where=c("root"))

```

Maybe replace `root` with `random` to get it working.

2. alternative using `tree` instead

```
genus.to.species.tree(test_tree, "t1_polytome1")

```
> 1 votes

---
Tags: phylogenetics, phylogeny
---


# Question
Title: How do I keep a glob from copying a file in Cromwell / WDL
Using Crowell 52. When I do a glob to get matching files, Cromwell places a copy of each matching file in the glob-\* directory.

As these are multi-gigabyte Bam files, I'd really prefer for Cromwell to store soft links to the files, rather than full copies.

The code:

```
   output {
      Array[File] Outputs = glob("*.fastq.gz")
   }

```

Is it possible to make this happen?

# Answer
As stated in https://github.com/broadinstitute/cromwell/issues/4395

The problem may be because you can't hardlink folders on your system for various reasons (i.e. filesystem doesn't allow hardlinking folders, only files, or you are trying to link between different filesystems).

You can try overriding the default way glob works.

i.e. add this to your cromwell.config

```
glob-link-command = "ls -L GLOB_PATTERN 2> /dev/null | xargs -I ? ln -s ? GLOB_DIRECTORY"

```

Like this

```
backend
{
  default = "Slurm"
  providers
  {
    Htcondor {
      actor-factory = "cromwell.backend.impl.sfs.config.ConfigBackendLifecycleActorFactory"

      config
      {
        glob-link-command = "find GLOB_PATTERN 2> /dev/null | xargs -I ? ln -s ? GLOB_DIRECTORY"

```
> 1 votes

---
Tags: wdl
---


# Question
Title: Help with finding an unsupervised learning algorithm for RNA sequence/structure pattern identification
I have a bunch of RNA sequences (and their optimal secondary structures) and their corresponding energy values (measured by mean free energy) and I'm trying to find a way to identify features (patterns in their structures or sequences) common between samples of similar energy values. the overall goal is to try and identify sequence or structural elements in the RNA samples that tend to correlate with low mean free energy values. would be grateful if anyone could recommend me algorithms to look into using - I'm assuming this would be an unsupervised learning project and I only have experience with supervised stuff, but I'm looking into PCA right now and not sure if that'd be useful. should I be looking into something else?

Intuitively I'm imagining it as like a clustering + feature extraction problem where I have a bunch of dots, each representing an RNA sample, and then the axis represents energy so the dots could be clustered in energy value similarity and then within each energy value cluster patterns/relationships could be found between the sequences and structures of the samples. but not sure if an ML algo exists to handle this and PCA seems like not what I'm looking for because the axes would be principal components and not energy...

# Answer
It depends what is meant by "features". Features can be a particular structure component in the RNA structure. In ML/DL these "features" are called instances (or observations). "Features" is also a technical word used in machine/deep learning, which can be completely different, it means the data or statistical variables associated with an instance.

Feature weight extraction (ML/DL definition) using unsupervised learning is complicated, but doable. This is where the features corresponding to each entry are weighted, i.e. quantified. Under traditional unsupervised learning this is not achievable, unless you want to include/exclude "features" and assess cluster integrity. I personally considered it a lot of work for no real interpretable results.

If the features are components of the RNA structure (what would be technically called 'instances' in ML/DL) then the general traditional unsupervised learning following PCA analysis is tSNE. However, tSNE is often used to get better resolution the PCA feature reduction output, i.e. tightens the clusters. PCA-\>tSNE is a cool technique and works well. PCA and tSNE can give different clusters by the way. As a rule of thumb I consider PCA \> tSNE.

Multi-dimensional scaling etc ... needs a defined analytics aim.

---

Information on tSNE and PCA this online resource looks fine

https://www.datacamp.com/tutorial/introduction-t-sne

The above is if you want to do tSNE and compare the results to PCA this link is good.

In terms of using tSNE to sharpen PCA clusters I always used O'Reilly's "Introduction to Machine Learning with Python". There is a very clear example of Python code of the technique there, towards the end of the book and why it works.
> 0 votes

---
Tags: rna, clustering, features, unsupervised-learning
---


# Question
Title: Multiple SRR ID downloaded for single GSM IDs input issue
So here I'm trying to donwload multiple fastq files based on the GSM ID input the script works fine when then input is only SRR id but it runs into error when it is given GSM ID

My small GSM ID input set

```
cat srr_id.txt
GSM5369701
GSM4663476

```

# Answer
> 1 votes
Basically you're asking why a given GEO accession (GSM...) has many Runs (SRR...). That is because the submitters uploaded multiple fastq files for the same sample, which is the case if a library has been sequenced over multiple lanes or flowcells, often to increase depth. You can use `cat` to combine these "technical replicates" into a single fastq file for both R1 and R2 file in case of paired-end data, or one single fastq file in case of single-end data.

For example, if you have `file1.fq.gz`, `file2.fq.gz` and `file3.fq.gz` just do:

```
cat file1.fq.gz file2.fq.gz file3.fq.gz > combined.fq.gz

```

---
Tags: fastq, shell, sratoolkit
---


# Question
Title: Applying PCA or tSNE on non-numerical data
Reading through this t-sne blogpost and I'm a little unclear on how i'd use either of these algorithms here for my goals. i get it'd cluster my samples based on similarity in terms of first and second principal components/t-SNEs, but what could i then use identify structural or sequence patterns within each cluster? please let me know if I'm missing something - thanks!

I'm also curious as to how I would apply PCA or t-SNE on my data considering its non-numerical nature. I was thinking of applying one-hot encoding, but is there a generally more suggested approach? I've attached a screenshot of a few of my dataset rows for context (for the sake of having the table fit, I chose rows with simpler structures/zero MFE values, but most rows have complex structures and non-zero MFE values). Thanks.

---

This is based on my previous post:

> I have a bunch of RNA sequences (and their optimal secondary structures) and their corresponding energy values (measured by mean free energy) and I'm trying to find a way to identify features (patterns in their structures or sequences) common between samples of similar energy values. the overall goal is to try and identify sequence or structural elements in the RNA samples that tend to correlate with low mean free energy values. would be grateful if anyone could recommend me algorithms to look into using - I'm assuming this would be an unsupervised learning project and I only have experience with supervised stuff, but I'm looking into PCA right now and not sure if that'd be useful. should I be looking into something else?

> Intuitively I'm imagining it as like a clustering + feature extraction problem where I have a bunch of dots, each representing an RNA sample, and then the axis represents energy so the dots could be clustered in energy value similarity and then within each energy value cluster patterns/relationships could be found between the sequences and structures of the samples. but not sure if an ML algo exists to handle this and PCA seems like not what I'm looking for because the axes would be principal components and not energy...

# Answer
I'd recommend PCA, tSNE the PCA (first)+ tSNE (second).

The sequence data - one hot key encoding is a good idea.

```
import pandas as pd
datadict = {'UTR':[.....],'Sequence':['AGAA','ACAG' ...], 'Optimal Secondary Structure':[....]}
df = pd.DataFrame(datadict)
df_onehot = pd.get_dummies(df['Sequence'])
df = pd.concat([df, df_onehot], axis=1)     # add the one hot encoded column to master dataframe
df.drop(...) # remove the original 'Sequence' column

```

---

The other issue is vectorising the "**dots**". This is really a separate question, here's what I'd do ...

For optimal secondary structure, maybe simply count the dots?

```
count = df_series.str.count(".")

```

This only works on a series so you'll need to do something like `df_list = df['Dots column'].to_list()` and then `df_series = pd.Series(df_list)`. Finally the `count` will need to be `concat` back onto the `df` and the original `Dots` deleted.

The thing I'm not clear about is your 'instances' that could be 'UTR\_info' thats the stuff you are calling you data or "features".
> 1 votes

---
Tags: rna, clustering, unsupervised-learning
---


# Question
Title: How to get coefficient tables from multiple regression model result!
I've done multiple regression modeling in R using the below script. Now I'm wondering how can I make a table for the coefficients of each model of all models from lm, glm, and rf separately.

```

all_variables <- names(PaC09c_w20_dat3)
response_variables <- all_variables[c(1:3)]
predictors <- all_variables[-c(1:3)]

lm_model <- lapply(
  response_variables,
  function(x) lm(reformulate(termlabels = predictors, response = x), data = PaC09c_w20_dat3)
) |>
  setNames(response_variables)

glm_model <- lapply(
  response_variables,
  function(x) glm(reformulate(termlabels = predictors, response = x), data = PaC09c_w20_dat3)
) |>
  setNames(response_variables)

rf_ranger <- lapply(
  response_variables,
  function(x) ranger(reformulate(termlabels = predictors, response = x), data = PaC09c_w20_dat3, importance="impurity")
) |>
  setNames(response_variables)

save.image(file="PaC09c_w20_Reg.RData")

```

Input data looks like (Original Table is very large),

```
  head(PaC09c_w20_dat3)
        EE87865ln1   EE87866ln1   EE87895ln1  blood_vessel_w20 adrenal_gland_w20 bone_element_w20 bronchus_w20
               <dbl>       <dbl>      <dbl>     <dbl>             <dbl>            <dbl>        <dbl>
        1    0.00391     0.00326    0.00332       0                 0                1            0
        2    0.00139     0.00116    0.00132       0                 0                0            0
        3    0.00360     0.00270    0.00469       1                 1                0            1
        4    0.00323     0.00348    0.00339       0                 0                1            0
        5    0.00323     0.00330    0.00382       0                 1                0            0
        6    0.00278     0.00208    0.00214       0                 0                1            0

```

The expected output would look like this,

```
                           EE87865_lmcof   EE87866_lmcof
blood_vessel_w20           6.332633e-06    5.916615e-06
adrenal_gland_w20         -2.159470e-06   -2.399176e-06
bone_element_w20           8.578804e-06    8.106153e-06
bronchus_w20              -2.956291e-05   -2.152231e-05

```

# Answer
> 1 votes
Well, I found a solution,

```
 coefs_lm<-as.data.frame(lapply(lm_model,function(x)coef(x)[2:5]))

 coefs_glm<-as.data.frame(lapply(glm_model,function(x)coef(x)[2:5]))

var_imp <- function(object, ...) {
  var_imp_r <- object$variable.importance
  var_imp_r
}

coefs_rf<-as.data.frame(lapply(rf_ranger,function(x)var_imp(x)[1:4]))

```

`2:5` for `lm` and `glm` is the number of samples `EE..`, where 1 is intercept.

# Answer
> 1 votes
To get the coefficients alone

```
lm_model$coefficients
glm_model$coefficients

```

Everything including coefficients its,

```
summary(lm_model)
summary(glm_model)

```

I assume this is the same for `summary(rf_ranger)`, albeit not sure whether `coefficients` would work here

---
Tags: r, linear-regression
---


title: Using bioconductor from Python

text:
# Question
Title: Using bioconductor from Python
Has anyone used bioconductor from Python? Is there any reason I would choose to use it from R instead of Python? It seems like there is a Python extension for it. Also, is there any reason to use bioconductor over Biopython?

# Answer
> 5 votes
Bioconductor packages are R packages that can only be used directly from python via rpy2 or similar. If you've ever had reason to use rpy2 you'd probably realize quickly that while it works, it's no where near as easy as doing the same thing directly in R.

Biopython is a package for doing a number of useful things, but it can in no way come remotely close to performing even a percent of what Bioconductor packages do. That would be an impossible feat for a single python package.

# Answer
> 1 votes
Recently I have done something similar for my work, I am using sangerseqR which is s Bioconductor package in Python Environment using rpy2

```
    import rpy2.robjects as r
    from rpy2.robjects.packages import importr
    utils = importr('utils')
    utils.install_packages('sangerseqR', repos="https://git.bioconductor.org/packages/sangerseqR")
    utils.chooseBioCmirror(ind=1) # select the first mirror in the list
    utils.install_packages('sangerseqR')

    sangerseqR = importr('sangerseqR')

   #now using the sangerseqR package to the read the sequcence
   Trace = sangerseqR.readsangerseq( file="1I1_F_P1815443_047.scf")

```

# Answer
> 1 votes
I had some trouble installing rpy2 because of binary incompatibility issues. To solve this problem, I took the following steps:

First, I wrote the necessary functions in R. Then, I created a straightforward REST-API using `library(plumber)` and started the server on my local machine. On the Python side, I utilized `requests` to call the functions and receive responses.

I know that this solution is somewhat ugly, but it has works for me

---
Tags: r, python, bioconductor
---


# Question
Title: Importing local files to a GATK docker image
I have pulled GATK docker on my computer but when running a GATK command, can not locate my input file

**My OS version**

```
cat /etc/os-release
NAME="Ubuntu"
VERSION="18.04.6 LTS (Bionic Beaver)"
ID=ubuntu
ID_LIKE=debian
PRETTY_NAME="Ubuntu 18.04.6 LTS"
VERSION_ID="18.04"
HOME_URL="https://www.ubuntu.com/"
SUPPORT_URL="https://help.ubuntu.com/"
BUG_REPORT_URL="https://bugs.launchpad.net/ubuntu/"
PRIVACY_POLICY_URL="https://www.ubuntu.com/legal/terms-and-policies/privacy-policy"
VERSION_CODENAME=bionic
UBUNTU_CODENAME=bionic

```

Files in my working directory

```
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ ls
GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz
NG-27280_CLTSS_LTS_001A_lib506241_7636_2_MarkedDup.bam
NG-27280_CLTSS_LTS_001A_lib506241_7636_2_MarkedDup.bam.bai
NG-27280_CLTSS_LTS_001A_lib506241_7636_2_MarkedDuplicates.txt
NG-27280_CLTSS_LTS_001B_lib506242_7915_2_MarkedDup.bam
NG-27280_CLTSS_LTS_001B_lib506242_7915_2_MarkedDup.bam.bai
NG-27280_CLTSS_LTS_001B_lib506242_7915_2_MarkedDuplicates.txt
NG-27280_CLTSS_LTS_002B_lib488377_7494_1_MarkedDup.bam
NG-27280_CLTSS_LTS_002B_lib488377_7494_1_MarkedDup.bam.bai
NG-27280_CLTSS_LTS_002B_lib488377_7494_1_MarkedDuplicates.txt
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$

```

I then pulled GATK docker

```
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ docker run -v ~/vcf:/gatk/vcf -it broadinstitute/gatk:latest
(gatk) root@8b2af4d7aadf:/gatk#

(gatk) root@8b2af4d7aadf:/gatk# ls
GATKConfig.EXAMPLE.properties  gatk-package-4.4.0.0-local.jar  gatkPythonPackageArchive.zip  scripts
README.md              gatk-package-4.4.0.0-spark.jar  gatkcondaenv.yml          vcf
gatk                   gatk-spark.jar              gatkenv.rc
gatk-completion.sh         gatk.jar                run_unit_tests.sh
(gatk) root@8b2af4d7aadf:/gatk#

```

As you see, folder seems to be empty because I can not connect my local working directory to container

```
(gatk) root@6d3d3ff69e21:/gatk# cd vcf
(gatk) root@6d3d3ff69e21:/gatk/vcf# ls
(gatk) root@6d3d3ff69e21:/gatk/vcf#

```

I have tried all of these but did not work

```
fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ docker run -v /data/Continuum/WES/vc:/gatk/ -it broadinstitute/gatk:latest

fi1d18@RBGO-Server2:/data/Continuum/WES/vcf$ docker run -v ~/data/Continuum/WES/vc:/gatk/ -it broadinstitute/gatk:latest

```

# Answer
@terdon answered that

```
docker run -v /data/Continuum/WES/:/gatk/vcf -it broadinstitute/gatk:latest

```
> 1 votes

---
Tags: gatk, docker, wes
---


# Question
Title: Invalid tag name: "1KG" and Invalid character '/' in 'SA' FORMAT field at chr1:6197766
I have a VCF file that is of the following format:

```
##fileformat=VCFv4.1
##fileDate=20210329-231145
##AnalysisPipelineName=NGS5
##AnalysisPipelineVersion=V5.2.1
##AnalysisConfigurationName=NGS5_Exome
##AnalysisConfigurationVersion=5.2.1
##NGSPanelName=SureSelect_ExomeV7_plus_720G
##NGSPanelVersion=V1.0
##BatchID=WES_032721_A
##Lane=Lane1
##reference=hg38
##PathogenicSilentFile=Pathogenic_Silent.20200730.vcf
##FlowcellName=210327_A01233_0050_AHYFJHDSXY
##VariantDefDateTime=20210323-070701
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=VF,Number=1,Type=Float,Description="Variant Frequency">
##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allele Depth">
##FORMAT=<ID=SB,Number=1,Type=String,Description="Strand Bias with read numbers (Fw/Rv)">
##FORMAT=<ID=SA,Number=4,Type=Integer,Description="Number of 1) forward ref alleles; 2) reverse ref; 3) forward non-ref; 4) reverse non-ref alleles, used in variant calling">
##FORMAT=<ID=SP,Number=1,Type=Float,Description="Phred of Chi-square p-value for strand bias">
##FORMAT=<ID=SO,Number=1,Type=String,Description="abs(log10( (RefFw*AltRv)/(RefRv*AltFw) )) for strand bias">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=TI,Number=.,Type=String,Description="Transcript ID">
##INFO=<ID=GI,Number=.,Type=String,Description="Gene ID">
##INFO=<ID=LO,Number=.,Type=String,Description="Genetic location">
##INFO=<ID=FC,Number=.,Type=String,Description="Functional Consequence">
##INFO=<ID=PC,Number=.,Type=String,Description="Protein level change">
##INFO=<ID=DC,Number=.,Type=String,Description="DNA level change">
##INFO=<ID=EXON,Number=.,Type=String,Description="Exon number">
##INFO=<ID=CI,Number=.,Type=String,Description="Clinical impact">
##INFO=<ID=dbSNP,Number=.,Type=String,Description="dbSNP version">
##INFO=<ID=1KG,Number=0,Type=String,Description="Existing in the 1000 Genome phase 3 data">
##FILTER=<ID=R8,Description="IndelRepeatLength is greater than 8">
##FILTER=<ID=R8.1,Description="IndelRepeatLength of a monomer is greater than 8">
##FILTER=<ID=R8.2,Description="IndelRepeatLength of a dimer is greater than 8">
##FILTER=<ID=sb,Description="Variant strand bias high">
##FILTER=<ID=sb.s,Description="Variant strand bias significantly high (only for SNV)">
##FILTER=<ID=rs,Description="Variant with rs (dbSNP) number in a non-core gene">
##FILTER=<ID=FP,Description="Possibly false positives due to high similarity to off-target regions">
##FILTER=<ID=NC,Description="Noncoding INDELs on non-core genes">
##FILTER=<ID=lowDP,Description="low depth variant">
##FILTER=<ID=Benign,Description="Benign variant">
##contig=<ID=chrM,length=16569>
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=242193529>
##contig=<ID=chr3,length=198295559>
##contig=<ID=chr4,length=190214555>
##contig=<ID=chr5,length=181538259>
##contig=<ID=chr6,length=170805979>
##contig=<ID=chr7,length=159345973>
##contig=<ID=chr8,length=145138636>
##contig=<ID=chr9,length=138394717>
##contig=<ID=chr10,length=133797422>
##contig=<ID=chr11,length=135086622>
##contig=<ID=chr12,length=133275309>
##contig=<ID=chr13,length=114364328>
##contig=<ID=chr14,length=107043718>
##contig=<ID=chr15,length=101991189>
##contig=<ID=chr16,length=90338345>
##contig=<ID=chr17,length=83257441>
##contig=<ID=chr18,length=80373285>
##contig=<ID=chr19,length=58617616>
##contig=<ID=chr20,length=64444167>
##contig=<ID=chr21,length=46709983>
##contig=<ID=chr22,length=50818468>
##contig=<ID=chrX,length=156040895>
##contig=<ID=chrY,length=57227415>
##contig=<ID=chrUn_GL000195v1,length=182896>
##contig=<ID=chrUn_GL000213v1,length=164239>
##contig=<ID=chrUn_GL000214v1,length=137718>
##contig=<ID=chrUn_GL000216v2,length=176608>
##contig=<ID=chrUn_GL000218v1,length=161147>
##contig=<ID=chrUn_GL000219v1,length=179198>
##contig=<ID=chrUn_GL000220v1,length=161802>
##contig=<ID=chrUn_GL000224v1,length=179693>
##contig=<ID=chrUn_GL000226v1,length=15008>
##contig=<ID=chrUn_KI270302v1,length=2274>
##contig=<ID=chrUn_KI270303v1,length=1942>
##contig=<ID=chrUn_KI270304v1,length=2165>
##contig=<ID=chrUn_KI270305v1,length=1472>
##contig=<ID=chrUn_KI270310v1,length=1201>
##contig=<ID=chrUn_KI270311v1,length=12399>
##contig=<ID=chrUn_KI270312v1,length=998>
##contig=<ID=chrUn_KI270315v1,length=2276>
##contig=<ID=chrUn_KI270316v1,length=1444>
##contig=<ID=chrUn_KI270317v1,length=37690>
##contig=<ID=chrUn_KI270320v1,length=4416>
##contig=<ID=chrUn_KI270322v1,length=21476>
##contig=<ID=chrUn_KI270329v1,length=1040>
##contig=<ID=chrUn_KI270330v1,length=1652>
##contig=<ID=chrUn_KI270333v1,length=2699>
##contig=<ID=chrUn_KI270334v1,length=1368>
##contig=<ID=chrUn_KI270335v1,length=1048>
##contig=<ID=chrUn_KI270336v1,length=1026>
##contig=<ID=chrUn_KI270337v1,length=1121>
##contig=<ID=chrUn_KI270338v1,length=1428>
##contig=<ID=chrUn_KI270340v1,length=1428>
##contig=<ID=chrUn_KI270362v1,length=3530>
##contig=<ID=chrUn_KI270363v1,length=1803>
##contig=<ID=chrUn_KI270364v1,length=2855>
##contig=<ID=chrUn_KI270366v1,length=8320>
##contig=<ID=chrUn_KI270371v1,length=2805>
##contig=<ID=chrUn_KI270372v1,length=1650>
##contig=<ID=chrUn_KI270373v1,length=1451>
##contig=<ID=chrUn_KI270374v1,length=2656>
##contig=<ID=chrUn_KI270375v1,length=2378>
##contig=<ID=chrUn_KI270376v1,length=1136>
##contig=<ID=chrUn_KI270378v1,length=1048>
##contig=<ID=chrUn_KI270379v1,length=1045>
##contig=<ID=chrUn_KI270381v1,length=1930>
##contig=<ID=chrUn_KI270382v1,length=4215>
##contig=<ID=chrUn_KI270383v1,length=1750>
##contig=<ID=chrUn_KI270384v1,length=1658>
##contig=<ID=chrUn_KI270385v1,length=990>
##contig=<ID=chrUn_KI270386v1,length=1788>
##contig=<ID=chrUn_KI270387v1,length=1537>
##contig=<ID=chrUn_KI270388v1,length=1216>
##contig=<ID=chrUn_KI270389v1,length=1298>
##contig=<ID=chrUn_KI270390v1,length=2387>
##contig=<ID=chrUn_KI270391v1,length=1484>
##contig=<ID=chrUn_KI270392v1,length=971>
##contig=<ID=chrUn_KI270393v1,length=1308>
##contig=<ID=chrUn_KI270394v1,length=970>
##contig=<ID=chrUn_KI270395v1,length=1143>
##contig=<ID=chrUn_KI270396v1,length=1880>
##contig=<ID=chrUn_KI270411v1,length=2646>
##contig=<ID=chrUn_KI270412v1,length=1179>
##contig=<ID=chrUn_KI270414v1,length=2489>
##contig=<ID=chrUn_KI270417v1,length=2043>
##contig=<ID=chrUn_KI270418v1,length=2145>
##contig=<ID=chrUn_KI270419v1,length=1029>
##contig=<ID=chrUn_KI270420v1,length=2321>
##contig=<ID=chrUn_KI270422v1,length=1445>
##contig=<ID=chrUn_KI270423v1,length=981>
##contig=<ID=chrUn_KI270424v1,length=2140>
##contig=<ID=chrUn_KI270425v1,length=1884>
##contig=<ID=chrUn_KI270429v1,length=1361>
##contig=<ID=chrUn_KI270435v1,length=92983>
##contig=<ID=chrUn_KI270438v1,length=112505>
##contig=<ID=chrUn_KI270442v1,length=392061>
##contig=<ID=chrUn_KI270448v1,length=7992>
##contig=<ID=chrUn_KI270465v1,length=1774>
##contig=<ID=chrUn_KI270466v1,length=1233>
##contig=<ID=chrUn_KI270467v1,length=3920>
##contig=<ID=chrUn_KI270468v1,length=4055>
##contig=<ID=chrUn_KI270507v1,length=5353>
##contig=<ID=chrUn_KI270508v1,length=1951>
##contig=<ID=chrUn_KI270509v1,length=2318>
##contig=<ID=chrUn_KI270510v1,length=2415>
##contig=<ID=chrUn_KI270511v1,length=8127>
##contig=<ID=chrUn_KI270512v1,length=22689>
##contig=<ID=chrUn_KI270515v1,length=6361>
##contig=<ID=chrUn_KI270516v1,length=1300>
##contig=<ID=chrUn_KI270517v1,length=3253>
##contig=<ID=chrUn_KI270518v1,length=2186>
##contig=<ID=chrUn_KI270519v1,length=138126>
##contig=<ID=chrUn_KI270521v1,length=7642>
##contig=<ID=chrUn_KI270522v1,length=5674>
##contig=<ID=chrUn_KI270528v1,length=2983>
##contig=<ID=chrUn_KI270529v1,length=1899>
##contig=<ID=chrUn_KI270530v1,length=2168>
##contig=<ID=chrUn_KI270538v1,length=91309>
##contig=<ID=chrUn_KI270539v1,length=993>
##contig=<ID=chrUn_KI270544v1,length=1202>
##contig=<ID=chrUn_KI270548v1,length=1599>
##contig=<ID=chrUn_KI270579v1,length=31033>
##contig=<ID=chrUn_KI270580v1,length=1553>
##contig=<ID=chrUn_KI270581v1,length=7046>
##contig=<ID=chrUn_KI270582v1,length=6504>
##contig=<ID=chrUn_KI270583v1,length=1400>
##contig=<ID=chrUn_KI270584v1,length=4513>
##contig=<ID=chrUn_KI270587v1,length=2969>
##contig=<ID=chrUn_KI270588v1,length=6158>
##contig=<ID=chrUn_KI270589v1,length=44474>
##contig=<ID=chrUn_KI270590v1,length=4685>
##contig=<ID=chrUn_KI270591v1,length=5796>
##contig=<ID=chrUn_KI270593v1,length=3041>
##contig=<ID=chrUn_KI270741v1,length=157432>
##contig=<ID=chrUn_KI270742v1,length=186739>
##contig=<ID=chrUn_KI270743v1,length=210658>
##contig=<ID=chrUn_KI270744v1,length=168472>
##contig=<ID=chrUn_KI270745v1,length=41891>
##contig=<ID=chrUn_KI270746v1,length=66486>
##contig=<ID=chrUn_KI270747v1,length=198735>
##contig=<ID=chrUn_KI270748v1,length=93321>
##contig=<ID=chrUn_KI270749v1,length=158759>
##contig=<ID=chrUn_KI270750v1,length=148850>
##contig=<ID=chrUn_KI270751v1,length=150742>
##contig=<ID=chrUn_KI270752v1,length=27745>
##contig=<ID=chrUn_KI270753v1,length=62944>
##contig=<ID=chrUn_KI270754v1,length=40191>
##contig=<ID=chrUn_KI270755v1,length=36723>
##contig=<ID=chrUn_KI270756v1,length=79590>
##contig=<ID=chrUn_KI270757v1,length=71251>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  TN21-116928
chr1    6197766 rs2294714   T   C   0.0 rs  DP=1931;TI=NM_000983.3;GI=RPL22;FC=Silent;DC=c.13-10A>G;LO=INTRON;EXON=2    GT:VF:AD:SB:SA:SP:SO    0/1:0.483:999,932:883/1048:451/548,432/500:2:0.0
chr1    7677739 rs3011926   A   G   0.0 rs  DP=851;TI=NM_015215.3;GI=CAMTA1;FC=Silent;DC=c.2914+6A>G;LO=INTRON;EXON=11  GT:VF:AD:SB:SA:SP:SO    1/1:0.999:1,850:309/542:0/1,309/541:0:INF
chr1    9724365 rs11121484  C   T   0.0 rs  DP=1664;TI=NM_005026.4;GI=PIK3CD;FC=Synonymous;PC=Y936=;DC=c.2808C>T;LO=EXON;EXON=22    GT:VF:AD:SB:SA:SP:SO    0/1:0.455:907,757:842/822:448/459,394/363:6:0.0
chr1    10258594    rs12402052  C   G   0.0 rs  DP=1400;TI=NM_015074.3;GI=KIF1B;FC=Synonymous;PC=A95=;DC=c.285C>G;LO=EXON;EXON=4    GT:VF:AD:SB:SA:SP:SO    0/1:0.462:751,647:678/722:375/376,302/345:6:0.1
chr1    10295776    rs3753037   C   T   0.0 rs  DP=1050;TI=NM_015074.3;GI=KIF1B;FC=Silent;DC=c.1639+10C>T;LO=INTRON;EXON=17 GT:VF:AD:SB:SA:SP:SO    0/1:0.477:549,501:395/655:211/338,184/317:2:0.0
chr1    10375266    rs11121552  C   A   0.0 rs  DP=1487;TI=NM_015074.3;GI=KIF1B;FC=Synonymous;PC=T1721=;DC=c.5163C>A;LO=EXON;EXON=46    GT:VF:AD:SB:SA:SP:SO    0/1:0.484:767,719:870/617:443/324,427/292:3:0.0
chr1    11121270    rs11121691  C   T   0.0 .   DP=1092;TI=NM_004958.3;GI=MTOR;FC=Synonymous;PC=L2303=;DC=c.6909G>A;LO=EXON;EXON=49 GT:VF:AD:SB:SA:SP:SO    0/1:0.483:565,527:664/428:346/219,318/209:1:0.0
chr1    11130589    rs2275527   G   A   0.0 .   DP=1048;TI=NM_004958.3;GI=MTOR;FC=Synonymous;PC=S1851=;DC=c.5553C>T;LO=EXON;EXON=39 GT:VF:AD:SB:SA:SP:SO    0/1:0.509:515,533:572/476:275/240,297/236:3:0.0
chr1    11145001    rs1057079   C   T   0.0 .   DP=1392;TI=NM_004958.3;GI=MTOR;FC=Synonymous;PC=A1577=;DC=c.4731G>A;LO=EXON;EXON=33 GT:VF:AD:SB:SA:SP:SO    0/1:0.511:681,711:703/689:352/329,351/360:4:0.0
chr1    11228701    rs1064261   G   A   0.0 .   DP=1859;TI=NM_004958.3;GI=MTOR;FC=Synonymous;PC=N999=;DC=c.2997C>T;LO=EXON;EXON=19  GT:VF:AD:SB:SA:SP:SO    0/1:0.490:948,911:1059/800:536/412,523/388:2:0.0
chr1    11241657    rs1135172   A   G   0.0 .   DP=1643;TI=NM_004958.3;GI=MTOR;FC=Synonymous;PC=D479=;DC=c.1437T>C;LO=EXON;EXON=10  GT:VF:AD:SB:SA:SP:SO    0/1:0.512:801,842:793/850:376/425,417/425:5:0.0

```

Unfortunately, when I try to utilize the following Python code to validate my VCF file:

```
import pysam

def validate_vcf(vcf_file_path):
    try:
        # Open the VCF file
        vcf_file = pysam.VariantFile(vcf_file_path)

        # Iterate through the records to check if they can be read
        for record in vcf_file:
            pass

        print("VCF file is valid.")
    except Exception as e:
        print(f"Error validating VCF file: {str(e)}")

# Path to the VCF file you want to validate
vcf_file_path = "path/to/yourfile.vcf"

validate_vcf(vcf_file_path)

```

I receive the following error:

```
[W::bcf_hrec_check] Invalid tag name: "1KG"
[E::vcf_parse_format] Invalid character '/' in 'SA' FORMAT field at chr1:6197766

```

Does anyone know what I can do to fix this?

# Answer
1. Given issues I've had in various programming languages that use numbers at the start of a variable, the tag issue might be related to it starting with a numeric value. Is it possible to change that from `1KG` to something like `ph3_1KG`?
2. According to the information in the VCF header, the SA field should have 4 values:

```
##FORMAT=<ID=SA,Number=4,Type=Integer,Description="Number of 1) forward ref alleles; 2) reverse ref; 3) forward non-ref; 4) reverse non-ref alleles, used in variant calling">

```

While that location does have four numerical values in the actual data for that variant location, it seems like the validation code isn't expecting the forward and reverse alleles to be separated by `/`:

```
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  TN21-116928
chr1    6197766 rs2294714   T   C   0.0 rs  DP=1931;TI=NM_000983.3;GI=RPL22;FC=Silent;DC=c.13-10A>G;LO=INTRON;EXON=2    GT:VF:AD:SB:SA:SP:SO    0/1:0.483:999,932:883/1048:451/548,432/500:2:0.0
                                                                                ^^                                     ^^^^^^^^^^^^^^^

```

Whether or not these errors are actually important will depend on the program that this VCF file is used for.
> 1 votes

---
Tags: python, ngs, vcf
---


# Question
Title: What are the fastest way of feature selection in R for large dataset in Machine Learning?
I have about 100 data sets with over 100,000 rows and 27 predictor variables and multiple response variable on each dataset. I want to select the most essential Features for these datasets. I used Boruta() from the package Boruta, but it takes an extremely long time (about 3 hours for a single sample). Could someone kindly tell me the quickest technique to select key features from a large dataset?

This is the extension from my previous questions:

1. https://stackoverflow.com/questions/76824529/how-to-do-regression-modelling-in-r-multiple-one-at-a-time-to-multiple-column
2. How to get coefficient tables from multiple regression model result!

# Answer
> 0 votes
Feature selection: the fastest thing moving (i.e. the simplest) is naive Bayes ML. I don't use it but KNN (K-nearest neighbor) I do use, which is also computationally simple, i.e. fast. I think it gives better results, i.e. the best results for a simple algorithm.

R ML has improved decision trees BTW - quite popular, but you'd never use this in Python.

---
Tags: r, modelling, linear-regression, features
---


# Question
Title: How can I run the fdnadist command from the EMBOSS package?
My operating system is Mac OS X.

I would like to compute DNA distances between different sequences. For this, the dnadist package from the Phylip software would be good but hard to automate. EMBOSS seems to have an fdnadist command for this. I have installed emboss version 6.6.0 from the bioconda package but if I type `fdnadist`, the command is not found. How can I use `dnadist` or `fdnadist` on the command line?

# Answer
To find `emboss` packages type `wossname`, then hit enter twice. You will see all the programs. There is one called `distmat`. Then do `tfm distmat`, to get `the fkn manual` for `distmat`.
> 2 votes

# Answer
I found fdnadist in embassy-phylip package https://anaconda.org/bioconda/embassy-phylip
> 2 votes

---
Tags: sequence-analysis, phylogenetics, software-installation, emboss
---


# Question
Title: What algorithm should I use to find what RNA sequence/structural elements lead to different labelled MFE values?
To put it super simply: I have RNA sequences, their optimal secondary structures, and the MFE values of those structures. I'm trying to computationally identity what structural elements and/or sequence motifs tend to be present in RNAs with lower MFEs.

Was thinking about using unsupervised learning for this, but is this the wrong approach?

Past version of this question

# Answer
> 0 votes
I would maintain traditional unsupervised learning is a very good method here (PCA and/or tSNE). The number 1 question is which column is the observation (label)? That depends on the aim of the study. The absence of a "target" column is the thing that separate unsupervised learning from supervised learning.

The "observation" (label), i.e. the column that provides the name, BUT doesn't take part in the calculation, will either be a gene or an isolate/species/taxonomic-group and the aim is to see which genes or isolates cluster together. The implication is that the MFE (MLE?) is the observation (label) ... can be I guess. Again its about the aim of the study.

The trick with unsupervised learning is to get as many "features" (columns of data) for each "observation" (probably gene name or taxonomic group). MLE (MFE?) energy (entropy?) has a value, but there are other values (which I forget but other entropy values), those can be included. Dumbie variables are used for e.g. sequence data, but if there are counts will work and experimental data measurements can be included. There are a few rules:

1. When you have a cateogorical variable as a column in a dataframe it must be designated as such, pandas will allow this. Similarly if you have continuous variable, MLE that must be designated as a float and counts will be designated as integers. You don't want a float to appear as a categorical variable, that will be disastrous. You need to set the data type of each column (see below).
2. In PCA a reasonable number of components must comprise a large amount (sum) of the variance, e.g. \>70%.
3. You can only see realistically two components at any one time on a graph, just plot them out, e.g. PC1 versus PC2, then PC1 vs PC3. Look at which give nice clusters. Don't just look at PC1 and PC2 and believe thats the only answer (which is what everyone does).

**Point 1** is important if the PCA or tSNE may not permit mixed data types (I can't remember), but if so it will complain.

**Point 2** If the sum of the first three principle components exceed 70% of the variance - thats very, very good. If its alot of then you might consider performing tSNE on the PCs to tighten the result.

Thats about it, what you are looking for is clusters that give some indication about the way the data is behaving. This can be the end point of the analysis or it can be the hypothesis for further testing.

Very clear distinct cluster which make biological sense is a good result with unsupervised learning. Obviously the issue is that the clusters are not significant, but if they look good and make sense that's the key to the analysis.

---
Tags: rna
---