# Question
Title: Find protein in DNA sequences with arbitrary encoding and possible frameshifts
I have a dataset with lots of DNA sequences (~10 M seqs, each ~6 KBp long). I want to detect presence of a fixed given protein (~600 Bp / 200 aa) within each sequence and, if detected, obtain the frame and position of the protein within the sequence. However, the protein may be coded arbitrarily and differently in each sequence. In addition, I need to account for frameshifts induced by sequencing errors.

I've tried to use `tfasty36` from the FASTA36 package. This is basically a `tblastn` equivalent (protein query search against DNA database in all frames). Due to the statistical filtering built into these tools, I've come to the conclusion that `tfasty36` and `tblastn` are not the right tools for this job.

`tfasty36` didn't work because it requires fiddling with statistical tolerances to show all matches. Also, the matches are quite fragmented (I guess, it uses heuristics instead of dynamic programming). For `tblastn`, building the database seems to be unnecessary (given that it may have the same problems as `tfasty36` afterwards).

I've looked into alignment tools, but couldn't find any that can align a protein sequence to a DNA sequence while taking frameshifts into account.

How can I solve this problem?

# Answer
If you need to model introns and splicing, then you should use a tool that is designed for that. One of the limitations of blast and its cousins is that they only look at sequence similarity and do not attempt to join HSPs into possible transcript models.

I haven't worked on this in a long time, but back when I did, the programs I used were:

* exonerate

> Exonerate is a generic tool for pairwise sequence comparison. It allows you to align sequences using a many alignment models, either exhaustive dynamic programming or a variety of heuristics.

In this case, you would use your protein sequence as the query and the set of DNA sequences as the target and try to align using the protein2genome model. Something like this:

```
  exonerate --model protein2genome protein.pep dna.fa > out

```

This will result in (many) alignments of your query protein agains the various genomic sequences. Here's an example I had handy:

```
$ exonerate --model protein2genome -n 1 frog.histone.pep c.elegans.chrI.fa
Command line: [exonerate --model protein2genome -n 1 frog.histone.pep c.elegans.chrI.fa]
Hostname: [oregano]

C4 Alignment:
------------
         Query: frog
        Target: CM041006.1 Caenorhabditis elegans strain DL226 chromosome I, whole genome shotgun sequence:[revcomp]
         Model: protein2genome:local
     Raw score: 732
   Query range: 32 -> 300
  Target range: 3039620 -> 3034073

      33 : LeuGlnLysLysLeuGluGluLeuGluLeuAspGluGlnGlnArgLysArgLeuGl :      51
           |||:!!:!!:!!|||:!!:!!|||! !|||   |||! !|||! !:!!||||||..
           LeuGluGluGlnLeuLysLysLeuGlyLeuThrGluProGlnThrGlnArgLeuSe
 3039620 : CTAGAAGAACAGCTGAAAAAGCTCGGGCTCACGGAGCCACAAACTCAACGGCTCAG : 3039566

      52 : uAlaPheLeuThrGlnLysGlnLysValGlyGluLeuLysAspAspAspPheGluL :      70
           .! !||||||  !   |||:!!   :!!  !||||||!..!!:|||   .!.:!!!
           rGluPheLeuGlnValLysGluGlyIleLysGluLeuSerGluAspMetLeuGlnT
 3039565 : CGAATTCTTGCAGGTCAAGGAGGGCATCAAGGAGTTGAGTGAGGATATGCTACAGA : 3039509

      71 : ysValSerGluLeuGlyAlaGlyAsnGlyGlyValValPheLys  >>>> Targe :      85
             ! !.!!|||||||||  !||||||||||||||||||  !|||
           hrGluGlyGluLeuGlyHisGlyAsnGlyGlyValValAsnLys++
 3039508 : CTGAAGGAGAACTTGGCCATGGAAATGGAGGAGTTGTTAATAAGgt.......... : 3039462

      86 : t Intron 1 >>>>  ValSerHisLysProThrSerLeuIleMetAlaArgLys :      97
           55 bp              !  !|||::!  !|||.!!:!!|||||||||||||||
                          ++CysValHisArgLysThrGlyValIleMetAlaArgLys
 3039461 : ...............agTGTGTTCATCGAAAAACCGGTGTAATTATGGCTCGGAAG : 3039373

      98 : LeuIleHisLeuGluIleLysProAlaIleArgAsnGlnIleIleArgGluLeuGl :     116
           |||:!!||||||||||||||||||:!!:!!|||.!.||||||:!!!:!||||||
           LeuValHisLeuGluIleLysProSerValArgGlnGlnIleValLysGluLeuAl
 3039372 : CTTGTACACCTCGAGATCAAGCCGTCAGTCCGACAACAGATCGTCAAAGAGCTCGC : 3039316

     117 : nValLeuHisGluCysAsnSerProTyrIleValGlyPheTyrGlyAlaPheTyrS :     135
            |||||||||:!!||||||||||||!:!||||||||||||||||||||||||   .
           aValLeuHisLysCysAsnSerProPheIleValGlyPheTyrGlyAlaPheValA
 3039315 : CGTACTTCACAAGTGCAATAGTCCATTCATTGTTGGATTCTACGGGGCTTTTGTGG : 3039259

     136 : erAspGlyGluIleSerIleCysMetGluHisMetAspGlyGlySerLeuAspGln :     153
           .!:!!..!!!:||||||||||||||||||:!!|||||||||  !|||||||||  !
           spAsnAsnAspIleSerIleCysMetGluTyrMetAspGlyLeuSerLeuAspIle
 3039258 : ACAATAACGATATATCAATTTGTATGGAGTATATGGATGGGCTTTCACTTGATATA : 3039205

     154 : ValLeuLysLys  >>>> Target Intron 2 >>>>  AlaGlyLysIlePro :     162
           ||||||||||||           3837 bp           !.!|||::::!!|||
           ValLeuLysLys++                         ++ValGlyArgLeuPro
 3039204 : GTTTTGAAAAAGgt.........................agGTTGGCCGCCTCCCG : 3035341

     163 : GluLysIleLeuGlyLysValSerIleAlaValIleLysGlyLeuThrTyrLeuAr :     181
           ||||||.!!:!!|||::!:!!|||:!!||||||:!!:::|||||||||||||||!:
           GluLysPheValGlyArgIleSerValAlaValValArgGlyLeuThrTyrLeuLy
 3035340 : GAAAAGTTTGTCGGACGAATCTCGGTTGCAGTTGTCCGTGGATTGACGTACTTGAA : 3035284

     182 : gGluLysHisLysIleMetHisArg{A}  >>>> Target Intron 3 >>>>  :     190
           !!!::!!  !||||||:!!||||||{|}            50 bp
           sAspGluIleLysIleLeuHisArg{A}++                         +
 3035283 : AGATGAAATTAAGATTTTGCATAGA{G}gt.........................a : 3035208

     191 :  {sp}ValLysProSerAsnIleLeuValAsnSerArgGlyGluIle--LysLeuC :     207
            {||}|||||||||||||||!!:||||||||||||!.!|||||||||##   ||||
           +{sp}ValLysProSerAsnMetLeuValAsnSerAsnGlyGluIle##---LeuC
 3035207 : g{AT}GTTAAACCATCGAATATGCTTGTCAATAGTAACGGAGAAATTAA---TTAT : 3035157

     208 : ysAspPheGlyValSerGlyGlnLeuIleAspSerMetAlaAsnSerPheValGly :     225
           ||||||||||||||||||||..!|||||||||||||||||||||||||||||||||
           ysAspPheGlyValSerGlyMetLeuIleAspSerMetAlaAsnSerPheValGly
 3035156 : GCGATTTTGGAGTCTCTGGAATGTTGATTGATTCGATGGCCAACTCATTTGTCGGA : 3035103

     226 : ThrArgSerTyrMetSer  >>>> Target Intron 4 >>>>  ProGluArg :     234
           |||||||||||||||:!!            799 bp           |||||||||
           ThrArgSerTyrMetAla++                         ++ProGluArg
 3035102 : ACGAGGAGTTATATGGCGgt.........................agCCCGAACGA : 3034277

     235 : LeuGlnGlyThrHisTyrSerValGlnSerAspIleTrpSerMetGlyLeuSerLe :     253
           |||  !|||:!!||||||!!!:!!..!|||||||||||||||.!.|||||||||||
           LeuThrGlySerHisTyrThrIleSerSerAspIleTrpSerPheGlyLeuSerLe
 3034276 : CTCACAGGATCCCACTATACAATTTCATCGGATATCTGGTCATTTGGATTATCACT : 3034220

     254 : uValGluMetAlaIleGlyArgTyrProIleProProProAspAlaLysGluLeuG :     272
           |||||||:!:  !|||||||||||||||:!:||| !!|||..!  !  !|||! !!
           uValGluLeuLeuIleGlyArgTyrProValProAlaProSerGlnAlaGluTyrA
 3034219 : TGTCGAACTTCTTATTGGACGATATCCAGTGCCTGCACCGTCTCAAGCAGAATATG : 3034163

     273 : luLeuIlePheGlyCysSerVal<->GluGlyAspProAlaSerSerGluLeuAla :     289
               !!!:|||..!   :!!! !   |||  !!!:! !|||..!|||     !! !
           laThrMetPheAsnValAlaGluAsnGluIleGluLeuAlaAspSerLeuGluGlu
 3034162 : CTACAATGTTTAATGTGGCTGAAAATGAGATTGAACTGGCTGATAGTCTTGAAGAG : 3034109

     290 : ProArgProArgProProGlyArgProIleSer :     300
           |||!.!  !!.!||||||.!!!.!|||   |||
           ProAsnTyrHisProProSerAsnProAlaSer
 3034108 : CCCAACTATCATCCGCCGAGCAATCCGGCGTCA : 3034074

vulgar: frog 32 300 . CM041006.1 3039620 3034073 - 732 M 52 156 5 0 2 I 0 51 3 0 2 M 73 219 5 0 2 I 0 3833 3 0 2 M 32 96 S 0 1 5 0 2 I 0 46 3 0 2 S 1 2 M 14 42 F 0 2 G 1 0 M 26 78 5 0 2 I 0 795 3 0 2 M 48 144 G 0 3 M 21 63
-- completed exonerate analysis

```

In addition, exonerate lets you control the output you want. Check out the `--ryo` (roll your own output) option in `man exonerate`.
* genewise

> GeneWise compares a protein sequence to a genomic DNA sequence, allowing for introns and frameshifting errors.

Unlike exonerate, genewise *only* does protein to DNA alignment, but like exonerate it models introns and splicing. 20 years ago, when I did this every day, my impression was that genewise was slower but slightly more sensitive than the default exonerate options (exonerate can be tweaked to be made more sensitive). So I would usually try exonerate first and fall back to genewise if it failed to find a good match.
> 2 votes

---
Tags: sequence-alignment
---


# Question
Title: How to BLAST a protein against a large set of organisms?
I have a list of organisms with their full scientific names and Taxonomy IDs, total ~1500. I want to BLAST my protein against this 1500 genome or proteome. How can I do that?

Here a sample protein that I want to search the following database: http://togodb.org/db/thermobase Please note that database can be downloades as CSV file. I am very comfortable using web tools or programmatic tools such as Perl, Python, bash, or Entrez Programming Utilities.

```
>QCC44558.1 chitinase [Halobacterium salinarum] TaxID 641522631
MPHDRRSYLRTSSAVIASLLAASTPTSAADTPPEWDPDTVYTDGDKATFDGYVWEAKWWTKGDKPGADEW
GPWTQRRPVDGSPGGPTAAFTVTDTTVEPGTAVTVDASDTTGSVDSYEWAFGDGTTAAGVTASHTYDAAG
EYTIELTVTTGDGTTDTATQQITATTSPGDDEFKVVGYYPSWKGTDDYDFYPADVPFDQVTDVLYAFLDV
QPDGTVVLPDSDVDHESLLASFASLKQDRAADTRLKLSIGGWGLSPGFEDAAADQASRERFATTAVDLMR
TYDFDGIDVDWEHPGPRRGKCECGSAQGPANHVALLETVRDRLNDAEVEDGRTYDLSVANGGSDWNAALI
NHREVAAVVDDIYMMAYDFTGVWHGTAGLNAPIYGTPPDYPPSGDAQQYTLETTLTIWKEQGYWVDWMEW
EDHGDPVDDPGTLVLGMPFYGRGCNVENGIWDTFSLSEWQQGDPKYQNDVIPPGTWNDLRGPDDANTGAF
DYGDLAANYEGADAWTKQRNEQGGVPYLWNDSEGVFISYDDPTSIAAKVELAVAEDLGGVMIWELSQDYD
GTLLETINQHTP

```

# Answer
> 1 votes
My solution is simply.

1. perform a web-based Blastp using the following parameters:

* most dissimilar
* 5000 hits (maximum limit)
* Restrict to Archaea (which is what you samples appear to be)\*

If you hit the limit of 5000 results, I would repeat the search for each class or phylum within Archaea. I can't remember if it was re-classified as a kingdom or was just another phylum. Concatenate all hits into a single fasta file.

2. Reg-ex your extremophiles database (personally I'd use pandas), so that the species name of each bacterium is recovered no other information is present. *Or* this is an easily accessible value within a hash/database (this is a better solution).
3. Use Biopython to loop through your blast fasta file and filter using `re` search against the species from your database. If you get a hit save that sequence.

Ideally, the point 2 and 3 would be a pandas merge operation, but looping is easier.

\*, If there are extremophiles which are not Archaea (I don't know much about them), then take each general taxonomic group and repeat the operation for the search criteria "restrict the blastp to ... *specific group*".

---

**Issues** The issues with command line blastp, not via the web, is that access is greatly restricted during US working hours. If you are searching more than 100 hits and carrying out the search during the working day US time - they will not let it through. 100 hits or less is fine however. For \>100 hits (e.g. 5000) it's got to be done well outside US working hours/ weekends. Web-based blastp doesn't carry this restriction of access during US working hours. Basically, the API route is deprioritised for the web-based approach. I assume this is to stop bots - albeit of course a bot can automate through NCBI's webpage.

**Note** The original solution by @terdon might easily work as a remote blastp because each query is restricted to \<100 - so it will always get through during working hours. What I am saying is if you're using Blast API it needs a carefully defined strategy, else shift to the web for blanket approaches.

Entrez operations from a commandline do work ok during US working day time - albeit it's got to be coded to deal with "random" server refusals otherwise it will break the code (sporadic queries are refused). This is just my recent experience nothing more.

Just to be clear - this is potentially a medium sized project. Once the strategy has been decided then specific code that supports this approach can be given in separate questions.

# Answer
> 1 votes
In the interface of blastp, you can limit the search space to organisms you want to. Yo can do this by introducing name species or Tax IDs, you can add more than one, but I bet you don't wanna introduce the 1,500 one by one.

The only other way to do this, is by using blastp in its command-line version. You would need to build a custom database to be used as the searching space, and then run blast on your computer or server using this database. Bulding it is a pain in the ass, as far as I know. You need first to download all the proteins in all organisms in fasta format and then use this compilation as a database (makeblastdb command). I am not aware of a way to download all proteins from various organism given a list of taxonomic IDs, but if any, it has to be NCBI, Uniprot or similar. Then you can use blastp locally together with your custom database.

But why limiting the search space is such a way? You are biasing your results to what you want to get, instead of finding the stuff that actually is closest to your protein of interest. My suggestion is to use the online version of blastp with one of the various databases available and rather than pre-selecting organisms, you can exclude the ones from which you definitely don't want to get results from, and you can also exclude entire families/genera such as primates (taxid:9443), to name an example.

---
Tags: blast, ncbi, blastp, makeblastdb, e-utilities
---


# Question
Title: Is it possible to, all in memory without writing a file, given a list of AlignedSegments and a bam header, build a bam, index and get coverage?
My current job is rate-limited chiefly by the IO overhead of writing reads to their corresponding bam files (they are split up by both cell barcode and contig). The only reason I need to write these mini-split bams is so that I can then use the pysam coverage function to assess read depth at certain locations. I already have all of these reads in memory as a list of AlignedSegments. Is it necessary to write them to a .bam file in order to be able to index and use the coverage function? I feel like there should be a way to have the AlignedSegments be treated as a bam object in memory, possibly even indexed in memory, and then treated as if they were just streamed from a bamfile with an index in order to be able to use the coverage method. Thanks in advance, if this is possible it would increase my processing speed x10. Or, given a list of AlignedSegments, is there a way to get coverage at specific positions across all those AlignedSegments in my list (that would get me what I need).

# Answer
> 3 votes
Assuming you are on an operating system that supports them, like Linux, you can try using named pipes in a ramdisk or tmpfs. Not the most elegant approach, but I don't know if what you are asking for is possible (it may very well be, I just don't know) and this can be a decent workaround.

First, check that you do have a mounted tmpfs. On my system, I have a few:

```
$ mount | awk '$1=="tmpfs"'
tmpfs on /dev/shm type tmpfs (rw,nosuid,nodev,inode64)
tmpfs on /tmp type tmpfs (rw,nosuid,nodev,size=16278892k,nr_inodes=1048576,inode64)
tmpfs on /run/user/1000 type tmpfs (rw,nosuid,nodev,relatime,size=3255776k,nr_inodes=813944,mode=700,uid=1000,gid=1000,inode64)

```

These filesystems are actually stored in RAM and should be significantly faster. For instance, on my machine, which has an SSD as the main storage, I still get significant speedup when using `/dev/shm`:

```
## SSD
$ time seq 1000000000 > file
real    0m13.687s
user    0m5.623s
sys     0m7.516s

$ time wc -l < file
1000000000

real    0m4.972s
user    0m0.101s
sys     0m2.248s

## tmpfs, RAM
$ time seq 1000000000 > /dev/shm/file

real    0m6.755s
user    0m4.982s
sys     0m1.763s

$ time wc -l < /dev/shm/file
1000000000

real    0m1.024s
user    0m0.120s
sys     0m0.902s

```

As you can see, using `/dev/shm` is significantly faster. And using named pipes can give you a little boost too:

```
$ time ( seq 1000000000 > /dev/shm/file;  wc -l < /dev/shm/file)
1000000000

real    0m8.341s
user    0m5.277s
sys     0m3.044s

$ rm /dev/shm/file && mkfifo /dev/shm/file
$ time ( seq 1000000000 > /dev/shm/file &  wc -l < /dev/shm/file)
1000000000

real    0m6.564s
user    0m5.426s
sys     0m3.965s

```

So, this is essentially a way of doing everything in RAM without needing to code it as doing everything in RAM. You can create these files (the FIFOs, the named pipes) and treat them as files, but since you're creating them in a `tmpfs` filesystem, they're actually in RAM and by using FIFOs you not only get a further speed jump but also minimise the actual data that need to be stored.

---
Tags: coverage, pysam, indexing
---


title: Error with eggNOG-mapper jobs

text:
# Question
Title: Error with eggNOG-mapper jobs
I recently tried using eggNOG-mapper for pan-genome and core genome analysis of files from OrthoFinder, specifically the 'Orthogroup\_Sequences' which contains a huge list of protein files in fasta format. As instructed, the folder was compressed into a tar.gz zip file and uploaded. If the upload was successful, an email is sent to your mail ID where you can check the details of your job and decide whether to START the job or not. Every time I START the job, it fails within a few minutes. Initially I presumed that it was due to the file size and hence reduced the number of sequences submitted, but in vain. Can anyone help me?

eggNOG-mapper : http://eggnog-mapper.embl.de/

# Answer
There is a known bug between eggNOG-mapper2 and eggNOG6. This is described here: https://github.com/eggnogdb/eggnog-mapper/issues/428

In your case, I'm not sure on what strategy you are deploying after OrthoFinder, but there's a chance it's not solvable at present and the best place for information is the eggnog GitHub.
> 2 votes

# Answer
I'd like to answer my own question as I have found the solution for this and thought it might help anyone else trying to use eggNOG-mapper. The prerequisite for uploading the input file with the required collection of .faa files, is that, the file extension has to end with .gz. The site mentions that too. I had made the mistake of uploading a .tar.gz which is what my system produced as a compressed file of my folder of sequences. Although the upload was successful, the runs failed as the file had some unwanted lines owing to the .tar.gz format. I randomly uploaded a concatenated .faa file which worked well and got the results too! The trick is that, the input file has to be strictly a .gz and NOT .tar.gz or it could also be a .faa file with all the sequences.
> 2 votes

---
Tags: orthologues, orthofinder, embl, eggnog, pangenome
---


# Question
Title: How can I programmatically add a Hydrogen 'Atom' to a 'Residue' object?
I know the algorithm for creating a Hydrogen atom and adding to a residue:

```
Point3d  create_hydrogen(Point3d C, Point3d N, Point3d CA, Point3d H)
{
  H.set(N);
  H -= C;
  H.norm();

  Point3d tmp2(N);
  tmp2 -= CA;
  tmp2.norm();

  H += tmp2;
  H.norm(N_H_BOND_LENGTH);
  H += N;

  return H;
}

Point3d  create_hydrogen(Residue prev_residue, Residue residue, Point3d H)
{
  return create_hydrogen(prev_residue.find_atom(" C  "),
                              residue.find_atom(" N  "),
                              residue.find_atom(" CA "),
                                   H);
}

Atom create_hydrogen(Residue prev_residue, Residue residue)
{
  Atom H = make_shared_atom(0, " H  ", 1);
  create_hydrogen(prev_residue, residue, H);
  residue.add(H);
  return H;
}

```

I didn't find this implementation in BioPython.

Is this already implemented in BioPython? If YES, in which module?

If NO, what would be the implementation like in BioPython?

Added image of peptide bond too try to depict the problem at hand, I believe the hydrogen in cause is the one on the peptide bond plane(pale yellow) indentifiedin purple

# Answer
Not sure but here https://biopython.org/docs/1.75/api/Bio.PDB.Residue.html Bio.PDB.Residue module Your can find :

add(self, atom) :

```
      Add an Atom object.

```

Checks for adding duplicate atoms, and raises a PDBConstructionException if so.

Can’t try it just now but who knows it may work

Googling around I found out that if you plan to add Hydrogens to PDB files \[I believe nowadays PDBx is the repository grade standard\] you could do thar following PDB format rules. I am copying them from UCSF Chimera website (https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html#hydrogens):

Hydrogen Atoms

In brief, conventions for hydrogen atoms in version 3.0 PDB format are as follows:

Hydrogen atom records follow the records of all other atoms of a particular residue.
A hydrogen atom name starts with H. The next part of the name is based on the name of the connected nonhydrogen atom. For example, in amino acid residues, H is followed by the remoteness indicator (if any) of the connected atom, followed by the branch indicator (if any) of the connected atom; if more than one hydrogen is connected to the same atom, an additional digit is appended so that each hydrogen atom will have a unique name. Hydrogen atoms in standard nucleotides and amino acids (other than the rarely seen HXT) are named according to the IUPAC recommendations (Markley et al., Pure Appl Chem 70:117 (1998)). Names of hydrogen atoms in HETATM residues are determined in a similar fashion.
If the name of a hydrogen has four characters, it is left-justified starting in column 13; if it has fewer than four characters, it is left-justified starting in column 14.

In the following excerpt from entry 1vm3, atom H is attached to atom N. Atom HA is attached to atom CA; the remoteness indicator A is the same for these atoms. Two hydrogen atoms are connected to CB, one is connected to CG, three are connected to CD1, and three are connected to CD2.

ATOM     10  N   LEU A   2       4.595   6.365   3.756  1.00  0.00           N
ATOM     11  CA  LEU A   2       4.471   5.443   2.633  1.00  0.00           C
ATOM     12  C   LEU A   2       5.841   5.176   2.015  1.00  0.00           C
ATOM     13  O   LEU A   2       6.205   4.029   1.755  1.00  0.00           O
ATOM     14  CB  LEU A   2       3.526   6.037   1.578  1.00  0.00           C
ATOM     15  CG  LEU A   2       2.790   4.919   0.823  1.00  0.00           C
ATOM     16  CD1 LEU A   2       3.803   3.916   0.262  1.00  0.00           C
ATOM     17  CD2 LEU A   2       1.817   4.196   1.769  1.00  0.00           C
ATOM     18  H   LEU A   2       4.169   7.246   3.704  1.00  0.00           H
ATOM     19  HA  LEU A   2       4.063   4.514   2.992  1.00  0.00           H
ATOM     20  HB2 LEU A   2       2.804   6.675   2.065  1.00  0.00           H
ATOM     21  HB3 LEU A   2       4.099   6.623   0.873  1.00  0.00           H
ATOM     22  HG  LEU A   2       2.234   5.353   0.004  1.00  0.00           H
ATOM     23 HD11 LEU A   2       4.648   4.447  -0.148  1.00  0.00           H
ATOM     24 HD12 LEU A   2       3.334   3.331  -0.516  1.00  0.00           H
ATOM     25 HD13 LEU A   2       4.137   3.260   1.052  1.00  0.00           H
ATOM     26 HD21 LEU A   2       0.941   3.892   1.216  1.00  0.00           H
ATOM     27 HD22 LEU A   2       1.522   4.860   2.568  1.00  0.00           H
ATOM     28 HD23 LEU A   2       2.296   3.323   2.188  1.00  0.00           H

```
</code></pre>
> 2 votes

# Answer
The other answer is correct. But I thought I'd give a few pointers in how one can understand how to do something with a given Python module —teach a man to fish kind of thing...

In a Jupyter Notebook or IPython shell (don't), you can do:

* `help(function)` print out the docstring
* `dir(object)` shows all the public and magic attributes and methods —note a class will not show the instance methods
* `type(object)` or `object.__class__` gives you what the class is, but `object.__class__.__mro__` tells you what it inherits.
* in the standard library module `inspect` is the function `inspect.getsource(function)`, which for a notebook is messy due to newlines so needs to be `print_code = lambda fun: print(inspect.getsource(fun))` and `print_code(fun)`.

In your case you read a file...

```
from Bio import PDB
parser = PDB.PDBParser()
structure : PDB.Structure.Structure = parser.get_structure("tryp-cage", "1L2Y.first.pdb")

```

So `parser.get_structure` "magically" makes a structure, eh? This means it can add stuff! As mentioned `print_code(parser.get_structure)` will reveal its secrets, or within it call some `self._private_fun()`, which can be inspected (`self` = `parser`) all the way down the rabbit hole. Until you find a nice constructor that declares nicely an atom instance or whatever you ever fancy searching for —that is the great thing of open source code.

In the specific case of `PDB` module, everything inherits `Entity`, which has the attribute `child_list`, which is what the iterators (`__iter__` is the topmost) read without making a copy. The end of the above rabbit hole may just add a child `Entity`-inheriting class instance to this list, but it also do other things that is worth keeping an eye out.

In the case of `residue.add(atom)`, creating the `Bio.PDB.Atom.Atom` instance is slightly problematic as there's no typehinting or `:type xxx:` declarations in the docstring. But doing `type(atom.bfactor)` etc. on the attributes of an existing atom will reveal how to do it. One the serial number is a PDB thing, which is sequential for the model.

```
last_serial = list(structure.get_residues())[-1].child_list[-1].serial_number + 1
residue.add(PDB.Atom.Atom(name=' H  ', coord=np.array([1,1,1]), bfactor=0., occupancy=1., altloc=' ', fullname=' H  ', serial_number=last_serial+1,element='H'))

```

The coordinates will not be `[1,1,1]`, but the point that you have calculate. Say for `H`, this would be coplanar with `C`, `CA` and `N`, that is 1 and smidge Å away from `N` at an angle not quite in the middle as found in amide bonds.
> 2 votes

# Answer
OK my attempt, input pdb created in **pymol** with `fab AVGLWPT , AVGLWPT , hydro=0` and saved as `avgwpt.pdb` :

```
ATOM      1  N   ALA     1      67.846  51.727  15.161  1.00  0.00           N
ATOM      2  CA  ALA     1      68.522  53.021  15.161  1.00  0.00           C
ATOM      3  C   ALA     1      70.022  52.847  15.161  1.00  0.00           C
ATOM      4  O   ALA     1      70.553  51.730  15.150  1.00  0.00           O
ATOM      5  CB  ALA     1      68.014  53.813  16.378  1.00  0.00           C
ATOM      6  N   VAL     2      70.773  53.896  15.173  1.00  0.00           N
ATOM      7  CA  VAL     2      72.223  53.728  15.173  1.00  0.00           C
ATOM      8  C   VAL     2      72.922  55.066  15.188  1.00  0.00           C
ATOM      9  O   VAL     2      72.296  56.134  15.198  1.00  0.00           O
ATOM     10  CB  VAL     2      72.665  52.877  13.916  1.00  0.00           C
ATOM     11  CG1 VAL     2      74.193  52.762  13.664  1.00  0.00           C
ATOM     12  CG2 VAL     2      72.143  51.420  13.932  1.00  0.00           C
ATOM     13  N   GLY     3      74.212  55.098  15.190  1.00  0.00           N
ATOM     14  CA  GLY     3      74.888  56.392  15.204  1.00  0.00           C
ATOM     15  C   GLY     3      76.388  56.218  15.203  1.00  0.00           C
ATOM     16  O   GLY     3      76.917  55.104  15.192  1.00  0.00           O
ATOM     17  N   LEU     4      77.138  57.267  15.215  1.00  0.00           N
ATOM     18  CA  LEU     4      78.589  57.099  15.215  1.00  0.00           C
ATOM     19  C   LEU     4      79.288  58.437  15.230  1.00  0.00           C
ATOM     20  O   LEU     4      78.660  59.507  15.241  1.00  0.00           O
ATOM     21  CB  LEU     4      79.040  56.251  13.992  1.00  0.00           C
ATOM     22  CG  LEU     4      78.445  54.825  13.832  1.00  0.00           C
ATOM     23  CD1 LEU     4      79.094  54.113  12.637  1.00  0.00           C
ATOM     24  CD2 LEU     4      78.610  53.968  15.098  1.00  0.00           C
ATOM     25  N   TRP     5      80.577  58.469  15.232  1.00  0.00           N
ATOM     26  CA  TRP     5      81.254  59.763  15.246  1.00  0.00           C
ATOM     27  C   TRP     5      82.754  59.589  15.245  1.00  0.00           C
ATOM     28  O   TRP     5      83.281  58.478  15.233  1.00  0.00           O
ATOM     29  CB  TRP     5      80.817  60.551  16.512  1.00  0.00           C
ATOM     30  CG  TRP     5      79.341  60.958  16.536  1.00  0.00           C
ATOM     31  CD1 TRP     5      78.304  60.258  17.187  1.00  0.00           C
ATOM     32  CD2 TRP     5      78.731  62.003  15.872  1.00  0.00           C
ATOM     33  CE2 TRP     5      77.340  61.929  16.130  1.00  0.00           C
ATOM     34  CE3 TRP     5      79.257  63.017  15.031  1.00  0.00           C
ATOM     35  NE1 TRP     5      77.049  60.848  16.950  1.00  0.00           N
ATOM     36  CZ2 TRP     5      76.464  62.876  15.556  1.00  0.00           C
ATOM     37  CZ3 TRP     5      78.372  63.946  14.482  1.00  0.00           C
ATOM     38  CH2 TRP     5      76.996  63.879  14.742  1.00  0.00           C
ATOM     39  N   PRO     6      83.505  60.638  15.257  1.00  0.00           N
ATOM     40  CA  PRO     6      84.959  60.560  15.257  1.00  0.00           C
ATOM     41  C   PRO     6      85.498  60.036  16.599  1.00  0.00           C
ATOM     42  O   PRO     6      84.724  59.735  17.507  1.00  0.00           O
ATOM     43  CB  PRO     6      85.354  62.018  14.978  1.00  0.00           C
ATOM     44  CG  PRO     6      84.201  62.829  15.588  1.00  0.00           C
ATOM     45  CD  PRO     6      82.967  61.980  15.271  1.00  0.00           C
ATOM     46  N   THR     7      86.770  59.913  16.774  1.00  0.00           N
ATOM     47  CA  THR     7      87.263  59.414  18.054  1.00  0.00           C
ATOM     48  C   THR     7      88.771  59.342  18.062  1.00  0.00           C
ATOM     49  O   THR     7      89.436  59.668  17.100  1.00  0.00           O
ATOM     50  CB  THR     7      86.639  58.014  18.376  1.00  0.00           C
ATOM     51  CG2 THR     7      87.080  57.358  19.702  1.00  0.00           C
ATOM     52  OG1 THR     7      85.224  58.110  18.479  1.00  0.00           O
TER
END

```

pic :

code :

```
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
molecule create in pymol with fab AVGLWPT , AVGLWPT , hydro=0

"""

import Bio
print('########## Biopython  VERSION ##########################################')
print('\n\n BIOPYTHON_VERSION : ', Bio.__version__)
print('###################################################################')

from Bio.PDB import (
                    PDBParser,
                    Atom,
                    PDBIO
                    )

structure_code = 'avgwpt'

structure_name = structure_code

parser_pdb= PDBParser(QUIET=False)    ### shows warnings

structure_whole = parser_pdb.get_structure(structure_name, structure_name+".pdb")

last_serial = list(structure_whole.get_atoms())[-1].serial_number

print(" last_serial : ",last_serial)

for model in structure_whole:

    for chain in model:

        for resi in chain:

            print('1 : ', resi)

            print('2 : ',resi.id)

            try :

                print('3 : ',chain[resi.id[1]-1])

                doable = True

            except :

                print('3 : ','residue -1 missing !!!!!!!!!!!!!!!!!!!!!!1')

                doable = False

            print('4 : ',[i for i in resi.get_atoms()])

            if doable:

                print('5 : ','ok')

            else:

                print('5 : ','NO')

            if doable:

                print(chain[resi.id[1]-1]['C'])

                print(chain[resi.id[1]-1]['C'].get_vector(), type(chain[resi.id[1]-1]['C'].get_vector()))

                print(chain[resi.id[1]]['N'].get_vector())

                c = chain[resi.id[1]-1]['C'].get_vector()

                n = chain[resi.id[1]]['N'].get_vector()

                ca = chain[resi.id[1]]['CA'].get_vector()

                nc = n-c

                nca = n-ca

                print('vector BN : ', nc)

                nc.normalize()

                print('vector BN normalize : ', nc)

                print('vector B2B : ', nca)

                nca.normalize()

                print('vector B2B normalize : ', nca)

                print('sum of nc + nca : ' , nsum := nc + nca)

                nsum.normalize()

                nsum**1.023

                nsum += n

                print(nsum)

                print(nsum.get_array(), type(nsum.get_array()))

                resi.add(Atom.Atom(name=' H  ', coord = nsum.get_array(), bfactor=0., occupancy=1., altloc=' ', fullname=' H  ', serial_number = last_serial ,element='H'))

                last_serial += 1

io=PDBIO()

io.set_structure(structure_whole )

io.save(structure_name+"_H_added.pdb" , preserve_atom_numbering = True)

```

output , saved as `avgwpt_H_added.pdb`:

```
ATOM      1  N   ALA     1      67.846  51.727  15.161  1.00  0.00           N
ATOM      2  CA  ALA     1      68.522  53.021  15.161  1.00  0.00           C
ATOM      3  C   ALA     1      70.022  52.847  15.161  1.00  0.00           C
ATOM      4  O   ALA     1      70.553  51.730  15.150  1.00  0.00           O
ATOM      5  CB  ALA     1      68.014  53.813  16.378  1.00  0.00           C
ATOM      6  N   VAL     2      70.773  53.896  15.173  1.00  0.00           N
ATOM      7  CA  VAL     2      72.223  53.728  15.173  1.00  0.00           C
ATOM      8  C   VAL     2      72.922  55.066  15.188  1.00  0.00           C
ATOM      9  O   VAL     2      72.296  56.134  15.198  1.00  0.00           O
ATOM     10  CB  VAL     2      72.665  52.877  13.916  1.00  0.00           C
ATOM     11  CG1 VAL     2      74.193  52.762  13.664  1.00  0.00           C
ATOM     12  CG2 VAL     2      72.143  51.420  13.932  1.00  0.00           C
ATOM     52  H   VAL     2      70.368  54.810  15.182  1.00  0.00           H
ATOM     13  N   GLY     3      74.212  55.098  15.190  1.00  0.00           N
ATOM     14  CA  GLY     3      74.888  56.392  15.204  1.00  0.00           C
ATOM     15  C   GLY     3      76.388  56.218  15.203  1.00  0.00           C
ATOM     16  O   GLY     3      76.917  55.104  15.192  1.00  0.00           O
ATOM     53  H   GLY     3      74.741  54.249  15.182  1.00  0.00           H
ATOM     17  N   LEU     4      77.138  57.267  15.215  1.00  0.00           N
ATOM     18  CA  LEU     4      78.589  57.099  15.215  1.00  0.00           C
ATOM     19  C   LEU     4      79.288  58.437  15.230  1.00  0.00           C
ATOM     20  O   LEU     4      78.660  59.507  15.241  1.00  0.00           O
ATOM     21  CB  LEU     4      79.040  56.251  13.992  1.00  0.00           C
ATOM     22  CG  LEU     4      78.445  54.825  13.832  1.00  0.00           C
ATOM     23  CD1 LEU     4      79.094  54.113  12.637  1.00  0.00           C
ATOM     24  CD2 LEU     4      78.610  53.968  15.098  1.00  0.00           C
ATOM     54  H   LEU     4      76.733  58.181  15.224  1.00  0.00           H
ATOM     25  N   TRP     5      80.577  58.469  15.232  1.00  0.00           N
ATOM     26  CA  TRP     5      81.254  59.763  15.246  1.00  0.00           C
ATOM     27  C   TRP     5      82.754  59.589  15.245  1.00  0.00           C
ATOM     28  O   TRP     5      83.281  58.478  15.233  1.00  0.00           O
ATOM     29  CB  TRP     5      80.817  60.551  16.512  1.00  0.00           C
ATOM     30  CG  TRP     5      79.341  60.958  16.536  1.00  0.00           C
ATOM     31  CD1 TRP     5      78.304  60.258  17.187  1.00  0.00           C
ATOM     32  CD2 TRP     5      78.731  62.003  15.872  1.00  0.00           C
ATOM     33  CE2 TRP     5      77.340  61.929  16.130  1.00  0.00           C
ATOM     34  CE3 TRP     5      79.257  63.017  15.031  1.00  0.00           C
ATOM     35  NE1 TRP     5      77.049  60.848  16.950  1.00  0.00           N
ATOM     36  CZ2 TRP     5      76.464  62.876  15.556  1.00  0.00           C
ATOM     37  CZ3 TRP     5      78.372  63.946  14.482  1.00  0.00           C
ATOM     38  CH2 TRP     5      76.996  63.879  14.742  1.00  0.00           C
ATOM     55  H   TRP     5      81.105  57.620  15.224  1.00  0.00           H
ATOM     39  N   PRO     6      83.505  60.638  15.257  1.00  0.00           N
ATOM     40  CA  PRO     6      84.959  60.560  15.257  1.00  0.00           C
ATOM     41  C   PRO     6      85.498  60.036  16.599  1.00  0.00           C
ATOM     42  O   PRO     6      84.724  59.735  17.507  1.00  0.00           O
ATOM     43  CB  PRO     6      85.354  62.018  14.978  1.00  0.00           C
ATOM     44  CG  PRO     6      84.201  62.829  15.588  1.00  0.00           C
ATOM     45  CD  PRO     6      82.967  61.980  15.271  1.00  0.00           C
ATOM     56  H   PRO     6      83.072  61.539  15.267  1.00  0.00           H
ATOM     46  N   THR     7      86.770  59.913  16.774  1.00  0.00           N
ATOM     47  CA  THR     7      87.263  59.414  18.054  1.00  0.00           C
ATOM     48  C   THR     7      88.771  59.342  18.062  1.00  0.00           C
ATOM     49  O   THR     7      89.436  59.668  17.100  1.00  0.00           O
ATOM     50  CB  THR     7      86.639  58.014  18.376  1.00  0.00           C
ATOM     51  CG2 THR     7      87.080  57.358  19.702  1.00  0.00           C
ATOM     52  OG1 THR     7      85.224  58.110  18.479  1.00  0.00           O
ATOM     57  H   THR     7      87.409  60.156  16.044  1.00  0.00           H
TER      58      THR     7
END

```

pic :

Not sure I got the code above (in question) right, can't read it , had to use comments as trace.

Spent a while banging my head with : HAAD: A Quick Algorithm for Accurate Prediction of Hydrogen Atoms in Protein Structures

and trying to figure out their alghoritm, \[code (https://seq2fun.dcmb.med.umich.edu//HAAD/haad.f90) is in FORTRAN90\] from the paper , but no luck, if anybody wants to help eplaining this bit:

> For constructing the sp2H2 and sp2H1 H-atoms, we first decide on the orientation of the conjugated plane or the aromatic ring with respect to the neighboring heavy atoms; the normal vector of the conjugated plane is determined by taking the cross product of two vectors between the heavy atoms. For the sp2H2 H-atoms (illustrated in Fig. 1b), the normal vector of the conjugated plane is the cross product of the unit vectors B→N and B2→B; then the two H-atoms are placed at positions A1 and H, which are within the conjugated plane respected to the B→N vector with the exact bond angle from CHARMM22 force field. For sp2H1 H-atoms, two conformations are possible. The first is to place the H-atom in the peptide plane as illustrated in Fig. 1b, where A1 and B1 represent the alpha carbon atoms. The position of H in this case is decided by using the same method as the one used to determine the position of the sp2H2 H-atoms while holding the trans-conformation.

it refers to Fig.1 of linked freely available article
> 1 votes

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


# Question
Title: which SNP effect software, such as VEP, annovar or snpEff, allows the user to submit their own annotations in gtf format?
I have a use case where I want to annotate the variants from running Oxford Nanopore's Medaka on a very simple reference, a single sequence starting on an ATG and finishing on a stop codon (no introns).

Which of the variant effect prediction tools would be the easiest to use to programmatically create the .gtf annotation file and build a reference upon which one can then call the variant annotation, and get the different types of variant consequences on the vcf file (synonymous, non-synonymous with the pair of aa's involved, frame-shifts, etc.)

So far I have looked at:

* snpEff: it seems to have a way of building the reference, although it's not straightforward to match the built java .jar with the one in the system.
* annovar: I couldn't find the documentation on how to supply your own annotations, only how to download the human GRCh38 or equivalents.
* VEP: I couldn't find the documentation on how to supply your own annotations, only how to download the human GRCh38 or equivalents.

Any ideas?

# Answer
> 0 votes
I managed to solve the issue by using a slightly different approach to @Ram RS, which is `bcftools csq` instead of `bcftools annotate`. The formatting is via GFF3, and the tool is picky about how to format it, but it turns out Torsten @tseemann had already gone through the pain of matching what the tool expects with the GFF3 flavour, and then it works:

https://github.com/samtools/bcftools/issues/530#issuecomment-377637623

An example gff3 for a single CDS of length=1463 is:

```

seq1    geneX   gene    1       1463    .       +       .       ID=gene:gene1;Name=gene1;biotype=protein_coding
seq1    geneX   mRNA    1       1463    .       +       .       ID=transcript:transcript1;Name=transcript1;Parent=gene:gene1;biotype=protein_coding
seq1    geneX   CDS     1       1463    .       +       0       ID=CDS:cds1;Name=CDS1;Parent=transcript:transcript1;biotype=protein_coding

```

Then call it with:

```
bcftools csq -p m -f $fastaref -g $gff3 $inputvcf -o $outvcf

```

Thank you all for commenting.

# Answer
> 2 votes
All of them have ways you can add your own annotation dataset but I've found the easiest way to be to use bcftools annotate. You can give it a tab-delimited (GTF is tab-delimited, but might need some pre-processing to get relevant values from column-9) or VCF file and have it pick specific columns/attributes, match the locus by any subset of CHR, POS, REF, ALT and add the annotation in a variety of custom manners. vcfanno is another option.

---
Tags: variant-calling, snpeff, annovar
---


# Question
Title: how to add a custom elements into PyMOL visualization?
I would like to add some 3D elements (balls and lines) to the PyMOL visualization to highlight some "holes" in the structure. Say I `fetch 1tqn` and I want to add a yellow ball somewhere (to specific coordinates, say `[-19, -23, -14]`) in the view. Or a bar from - to coordinates.

I was going through wiki, did a online search but have not found a way how to add these simple 3D entities into the scene. I guess it should be trivial and I'm just using wrong queries. Because I need to add lots of these, I need to be able to do that via scripting interface.

Thank you in advance for any help!

# Answer
> 2 votes
my attempt :

```
import pymol

from pymol.cgo import *

from pymol import (

                    cmd ,
                    )

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

pymol.finish_launching()

cmd.fetch('1tqn')

cmd.remove('solvent')

pos = cmd.get_position()

r = 3

print('###############################################################################\n')

print('### POSITION : ' , pos,'  ###\n')

print('###############################################################################\n')

spherelist = [
    [COLOR,    1.000,    1.000,    0.000,
    SPHERE,   pos[0] , pos[1] , pos[2]  , r ],

    ]

for i in spherelist:

    a = str(spherelist.index(i))

    cmd.load_cgo(i, 'sphere_'+a )

cmd.zoom()

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

```

result:

---
Tags: pymol, modelling, 3d-structure
---


# Question
Title: How can I export a pruned phylogenetic tree in nexus format in R?
I have a phylogenetic tree with a lot of different populations and I wanted to removes most of them and apply new names to the labels. I wrote this code:

```
library(ape)
library(phylotools)
tree<- read.tree("Consensus_COMBINED_Cod3genesRatesC_trees.trees")

tree<- drop.tip(tree, tip = c("CeuMist2","CeuMist1","Ceu_PAN_2","Ceu_PAN_1","Ceu_PAN_3","Ceu_PAN_4"))
tree<- drop.tip(tree, tip= c("CcNav1","CcNav2","Cc_Mal_1","Cc_KalR2_2","CcRod1","CcRod2","CcRod3","CcRod4","CcHr1"))
tree<- drop.tip(tree, tip= c("CeST6_4","CeST6_5","CeAA7_1","CeST6_3","CeST6_2"))
tree<- drop.tip(tree, tip= c("Cg_MelR3_1","Cg_MelR6_2","Cg_MelR7_2","Cg_MelR6_1","Cg_MelR7_3","Cg_MelR5_1","Cg_MelR4_2","Cg_MelR1_1","Cg_MelR2_1","Cg_MelR4_1","Cg_MelR7_4","Cg_MelR5_2"))
tree<- drop.tip(tree, tip= c("ChTM5_1","ChTM4_3","ChTM4_1","Ch_MEN_5","Ch_MEN_13","Ch_MEN_8","Ch_LEV_3","Ch_LEV_5","Ch_LEV_4","Ch_MEN_6","Ch_MEN_10","Ch_LEV_2","Ch_MEN_4","Ch_MEN_14","Ch_LEV_1","Ch_MEN_11","Ch_MEN_3","Ch_MEN_9"))
tree<- drop.tip(tree, tip= c("Ci_MSP_2","Ci_MSP_1","Ci_DIA_1","Ci_DIA_2","CiPLAN1","CiPLAN3","CiPLAN2","Ci_MEN_1","Ci_MEN_2","Ci_MEN_12"))
tree<- drop.tip(tree, tip= c("CpGravR3","CpGravR4","CpGio_2","CpGravR5_2","CpGravR2","CpGravR7","CpGio_1","CpGio_3","CpGravR1","CpGravR6_2","CpGravR5_1","CpThe_1","CpSka_1"))
tree<- drop.tip(tree, tip= c("Cn_IG_2","Cn_IG_3","Cn_IG_1","LevantinaFC4779","Assyriella_naegelei"))

newlabels<- c("C. eucineta", "C. codringtonii", "C. elisabethae","C. gittenbergeri", "C. helenae", "C. intusplicata", "C. parnassia")
tree$tip.label<- newlabels

write.nexus(tree, "CodTree.trees")

```

When I use the `write.nexus` command I get the following error:

```
Error in y$tip.label : $ operator is invalid for atomic vectors

```

How can I circumvent this error and export the tree in nexus format?

I will also provide the tree here: https://drive.google.com/file/d/1lYmzW9htwaclC\_yYs3fQom-Z3zbjk8wH/view?usp=drive\_link

# Answer
So after digging I found the solution. What worked for me was adding the `file =`  in the `write.nexus` command. The final line should be `write.nexus(tree, file = "Cod.trees")`
> 3 votes

---
Tags: r, phylogenetics, file-formats
---


title: GWAS in python?

text:
# Question
Title: GWAS in python?
*A question similar to this has also been asked on Biostars*

I understand GWAS is done with R, but are any written in Python? I'm used to python and java, acclimating to R will take a bit of time and my advisor wants some data analyzed. I plan on using R but given the deadline, I want to finish the analysis in Python.

# Answer
> 5 votes
It doesn't (really) matter what language a piece of GWAS software is written in. SAIGE is written in a mix of R and C++. BOLT is written in C++. You shouldn't be writing a GWAS analysis from scratch, as it ends up being a very complex piece of software. You will almost certainly be interacting with the command line, when running a GWAS, so my advice is to find the software which is most suitable for your piece of analysis, not one which is written in a programming language you like.

---
Tags: python, gwas
---


# Question
Title: 1.What size PCR product will be generated using the primers (bold and underlined) in the sequence below?
I would like to know the steps of this question and how to determine the forward and the reverse primers and also the length of the sequence, without using NCBI. If you could explain and give an answer step by step I would appreciate it. Thank you

# Answer
> 3 votes
Without giving too much of the "what they want you to do yourself" bits of this away, here are a few explanations of the image that may help you answer this question:

* As mentioned in the text, the forward and reverse primers are underlined in the image.
* The `-` and `+` symbols below the forward-direction sequence help determine the positions of the bases above and below. A `+` appears every 10 bases.
* The first base in the sequence is at position 1441 (as indicated by the numbers) to the left of the sequence. The first `+` to the right of that first base represents position 1450.
* Primers are usually represented as a sequence in a "5-prime to 3-prime direction", which is almost always the left to right direction of the top sequence in an image like this, and the right to left direction of the bottom sequence.
* The length of a PCR product includes the primer sequence.

For help beyond that explanation... this looks like a canned question. Unless you give more context around what bioinformatics project you're working on and how this question relates to your project (i.e. explain why it's a real-world problem rather than a text-book question), I'm going to leave answers to the specific questions asked as an exercise to the reader.

---
Tags: homework, sequence, primer, pcr, dna
---


title: Making a DNAStringSet object

text:
# Question
Title: Making a DNAStringSet object
I wish to download a fasta sequence from NCBI using "rentrez" package and then save it in an object, then I want to have this object as a DNAStringSet object to work on the sequence using "Biostrings" capabilities.

I leave my incorrect code here so I would be grateful if you could help me out.

```
cov2 <- entrez_fetch(db="nucleotide" , rettype = "fasta", id = "NC_045512")
cov2 <- readDNAStringSet(cov2)
...
Error in .Call2("new_input_filexp", filepath, PACKAGE = "XVector") :
  cannot open file '>NC_045512.2

```

# Answer
Your issue is that `entrez_fetch()` returns the fasta record(s) as a single string with the definition lines and sequences concatenated. Also, newline characters (`\n`) appear in place of line breaks, with double newlines (`\n\n`) separating different records. Here's a function that creates a `DNAStringSet` given a list of nucleotide identifiers without writing each fasta to a file:

```
entrez2dss <-
  function(id_list) {

    require(rentrez)    # v1.2.3
    require(Biostrings) # v2.66.0
    require(stringr)    # v1.5.0

    # fetch all sequences, trim empty elements
    raw_ <-
      entrez_fetch(db = "nucleotide",
                   rettype = "fasta",
                   id = id_list) |>
      str_split_1(pattern = "\n\n") |>
      head(-1)

    # get sequence names
    names_ <-
      gsub(pattern = "\n.*",
           replacement = "",
           x = raw_) |>
      gsub(pattern = "^>",
           replacement  = "")

    # get sequences
    seqs_ <-
      sub(pattern = ".*?\n",
          replacement = "",
          x = raw_) |>
      gsub(pattern = "\n",
           replacement = "")

    # create named DNAStringSet object
    dss <- DNAStringSet(x = seqs_)
    names(dss) <- names_
    dss

  }

seqlist <- c("NC_045512", "OR558592")
cov2 <- entrez2dss(seqlist)

```
> 1 votes

# Answer
The previous error message was because the `entrez_fetch()` function returns the sequence as a character string, not a file path. To save the sequence to a file, I can use the `writeLines()` function. For example, the following code will save the sequence to a file named `cov2.fasta`:

```
# Download the FASTA sequence
cov2 <- entrez_fetch(db = "nucleotide", rettype = "fasta", id = c("NC_045512", "OR558592"))
# Save the sequence to a file
writeLines(cov2, "cov2.fasta")
# Read the FASTA sequence into a DNAStringSet object
cov2_dna <- readDNAStringSet("cov2.fasta")

```
> 1 votes

---
Tags: r, fasta, bioconductor, biostrings
---


# Question
Title: How can I incorporate wildcards in Snakemake in an R script?
I am facing the following issue:

I have a rule in Snakemake that looks something like this:

```
rule somerule:
input:
    tables = expand("results/{tables}_table.txt", tables = ["1", "2"]),
output:
    edit= expand("results/{tables}_edited.txt", tables = ["1", "2"]),
conda:
    "../envs/r.yaml"
script:
    "../scripts/script.R"

```

The script used in this rule is as follows:

```
library()

tab <- read.table(snakemake@input[["tables"]], sep = "\t", header = TRUE)
edit <- some_function(tab)
write.table(edit, file = snakemake@output[["edit"]], sep = "\t", quote = FALSE, row.names = TRUE, col.names = TRUE)

```

When I execute the code, I get this error:

```
Error in file(file, 'rt'): invalid 'description' argument
Called from: read.table -> file
Execution halted

```

However, when I input the data individually (not as wildcards), the script works fine and produces the correct output files. Where could the problem be? My suspicion is that Snakemake is passing all the input data to the R script simultaneously instead of running them individually, but I'm not sure if that's the case and even less sure how to fix it.

# Answer
> 4 votes
Your problem is that Snakemake is passing all the input data to the R script simultaneously instead of running them individually :)

As written your rule takes a list of input files and gives a list of output files (`expand` is just a formatting helper function that returns a list based on its arguments). Instead maybe have one rule that calls the R script with one input and one output, and another rule that asks for all those outputs as an easy way run the whole thing at once.

Maybe something like this? I changed `tables` to `table` here just to emphasize it's one by one.

```
rule all:
    input: expand("results/{table}_edited.txt", table = ["1", "2"])

rule somerule:
    input:
        table = "results/{table}_table.txt"
    output:
        edit = "results/{table}_edited.txt"
    conda: "../envs/r.yaml"
    script: "../scripts/script.R"

```

---
Tags: r, snakemake
---


# Question
Title: Converting VCF format to text for use with PLINK and understanding column mapping
I successfully completed Nature PRS tutorial, which is based on PLINK. Turning to my real data, I downloaded ukb-d-20544\_1.vcf.gz. Now I'm facing the problem that I seem to be unable to use it in PLINK or find the correct data format to download at all, and I am a bit confused about the correct general summary statistics-download procedure.

### 1. Convert .vcf format to a .txt file, PLINK can understand

In the tutorial they use a `.txt.gz` base data file aka summary statistics that looks like this:

```
CHR BP  SNP A1  A2  N   SE  P   OR  INFO    MAF
1   756604  rs3131962   A   G   388028  0.00301666  0.483171    0.997886915712657   0.890557941364774   0.369389592764921
1   768448  rs12562034  A   G   388028  0.00329472  0.834808    1.00068731609353    0.895893511351165   0.336845754096289
...

```

Having this `.vcf.gz` (please see excerpt at the very bottom), I unzipped and, in R, used `MungeSumstats::format_sumstats("ukb-d-20544_1.vcf", ref_genome = "GRCh37")`\[1\], to get something out of it similar to what they use in the .txt file in the tutorial. However, the column names won't match exactly. I know that I can use *BETA* instead of *OR*, and *FRQ* might correspond to *MAF*, but I can't seem to find the imputation information score *INFO*. Here is what `MungeSumstats::format_sumstats` produced:

```
SNP CHR BP  A1  A2  END FILTER  FRQ BETA    SE  LP  N   P
rs12238997  1   693731  A   G   693731  PASS    0.116734    0.00058531  0.000756397 0.357493    117706  0.439042942120051
rs371890604 1   707522  G   C   707522  PASS    0.0983981   0.000793177 0.000848559 0.456023    117706  0.349926634613021
...

```

I need to have the columns corresponding to MAF and INFO as per the tutorial.

The next step in the tutorial involves filtering based on specific criteria, such as MAF and INFO values. However, I can't proceed because these columns are missing.

I would like to primarily use the terminal and only turn to R for statistical analysis. Could you guide me on the correct approach to convert the .vcf file to .txt so that I can continue following the tutorial with PLINK?

### 2. General approach

In addition to addressing the conversion issue, I'm also seeking guidance on the appropriate data sources for genetic analysis. I'm curious about whether OpenGWAS is a reliable source for publication or if there are other recommended sources like dbGaP. Additionally, what is the ideal data format to obtain?

---

The beginning of the received .vcf file looks like this:

```
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##FORMAT=<ID=ES,Number=A,Type=Float,Description="Effect size estimate relative to the alternative allele">
##FORMAT=<ID=SE,Number=A,Type=Float,Description="Standard error of effect size estimate">
##FORMAT=<ID=LP,Number=A,Type=Float,Description="-log10 p-value for effect estimate">
##FORMAT=<ID=AF,Number=A,Type=Float,Description="Alternate allele frequency in the association study">
##FORMAT=<ID=SS,Number=A,Type=Float,Description="Sample size used to estimate genetic effect">
##FORMAT=<ID=EZ,Number=A,Type=Float,Description="Z-score provided if it was used to derive the EFFECT and SE fields">
##FORMAT=<ID=SI,Number=A,Type=Float,Description="Accuracy score of summary data imputation">
##FORMAT=<ID=NC,Number=A,Type=Float,Description="Number of cases used to estimate genetic effect">
##FORMAT=<ID=ID,Number=1,Type=String,Description="Study variant identifier">
##META=<ID=TotalVariants,Number=1,Type=Integer,Description="Total number of variants in input">
##META=<ID=VariantsNotRead,Number=1,Type=Integer,Description="Number of variants that could not be read">
##META=<ID=HarmonisedVariants,Number=1,Type=Integer,Description="Total number of harmonised variants">
##META=<ID=VariantsNotHarmonised,Number=1,Type=Integer,Description="Total number of variants that could not be harmonised">
##META=<ID=SwitchedAlleles,Number=1,Type=Integer,Description="Total number of variants strand switched">
##META=<ID=TotalControls,Number=1,Type=Integer,Description="Total number of controls in the association study">
##META=<ID=TotalCases,Number=1,Type=Integer,Description="Total number of cases in the association study">
##META=<ID=StudyType,Number=1,Type=String,Description="Type of GWAS study [Continuous or CaseControl]">
##SAMPLE=<ID=ukb-d-20544_1,TotalVariants=9938650,VariantsNotRead=0,HarmonisedVariants=9938650,VariantsNotHarmonised=0,SwitchedAlleles=209,StudyType=Continuous>
##contig=<ID=1,length=249250621,assembly=HG19/GRCh37>
##contig=<ID=2,length=243199373,assembly=HG19/GRCh37>
...
##contig=<ID=22,length=51304566,assembly=HG19/GRCh37>
##contig=<ID=X,length=155270560,assembly=HG19/GRCh37>
##contig=<ID=Y,length=59373566,assembly=HG19/GRCh37>
##contig=<ID=MT,length=16569,assembly=HG19/GRCh37>
##contig=<ID=GL000207.1,length=4262,assembly=HG19/GRCh37>
##contig=<ID=GL000226.1,length=15008,assembly=HG19/GRCh37>
...
##contig=<ID=GL000192.1,length=547496,assembly=HG19/GRCh37>
##gwas_harmonisation_command=--json /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/ukb-d-import/processed/ukb-d-20544_1/ukb-d-20544_1_data.json --ref /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/QC/genomes/b37/human_g1k_v37.fasta; 1.1.1
##file_date=2019-11-25T15:12:05.258952
##bcftools_annotateVersion=1.9-74-g6af271c+htslib-1.9-64-g226b4a8
##bcftools_annotateCommand=annotate -a /mnt/storage/home/gh13047/mr-eve/vcf-reference-datasets/dbsnp/dbsnp.v153.b37.vcf.gz -c ID -o /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/ukb-d-import/processed/ukb-d-20544_1/ukb-d-20544_1.vcf.gz -O z /mnt/storage/private/mrcieu/research/scratch/IGD/data/dev/ukb-d-import/processed/ukb-d-20544_1/ukb-d-20544_1_data.vcf.gz; Date=Mon Nov 25 15:37:00 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  ukb-d-20544_1
1       692794  rs530212009     CA      C       .       PASS    AF=0.111693     ES:SE:LP:AF:SS:ID       0.000379398:0.000798869:0.197332:0.111693:117706:1_692794_CA_C
1       693731  rs12238997      A       G       .       PASS    AF=0.116734     ES:SE:LP:AF:SS:ID       0.00058531:0.000756397:0.357493:0.116734:117706:rs12238997

```

# Answer
For the first part of question 2: if there are other recommended sources like dbGaP.

Please refer to this Databases of GWAS summary statistics:

| Database | Content |
| --- | --- |
| GWAS Catalog | GWAS summary statistics and GWAS lead SNPs reported in GWAS papers |
| GeneAtlas | UK Biobank GWAS summary statistics |
| Pan UKBB | UK Biobank GWAS summary statistics |
| GWAS Atlas | Collection of publicly available GWAS summary statistics with follow-up in silico analysis |
| FinnGen results | GWAS summary statistics released from FinnGen, a project that collected biological samples from many sources in Finland |
| Pheweb.jp | GWAS summary statistics of Biobank Japan and cross-population meta-analyses |

`https://www.nature.com/articles/s43586-021-00056-9/tables/3`
> 2 votes

# Answer
> However, the column names won't match exactly. I know that I can use BETA instead of OR, and FRQ might correspond to MAF, but I can't seem to find the imputation information score INFO.

Yes. There is not a standard format that everyone reports GWAS summary statistics in, unfortunately. First, check the paper that you are sourcing it from to see if they already filtered on imputation quality. Do ctrl+f and look in the methods section for imputation information (you may also try searching for imputation tools, like minimac or Beagle).

The good news is that if the sumstats have many alleles, imputation accuracy doesn't matter much. Errors get averaged out and each error only has a small effect since there are many alleles. Usually it is expected that reported sumstats will be pre-QC'd in some way, especially if they lack INFO. It would be safe to proceed without the INFO filter step.

> The next step in the tutorial involves filtering based on specific criteria, such as MAF and INFO values. However, I can't proceed because these columns are missing.

Sumstat formatting is whatever the author wants, and can be a bit annoying to parse. However, FRQ is basically MAF.

See here (ctrl+f "frq"): https://bioconductor.riken.jp/packages/3.14/bioc/vignettes/MungeSumstats/inst/doc/MungeSumstats.html https://academic.oup.com/bioinformatics/article/37/23/4593/6380562

> I would like to primarily use the terminal and only turn to R for statistical analysis. Could you guide me on the correct approach to convert the .vcf file to .txt so that I can continue following the tutorial with PLINK?

Did you try unzipping and simply replacing .vcf with .txt in the filename?

> 2. General approach In addition to addressing the conversion issue, I'm also seeking guidance on the appropriate data sources for genetic analysis. I'm curious about whether OpenGWAS is a reliable source for publication or if there are other recommended sources like dbGaP. Additionally, what is the ideal data format to obtain?

This aligns with your first question. If you use a database, sumstats will be standardized (usually) to the database's common format, so there will be less confusion over inconsistent column names. The GWAS catalog pretty good in my opinion, and all of them are acceptable.
> 2 votes

---
Tags: vcf, public-databases, plink, format-conversion, polygenic-risk-score
---


# Question
Title: Cannot obtain alignment summary after running Bowtie2
I am aligning my Small RNA Seq data with Bowtie 2. Although the alignment performs well, the only information I obtain after finishing running the alignment is the following:

```
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:01
## Time searching: 00:00:01
## Overall time: 00:00:01
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00
## Time loading reference: 00:00:00
## Time loading forward index: 00:00:00
## Time loading mirror index: 00:00:00
## Multiseed full-index search: 00:00:00
## Time searching: 00:00:00
## Overall time: 00:00:00

```

This is my command line to build the index: bowtie2-build -o 3 -t 5 Nuevofastaorina human\_mirs

And this is my command line to do the alignment:

```
bowtie2 -x human_mirs -U 1.lane.clean.gz -N 1 --norc --no-unal --no-head --ignore-quals –-quiet > sr1prueba.sam

```

Does anyone know how to obtain an alignment summary instead of the above information? What I mean by alignment summary is the following:

```
20000 reads; of these:
  20000 (100.00%) were unpaired; of these:
    1247 (6.24%) aligned 0 times
    18739 (93.69%) aligned exactly 1 time
    14 (0.07%) aligned >1 times
93.77% overall alignment rate

```

# Answer
Remove the `--quiet` option, and capture STDERR as well as STDOUT:

```
bowtie2 -x human_mirs -U 1.lane.clean.gz -N 1 --norc --no-unal --no-head --ignore-quals 1>sr1prueba.sam 2>sr1prueba.summary.log

```

**See also:**

`bowtie2` manual: alignment summary
> 5 votes

# Answer
Get samtools.

Then run:

```
samtools stats file.bam

```

and/or

```
samtools flagstat file.bam

```

For example outputs see here.
> 2 votes

---
Tags: rna-seq, sequence-alignment, bowtie2
---


# Question
Title: Where to download a file with major and minor alleles at every position?
I want a list of all variants, i.e. sites which are known to vary between human to human. For example, it should ideally cover all sites in here, but without samples.

I don't want a giant reference panel with many samples, nor a .fa which has no major/minor allele information.

I just want a file with the following information at each locus:

* Location (chromosome position and chromosome)
* Bonus (not necessary): rsID
* Major allele (i.e., A, C, T, or G) at that position
* Minor allele at that position

I'm not picky about how major/minor alleles are defined. Usually they are defined as more common in the general population vs. less common (loosely). Conforming to b37 or b38 would be ideal.

I'm looking for one such file that doesn't have samples in it, but does contain all possible positions (not all VCF files have this)

Where can I download something like this? I like working with VCFs but any format is good.

# Answer
> 3 votes
Sounds like you're looking for the dbSNp files. You haven't specified a species, but you can download the data for human from the dbSNP FTP site. Specifically, you want:

```
wget https://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606/VCF/00-All.vcf.gz

```

This file has lines like:

```
1   10848   rs1271898990    C   G   .   .   RS=1271898990;RSPOS=10848;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP;TOPMED=0.99998407237512742,0.00001592762487257
1   10851   rs1340958514    G   A   .   .   RS=1340958514;RSPOS=10851;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP;TOPMED=0.99996018093781855,0.00003981906218144
1   10854   rs1222740568    G   A,C .   .   RS=1222740568;RSPOS=10854;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP;TOPMED=0.99967348369011213,0.00019113149847094,0.00013538481141692
1   10857   rs1487378566    C   G   .   .   RS=1487378566;RSPOS=10857;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;R5;ASP;TOPMED=0.99933900356778797,0.00066099643221202

```

However, it sounds like you might be underestimating the complexity here. There isn't just one minor allele, many sites have multiple minor alleles. The example above shows one such multiallelic site (you can split these into separate lines using `bcftools`). But you also have cases like:

```
1   10786   rs1464509424    CGGCGCA C   .   .   RS=1464509424;RSPOS=10787;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000200;GENEINFO=DDX11L1:100287102;WGT=1;VC=DIV;R5;ASP;TOPMED=0.99991239806320081,0.00008760193679918
1   10786   rs1258377526    CGGCGCAGGCGCAGACACATGCTAGCGCGTCGGGGTGGAGGCGT    C   .   .   RS=1258377526;RSPOS=10787;dbSNPBuildID=151;SSR=0;SAO=0;VP=0x050000020005000002000200;GENEINFO=DDX11L1:100287102;WGT=1;VC=DIV;R5;ASP

```

Don't expect just one major and one minor allele per location.

# Answer
> 2 votes
I agree with @terdon that perhaps dbSNP is what you may want to look at. The dbSNP155 RefSNP VCF file in GRCh37 version may help. It can be downloaded from https://ftp.ncbi.nlm.nih.gov/snp/latest\_release/vcf/GCF\_000001405.25.gz.

You may also fine the 1000G genome website helpful. https://www.internationalgenome.org/data-portal/data-collection contains much of its release both in `.fa` format and in VCF format. I understand that a large number of samples included increase the file size tremendously, but it is necessary as you're not just looking at what the major and minor allele is, but also their frequency, which may vary greatly between populations. It's important to check which population we're seeing when using their genomic information.

Hope my answer helps. Regards.

---
Tags: vcf, data-download, rsid
---


title: LD clump GRCh38 GWAS results

text:
# Question
Title: LD clump GRCh38 GWAS results
The vignette of R package `ieugwasr` describes a `plink` based wrapper function for LD clumping GWAS data using the *1000 genomes reference dataset* (http://fileserve.mrcieu.ac.uk/ld/1kg.v3.tgz). I would like LD clump GRCh38 based GWAS data. I am under the impression that the previous LD reference dataset is based on GRCh37. Therefore, I kindly ask:

* Is `1kg.v3.tgz` based on GRCh37?
* Can you use GRCh37 based reference to LD clump GRCh38 GWAS results?
* Is there a GRCh38 reference version of 1000 genomes that is split to different super populations (EUR etc)?

# Answer
> Is 1kg.v3.tgz based on GRCh37?

> Is there a GRCh38 reference version of 1000 genomes that is split to different super populations (EUR etc)?

Yes, the genome build for 1000 Genomes phase 3 release is in GRCh37. Luckily ISGR has reanalyzed their original sequencing data based on GRCh38. The data is provided both at the ISGR website https://www.internationalgenome.org/data-portal/data-collection/30x-grch38 and it is also made available by Christopher at PLINK2 website https://www.cog-genomics.org/plink/2.0/resources. However I suspect the data was not split into super populations.

> Can you use GRCh37 based reference to LD clump GRCh38 GWAS results?

I wouldn't recommend doing so. Some variants identified in GRCh37 has been merged into others in GRCh38. Some variants are assigned different rsIDs. The position for the variant may also differ in different genome builds. Therefore clumping GRCh38 GWAS results with GRCh37 based reference may lead to fewer variants correctly mapped, which may influence your furthur analysis.

Hope you find it helpful. Best regards.
> 2 votes

---
Tags: gwas, reference-genome, 1000genomes, ld
---


title: seeing full results on R

text:
# Question
Title: seeing full results on R
running this code on R, using the Biostrings package,

```
vmatchPattern("GAAGACCGAGGCCACGCGGAGTACGAACGAGGGTACAGTCTTC", mySeq)

```

in the console, it writes:

```
...
<22 more elements>

```

How can I see all the results? The `print` function returns the same result.

# Answer
It returns an IRanges object. You can transform it to a data.frame using `as.data.frame()` and then either write as a file with something like `write.table()` to view outside of R or use `View()` in RStudio to explore it.
> 4 votes

---
Tags: r, biostrings
---


# Question
Title: Discrepancy with featurecounts analysis using a forward stranded and reverse stranded protocol
My RNAseq analysis pipeline is as follows: fastqc (read quality is good, some overrepresentation of adaptor sequence) → trimmomatic (trimmed adaptor sequence, qc report after trimming suggests the overrepresented adaptors are gone) → HISAT2 (used a forward stranded alignment protocol) → featurecounts

My HISAT 2 command line:

```
hisat2 --rna-strandness F -x HISAT2/grch38/genome -U "$fastq_file" | samtools sort -o "HISAT2/${sample_name}_aligned.bam"

```

My featurecounts command line:

```
featureCounts -s 2 -a hg38/gencode.v44.chr_patch_hapl_scaff.annotation.gtf -o "quant/${sample_name}_featureCounts.txt" "$bam_file"

```

However, the strange thing is that if I use a forward stranded protocol in feature counts (with the forward stranded alignment I generated with HISAT2), the assignment rate is 5%, with the majority of unassigned reads being having no features. But if I use a reverse stranded protocol in feature counts, the assignment rate is 70%. Does this suggest my data is reverse stranded?

My library was prepared by the NEB directional RNA library prep kit that claims to produce forward stranded libraries, so I’m really confused by this result. How should I proceed with this data? Should I run an alignment with a reverse strand protocol and just proceed with a reverse strand assignment? Or should I proceed with the forward alignment and forward strand assignment data because that’s what the library is supposed to be?

Thank you so much for your help!

# Answer
This kit uses the dUTP method to preserve strandedness. That (to my knowledge) always produces a reverse library. Not sure where you have your information of a forward library from, but we used this kit for many years and it is reverse.
> 3 votes

---
Tags: rna-seq, rna-alignment, featurecounts, strandedness
---


title: Asterisk (*) calculation method

text:
# Question
Title: Asterisk (*) calculation method
Does anyone here know how to calculate a value for the asterisk (\*) code that appears in substitution matrices?
From my observation, to all pairs with a one asterisk, the lowest value from the matrix is set. For `**` it's `1`. Example: BLOSUM62. But, it's just my observation, I couldn't find any source or explanation how it was calculated (for the exemplary BLOSUM62 or any other substitution matrix).

I'd appreciate sources, especially scientific papers.

# Answer
> 2 votes
*Its really easy*. No formula is needed.

Protein matrices are based on observed frequencies from nucleotide alignments for a variety of genes. They are not based on *a priori* formulas (it's not like nucleotides). Later matrices incorporated phylogenetic criteria using maximum likelihood - it's really complicated how it works.

Blosum used odds-ratio (OR) I think without a phylogenetic tree. I suspect raw OR. Explaining with/without a phylogenetic tree is a full lecture in phylogenetic theory, just think it'll be a raw value calculated directly.

So lets just run with raw OR based matrix ...

*So* ... `*` is a stop-codon it carries the lowest weight in the matrix because each protein only has one of them. *So* a stop-codon frequency in comparison to any amino acid is going to be extremely low - because for that aligned position there's only stop-codons (i.e. no other amino acids). Thus its `-4` here (lowest value in the matrix).

**What about the 1?** *However*, when you compare a stop codon frequency ... between proteins, well its occurrence is always 1 to 1, right? Each protein has one stop codon, compare it to another protein ... thats got one stop-codon too. So one million proteins which are nicely aligned to be homologues (thats how its calculated) have one million stop codons ... thus its `1` no matter how many protein there are in the alignment. The alignment position is always the same because stop codons (at least in Blosum) are always homologous.

Once more, comparing a stop codon (in a protein), with the other stop codon in another protein some the same homologous "amino acid" site - no matter which protein it is - its always going to be `1` because its a universal feature of all proteins.

Thats how it works and I would simply cite

> Henikoff, S.; Henikoff, J.G. (1992). "Amino Acid Substitution Matrices from Protein Blocks". PNAS. 89 (22): 10915–10919

---

I've checked the paper and the authors didn't use stop codons, nor did they use any phylogenetic criteria (that was easy to guess). I know they didn't include stop codons because because,

1. The Sigma function for matrix operations has a `20` above it (i.e. 20 amino acids), there's no allowance for stop codons.
2. The matrices they present omit stop codons.

When NCBI are leveraging this matrix they need answers for all possibilities, the idea that `X` (any amino acid) or `*` is omitted would not be cool for them. Thus they could have approached the authors to fill in this information, this would probably mean redoing the alignments to include the stop codon (which the authors wouldn't have liked), or NCBI could have added it after the event (its a pretty safe addition) because the log-ratio of an invariant site is always `1`. The authors/NCBI would just need to be sure the stop codon was in the same position for all protein alignments in Blosum.

You might find the wiki more understandable https://en.wikipedia.org/wiki/BLOSUM

---
Tags: amino-acids, substitution
---


# Question
Title: How to create a phylogenetic tree from diverse mitochondrial genomes
I would like to create a phylogenetic tree for the most species in my dataset.

I'm starting with around 1200 species, but since it's not good practice to align short and long sequences I tried filtering only for the species with the \[15k-18k bp\] range in length, but this leaves me with around 350 species in total.

Those sequences vary a lot in size, here's an histogram of the number of species I have data for  As you can see I have a lot of sequences with length \[12k-15k bp\] and \[15k-18k bp\] (values on the x axis are scaled by a factor of 3000 and then rounded). I would like maximise the species genetic information.
To do so I would like to extract the part of genome that is shared the most among all the species.

Is this achievable by using a regex?
I am not sure on the length and type of the query to be sure of not extracting non-homologous regions as it would create a big problem later on for the phylogenetic tree creation.

I can't make an alignment for all species because of the lack of computational power. That is why I'm trying to subset a piece of DNA that is the most shared (prior to the alignment, but I don't know if that is doable or logic)

My question is: how would you proceed to get the most phylogenetic data out of this dataset?

If you would use a regex, what type of query would you use?

# More Info

### Alignment

When I try to align 10 sequences drawn from the \[15k-18k bp\] range I get this result (the numbers on the right are the TaxID of the species).
We can't see much here so I used an online MSA viewer \[https://alignmentviewer.org/\]

### Data

My data is one Fasta file containing mitochondrial genomes for all the species.

To do this preliminary analysis I only used 10 randomly chosen species with a sequence of length in the range \[15k-18k bp\].

# Answer
> 2 votes
Your overall approach is cool, given caveats (below) to adapt the method. I can fairly authoritatively say your specific methodological approach will not work regarding its phylogenetic strategy without changing it. This is due to a phenomenon where the number of mutations becomes so great it is impossible to recover the ancestral state, basically it randomises the phylogenetic signal, this especially occurs at the 3rd codon nucleotide positions, i.e. 1 in every 3 nucleotides, 33.3% of the entire data set. This violates a central tenant of phylogenetic theory, therefore I am afraid that you cannot build a tree using this approach, but again you can adapt it to work (below).

If you want only to use nucleotides then 16S is the only target available. However, 16S has been a good target for amphibians, therefore just aligning that gene and building a tree is cool. Its perfectly solid given prior work in vertebrates ...

> Chan, K.O., Hertwig, S.T., Neokleous, D.N. et al. Widely used, short 16S rRNA mitochondrial gene fragments yield poor and erratic results in phylogenetic estimation and species delimitation of amphibians. BMC Ecol Evo 22, 37 (2022).

If you wish to use protein coding genes you need to switch amino acids and align those. In terms of conservation there alot of established tools rather than "reinvent the wheel". One program you could consider is **AAcon**, available here. This would give you can empirical basis for deducing which gene is

* the most conserved
* the most represented, i.e. "read-depth", or species representation.

... after its aligned, thereby yielding a genetic marker for phylogenetics. It uses a large number of metrics therefore most scenarios would be accommodated.

COX1 is good target and likely gene to emerge from AAcon analysis.

Building protein phylogenetic trees is a separate question.

---

**Lacking computational power to align**. That is solvable, even on small computers, and you can place that as a separate question stating your computer capacity (RAM and CPUs). It's not difficult to align 1000s even on a small machine, it's simply the strategy to do this. Algorithm performance and accuracy change rapidly.

# Answer
> 2 votes
It looks like these are random mitochondrial sequences pulled from NCBI, which is probably why the alignments aren't working well.

As a starting point, I recommend you make sure that all the mitochondrial DNA sequences are rotated and reverse-complemented as necessary to make sure that the same gene matches position 1 of the sequence. I prefer doing this for COX1, because it's a well-conserved gene. After doing that, it'll be much easier to look at the alignment of that gene alone, which will give you an idea of the complexity of extending the task to the full length of the mitochondrial genome.

I don't expect you'll find a long exact-match sequence that is consistently conserved across all species, even within COX1. It wouldn't surprise me if the longest conserved sequence were less than 16bp, with a few repeats of that sequence throughout the genome. If that is the case, you need to think about what to do with that information. Given the complexity of this intermediate problem you've asked about, it would be a good idea to think first about what you would do after you discovered there were no such sequence.

---
Tags: phylogenetics, sequence-alignment, phylogeny, multiple-sequence-alignment
---


title: GO Enrichment Analysis

text:
# Question
Title: GO Enrichment Analysis
I want to annotate 25,000+ proteins with GO\_ID, TERM, and CATEGORY prior to doing enrichment analysis. Can I extract this information by submitting my sequences? Meanwhile, my species is a fish, and it's a non-model organism. Since my species is a non-model species, I have annotated the proteins from the whole genome by myself using interproscan. There I was able to find GO IDs but there was no term or category which is not helpful for Enrichment analysis.

Seq type:

```
>maker-ctg000210_np1212-snap-gene-104.0-mRNA-1 p/t AED:0.26 eAED:0.29 QI:38|0.66|0.5|1|0.33|0.25|4|0|342
MPCPSKRSRAAKKRKADLRSTEVPSKRSPVAQVDLVSCDEVGRPRKKAHTQSSDAIQSSDCQLQCQLQCQ
LQCQLQSSDAIQANSTNSQTNSTNSLNMCHDNNATDCCFLNESVDFPERSVLMGSFHQGDSRFGNVRNKQ
CGAISLTAVLKSKMKNVWSWETGDLDDVLIKGTVLYSTPRFEVTGVTAVTHATAQMTDSSSCATSSVNYT
VPGKDVKQNVCRPLYVHEVGVTLNRDCASKDENRSCECNVRKPETKQDKQSFNESTTDDVILVNVETHKQ
SLKRNLQATDEVIISDVRFHKQVFKPLTRHDQDALCTCLDLQNENVYLGLMATTDLDDMGCP

```

After annotation:

```
maker-ctg001000_np1212-snap-gene-144.133-mRNA-1       710fdf6795f6f6bfbb960a017c734103        1139    ProSiteProfiles PS50056 Tyrosine specific protein phosphatases domain profile.  354     408     12.21363        T       20-02-2023      IPR000387       Tyrosine-specific protein phosphatases domain   GO:0016311|GO:0016791

```

But I want to use my protein sequences and the following output I need:

```
Seq_id    GO_ID    GO_term    GO_categoty

```

I would appreciate any of your suggestions. Thanks.

# Answer
> 0 votes
It's unclear what you're asking. Can you please give some more information about what you will do with this output once you've found it.

The single annotation you have provided has GO IDs, and using AmiGO, I get the following results:

* GO:0016311
    + Name: dephosphorylation
    + Ontology: biological\_process
* GO:0016791
    + Name: phosphatase activity
    + Ontology: molecular\_function

Bear in mind that GO is a hierarchical annotation, and most GO enrichment analysis programs ignore this hierarchy when carrying out their analysis, which is even more important for a less studied non-model species, where annotation is unlikely to be complete at any level.

Using a standard enrichment program that ignores the hierarchy would require fully populating the GO tree for each protein. For example, using the inferred tree view from AmiGO for that one protein:

* GO:0016311 -\> GO:0008150|GO:0009987|GO:0008152|GO:0044237|GO:0006793|GO:0006796|GO:0016311
* GO:0016791 -\> GO:0003674|GO:0003824|GO:0016787|GO:0016788|GO:0042578|GO:0016791

Consider this when thinking about the type of output you want. This will not be a single GO term/category per sequence; each sequence will have multiple GO IDs, multiple terms, and multiple categories, with many overlaps and lots of gaps.

Further information about exploring the GO database (both automated and manual) can be found here:

https://amigo.geneontology.org/amigo/software\_list

---
Tags: go-enrichment
---


title: How to convert BED to GFF3

text:
# Question
Title: How to convert BED to GFF3
I would like to convert a `BED` format to `GFF3`.

The only useful tool that I could find via a google search seems to be Galaxy, and I do not feel very comfortable with online tools, plus the webserver is currenlty under maintenance.

Does anyone knows about a command-line tool that can handle this conversion?

Edit: here are some lines of my BED file:

```
$ head -4 last_minion-r7_sort.bed
211000022278137 175 211 8e5d0959-7cdb-49cf-9298-94ed3b2aedb5_Basecall_2D_000_2d 42  +
211000022279134 0   503 e8a9c6b8-bad2-4a7e-97d8-ca4acb34ff70_Basecall_2D_000_2d 69  -
211000022279134 24  353 e258783d-95a3-41f5-9ad5-bb12311dbaf4_Basecall_2D_000_2d 45  -
211000022279134 114 429 26601afb-581a-41df-b42b-b366148ea06f_Basecall_2D_000_2d 100 -

```

The bed file thus has 6 columns as for now: chromosome, start coordinate, end coordinate, read name, score, strand. This file was obtained from conversion of `MAF` format (as output of alignment of RNA-seq reads to reference genome, using `LAST`) converted to `SAM` using `maf-convert`, then to `BAM` using `samtools`, finally to `BED` using `bedtools`.

The aim of my conversion is basically to convert `SAM` -\> `GTF`, for post-processing. Since there is no straightforward way to do this, I am going through steps, the only way to do this in my knowledge is : `SAM` -\> `BAM` -\> `BED` -\> `GFF3` -\> `GTF` but for now I am stuck in the `BED` -\> `GFF3` part.

# Answer
> 8 votes
Galaxy has API and API-consuming libraries (such as BioBlend) that will allow you to interactively script against it without opening the graphical interface at all.

However you can also take almost any tool out of Galaxy and use it independently since everything is open source. The converter you mentioned is available as a Python script here and the tool 'wrapper' which you can use to understand how to invoke the python script is next to it.

# Answer
> 8 votes
To answer the question as asked, for people googling. For BED6, in python:

```
#contigs.tsv contians chromosome names and lengths in two columns
for line in open("contigs.tsv"):
   fields = line.strip().split("\t")
   print fields[0], ".", "contig","1",str(fields[1]), ".", "+", ".", "ID=%s" % fields[0]

for line in open("my_bed_file.bed"):
   fields = line.strip().split("\t")

   # note: BED is 0-based, half-open, GFF is 1-based, closed
   start = str(int(fields[1]) + 1)
   print fields[0], "bed", "interval", fields[1], fields[2], fields[4], fields[5], ".", "ID=%s;parent=%s" % (fields[3], fields[0])

```

For bed12, in python:

```
#contigs.tsv contians chromosome names and lengths in two columns
for line in open("contigs.tsv"):
   fields = line.strip().split("\t")
   print fields[0], ".", "contig","1",str(fields[1]), ".", "+", ".", "ID=%s" % fields[0]

for line in open("my_bed12.bed"):

   fields = line.strip().split("\t")
   contig = fields[0]
   # note: BED is 0-based, half-open, GFF is 1-based, closed
   start= int(fields[1]) + 1)
   end=fields[2]
   name=fields[3]
   score=fields[4]
   strand=fields[5]
   print contig, "bed", "interval", str(start), end, score, strand, ".", "ID=%s;parent=%s" % (name, contig)

   block_starts = map(int,fields[11].split(","))
   block_sizes = map(int, fields[10].split(","))

   for (block, (bstart, blen)) in enumerate(zip(block_starts, block_sizes)):
      bend = start + bstart + blen
      print contig, "bed", "block", str(start + bstart), str(bend), score, strand, ".", "ID=%s_%i;parent=%s" %(name, block, name)

```

# Answer
> 5 votes
To convert BAM to GTF, which is the best way to get a file to compare with cuffcompare:

```
import pysam

bamfile=pysam.AlignmentFile("my_bam_file.bam")

for alignment in bamfile.fetch():

    if alignment.is_unmapped:
        continue

    contig = bamfile.get_reference_name(alignment.reference_id)
    name = alignment.query_name

    if alignment.is_reverse:
        strand = "-"
    else:
        strand = "+"

    for start, end in alignment.getblocks():
        # note: BED is 0-based, half-open, GFF is 1-based, closed
        print contig, "BAM", "exon", str(start + 1), str(end), "0", strand, ".", 'gene_id "%s"; transcript_id "%s";' % (name, name)

```

NB, this will also work with a SAM file as long as it is headered.

# Answer
> 4 votes
Bioconductor makes this so easy. It does the coordinate conversion on import.

```
library(rtracklayer)

## import the bed file
bed.ranges <- import.bed('regions.bed')

## export as a gff3 file
export.gff3(bed.ranges,'regions.gff3')

```

And people wonder why R is so popular for bioinformatics...

Also if needed, you could go straight from **BAM** file to **gff3** in R as well.

# Answer
> 3 votes
I recently developed bed2gff to quickly convert .bed files to a gff3 format, a tool written in Rust. Could be of help here!

# Answer
> 2 votes
If anybody needs a simple bed to gff converter, simple in the sense that only takes a bed with a single type of feature, lets say peaks and translate it to GFF I recently wrote this, is Bash and as far as I know all the tools it uses should be present in most Linux / Mac operating systems.

https://github.com/UBMI-IFC/bed2gff/tree/main

It cannot deal, however with multi featured Bed files as it ignores everything from the 7th column and beyond.

# Answer
> 0 votes
You could also use `agat_convert_bed2gff`. `agat` is generally a very useful tool when needing to convert from / to / between gff formats. See here how to use it to convert from bed to gff. https://agat.readthedocs.io/en/latest/tools/agat\_convert\_bed2gff.html

---
Tags: software-recommendation, bed, format-conversion, gff3
---


# Question
Title: bcbio run not running samples in parallel
Sorry for the cross post but I'm not getting a response on bcbio github:

I'm trying to run bcbio on 20 samples and it's only running one at a time as it steps through the yaml seemingly. I'm using slurm. I was wondering if someone can look at this bcbio snippet and if they see anything wrong?

Slurm script

```
#!/bin/bash

#SBATCH --job-name=bcbiopipeline_GPC3_PDx
#SBATCH -t 1-00:00
#SBATCH -c 1
#SBATCH -p core -n 32
#SBATCH --mem-per-cpu=40G
#SBATCH --error=job.%J.err
#SBATCH --output=job.%J.out
#SBATCH --export=ALL

echo "Load bcbio..."
module use /projects/ngs_prodifx/local/software/modules/
module load bcbio-nextgen/1.2.9

bcbio_nextgen.py seqc.yaml -n 72 -t ipython -s slurm -q core -r t=0-72:00 -r conmem=20 --timeout 3000 --tag "rna"

```

yaml snippet:

```
resources:
  default:
    memory: 30G
    cores: 72
details:
- algorithm:
    aligner: star
    disambiguate: mm10
    expression_caller:
    - salmon
    - kallisto
    quality_format: standard
    strandedness: auto
    transcriptome_fasta: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/GRCh38.primary_assembly.genome_GPC3CART.fa
    transcriptome_gtf: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/gencode.v29.annotation_GPC3CART.gtf
  analysis: RNA-seq
  description: 10_LI6612_Baseline_8743
  files:
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-10_R1_001.fastq.gz
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-10_R2_001.fastq.gz
  genome_build: hg38
  metadata:
    panel: Baseline
- algorithm:
    aligner: star
    disambiguate: mm10
    expression_caller:
    - salmon
    - kallisto
    quality_format: standard
    strandedness: auto
    transcriptome_fasta: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/GRCh38.primary_assembly.genome_GPC3CART.fa
    transcriptome_gtf: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/gencode.v29.annotation_GPC3CART.gtf
  analysis: RNA-seq
  description: 11_LI6612_Baseline_8765
  files:
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-11_R1_001.fastq.gz
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-11_R2_001.fastq.gz
  genome_build: hg38
  metadata:
    panel: Baseline
- algorithm:
    aligner: star
    disambiguate: mm10
    expression_caller:
    - salmon
    - kallisto
    quality_format: standard
    strandedness: auto
    transcriptome_fasta: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/GRCh38.primary_assembly.genome_GPC3CART.fa
    transcriptome_gtf: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/gencode.v29.annotation_GPC3CART.gtf
  analysis: RNA-seq
  description: 12_LI6612_Baseline_8770
  files:
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-12_R1_001.fastq.gz
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-12_R2_001.fastq.gz
  genome_build: hg38
  metadata:
    panel: Baseline
fc_name: seqc
upload:
  dir: ../output/final

```

# Answer
> 0 votes
I understand (see **Slurm array job solution** below), what I can definitely help with is the automated generation of the `seqc.yaml` file. This avoids the repetition (and risk of error) in copy and pasting. The following one-liner exactly reproduces the `seqc.yaml` file and needs the `chunk.yaml` below.

```
for ((x=10;x<=13;x++)); do while IFS='' read -r a; do echo "${a//SM-745-10/SM-745-${x}}"; done < chunk.yaml >> seqc.yaml; done

```

This could be generalised to 3 command lines in bash,

```
y='SM-745-10'
z='M-745'
for ((x=10;x<=13;x++)); \
do while IFS='' read -r a; do echo "${a//${y}/${z}${x}}";\
done < chunk.yaml >> seqc.yaml; done

```

Here everything from `for` is a bash one-liner. Obviously `((x=10,x<=13;x++))` would be `((x=1,x<=20;x++))` for 20 fastq forward and reverse files. This will ensure the separate files are correctly named in the yaml file, in this example, viz.

```
# first pair
SM-745-10_R1_001.fastq.gz
SM-745-10_R2_001.fastq.gz
# second pair
SM-745-11_R1_001.fastq.gz
SM-745-11_R2_001.fastq.gz
# third pair
SM-745-12_R1_001.fastq.gz
SM-745-12_R2_001.fastq.gz

```

The `chunk.yaml` is,

```
details:
- algorithm:
    aligner: star
    disambiguate: mm10
    expression_caller:
    - salmon
    - kallisto
    quality_format: standard
    strandedness: auto
    transcriptome_fasta: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/GRCh38.primary_assembly.genome_GPC3CART.fa
    transcriptome_gtf: /swiftcache/ngs/oncology/analysis/todd_creasy/working/reference/gencode.v29.annotation_GPC3CART.gtf
  analysis: RNA-seq
  description: 10_LI6612_Baseline_8743
  files:
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-10_R1_001.fastq.gz
  - /swiftcache/ngs/oncology/analysis/todd_creasy/working/fastq/SM-745/SM-745-10_R2_001.fastq.gz
  genome_build: hg38
  metadata:
    panel: Baseline
-

```

---

### bcbio

To summarise the response from the Git support team is there's conflict occurring ...

1. The issue is here for the slurm script ...

```
#SBATCH -p core -n 32

```

I assume this is `-N 1` ,i.e. \>= 32 core node, if it was 16 cores per node then `-N 2`.

However, **bcbio** is taking care of this and the restriction `32` is conflicting with `72` from `bcbio`, i.e. it creates its own Slurm controller script.

2. in **bcbio** yaml

```
    default:
    memory: 30G
    cores: 72

```

2a. Firstly `72` is specified in Slurm's CLI here ... `-n 72`

2b. This isn't needed because,

```
#SBATCH --mem-per-cpu=40G

```

`--mem-per-cpu` This is a very recent feature BTW.

---

### Slurm array job solution

What I considered was using the `chunk.yaml` above with an array job - this solution isn't needed because **bcbio** is dynamically reading and parallelising `seqc.yaml` (maybe deploying a wrapper for `sinfo`, `squeue`). In the solution below what I'm doing is placing each fastq pair directly into the Slurm queue. This would then parallelise everything directly. Again whilst it would work, it's unnecessary.

Thus for all 20 fastq files assuming `SM-745-1_R1_001.fastq.gz` ... `SM-745-20_R1_001.fastq.gz`.

```
#!/bin/bash

#SBATCH --job-name=bcbiopipeline_GPC3_PDx
#SBATCH -t 1-00:00:00 # 1- is one day. 00:00:00 is a wall clock. Does "1-00:00" work?
#SBATCH -c 1
#SBATCH --array=1-20 # this is now in $SLURM_ARRAY_TASK_ID
#SBATCH --mem-per-cpu=40G
#SBATCH --error=job.%J.err
#SBATCH --output=job.%J.out
#SBATCH --export=ALL

echo "Load bcbio..."
module use /projects/ngs_prodifx/local/software/modules/
module load bcbio-nextgen/1.2.9

# SM-745-10 and SM-745 could be assigned as variables y and z as above

while IFS='' read -r a; \
do echo "${a//SM-745-10/SM-745-$SLURM_ARRAY_TASK_ID"; \
done < chunk.yaml > seqc.yaml; done
bcbio_nextgen.py seqc.yaml -n 72 -t ipython -s slurm -q core -r t=0-72:00 -r conmem=20 --timeout 3000 --tag "rna"

```

**Note** the `-n 72` parameter, i.e. number of cores deployed, will likely be reduced under this solution.

**Alternative** The alternative is to issue a slight modification of the bash one-liner described above to produce 3 (or 20) separate yaml files via bash `while` referring the different fastq paired files and use array job to cycle through them and delete each in turn when done.

---
Tags: python, slurm, yaml
---