# Question
Title: finding motifs in fasta file with multiple sequences
I am able to find matching motifs on a single sequence, as well as their position as shown below :

```
demo = "ABCLXLYLZLABC"
motif="L.L.L"
for element in re.finditer(motif,demo):
    print(element.start(),element.end,element.group())

```

This works and gives me the desired output which is:  `3 8 LXLYL`

Now I want to know how can I then do this same process but on a multifasta file. As done in this example script from https://www.biostars.org/p/14305/ .

The protein sequences on the file are :

```
> P1
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL
>P2
 MAEVRKFTKRLSKPGTAAELRQSVSEAVRGSVVLEKAKLVEPLDYENVITQRKTQIYSDP
 LRDLLMFPMEDISISVIGRQRRTVQSTVPEDAEKRAQSLFVKECIKTYSTDWHVVNYKYE
 DFSGDFRMLPCKSLRPEKIPNHVFEIDEDCEKDEDSSSLCSQKGGVIKQGWLHKANVNST
 ITVTMKVFKRRYFYLTQLPDGSYILNSYKDEKNSKESKGCIYLDACIDVVQCPKMRRHAF
 ELKMLDKYSHYLAAETEQEMEEWLIMLKKIIQINTDSLVQEKKDTVEAIQEEETSSQGKA
 ENIMASLERSMHPELMKYGRETEQLNKLSRGDGRQNLFSFDSEVQRLDFSGIEPDVKPFE
 EKCNKRFMVNCHDLTFNILGHIGDNAKGPPTNVEPFFINLALFDVKNNCKISADFHVDLN
 PPSVREMLWGTSTQLSNDGNAKGFSPESLIHGIAESQLCYIKQGIFSVTNPHPEIFLVVR
>P3
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

```

The script that was used to find the motif is below, the problem now is that it does not show the position where the motif is found. Could you assist in how one could do so.

```
import re
infile=open("seq.fasta",'r')
out=open("results.csv",'w')
pattern=re.compile(r"(P[A-Z]{2}P)")
for line in infile:
  line = line.strip("\n")
  if line.startswith('>'):
     name=line
else:
   s = re.findall(pattern,line)
   print '%s:%s' %(name,s)
   out.write('%s:\t%s\n' %(name,s))

```

Their script seems to have worked because the desired output was

```
>p1:
PGCP

>p1:
PHCP,
PKCP

```

The problem with their script is that, it does not provide the position where the sequence was found. Could you please assist on how that script can be modified? I want to try it but I want my output to have (1) Name Of Sequence, (2) The Motif Found , (3) The Position

# Answer
This is easy just replace `findall` with `finditer` in the above code. The index will then appear where the match first appears and where it stops. However, you will to parse the `finditer` object which in the code is `s`

Then replace

```
print '%s:%s' %(name,s)

```

with

```
for match_obj in s:
    # print each re.Match object
    print(match_obj)

    # extract each matching number
    print(match_obj.group())

```

That should work. When it does you simply pass the `match_obj` and `match_obj.group()` to a file by replacing `s` with `match_obj, match_obj.group()`

Thus replace

```
  out.write('%s:\t%s\n' %(name,s))

```

with

```
  out.write('%s:\t%s\t%s\n' %(name,match_obj,match_obj.group()))

```

You might need to do an append BTW.

thus replace

```
out=open("results.csv",'w')

```

with

```
import os
if os.path.exists("results.csv"):
    os.unlink ("results.csv")
out=open("results.csv",'a')

```
> 0 votes

---
Tags: proteins, python, motifs, multi-fasta
---


# Question
Title: Compare and Reorganize Fasta Headers Python
I want to compare the headers from the fasta file to file1, and if there's a match, reorganize the header and put the match first. If there's no match between fasta file and file1, look at file2 and if there's a match, put that name first on the header.

The goal is to organize the fasta file header with the similar match listed as the first entry on the header. How can I do this with this code or is there another way to write it in python or biopython?

`BioPython` solutions welcome.

Current files look like this

seqs.fasta:

```
>ant| bee| fly
AGCT...
>bird| eagle| hawk| vulture
GATCG...

```

file1.txt:

```
mouse
hawk
cow

```

file2.txt:

```
crane
fly

```

The output file should look like this, so that the order of the header starts with the match from file 1 or file 2:

```
> fly| ant| bee
AGCT...
> hawk| bird| eagle| vulture
AGCT...

```

My code is:

```
fasta= open('seq.fasta', 'r')
one= open('file1.txt' , 'r')
two= open('file2.txt' , 'r')
output= open('newseq.fasta', 'w')

for line in fasta:
    if line.startswith('>'): # getting just the header
        header= line.strip().split("|") # there are several entries on the header separated by the |
        for name in one:
           if name in header: # looking to see if file1.txt has match in fasta file
              header.insert(0,line.pop(header.index(name))) # insert the name from file1 as the first header entry
        else: # if no match in file1 look at file2
            for name in two:
               if name in header:
                    header.insert(0,line.pop(header.index(name))) # insert the name from file2 as first header entry
    else:
output.write(line)

```

# Answer
There is nothing wrong with the concept of the code, that is exactly how I would do it. The problem is here:

```
    else:
output.write(line)

```

Corrections

```
import os
output = 'newseq.fasta'
if os.path.exists(output):
    os.unlink(output)

for line in fasta:
    if line.startswith('>'): # getting just the header
        header= line.strip().split("|") # there are several entries on the header separated by the |
        for name in one:
           if name in header:
              line = header.insert(0,line.pop(header.index(name))) # insert the name from file1 as the first header entry
        else: # if no match in file1 look at file2
            for name in two:
               if name in header:
                    line = header.insert(0,line.pop(header.index(name)))

# The indentation must be absolutely correct
# otherwise the code falls over. Must be one indent from
# the first for statement
    with open(output, 'a') as fout:
        fout.write(line)

```

I assume that `fasta`, `one` and `two` are loaded successfully.
> 1 votes

---
Tags: python, fasta, phylogenetics, phylogeny
---


# Question
Title: Calling isoforms from long read data generated from partially degraded RNA
What will be the best tool to call isoforms from long read data generated from partially degraded RNA. By mistake we processed some samples with poor quality RNA to generate long read. Now we are looking for some tool which can merged the reads for the potential degraded transcripts. I think majority of the available tools like Squanti, Talon, bamboo, flair considered independent read as full length. Do the community has any suggestions?

Thanks

# Answer
I haven't used these long-read RNA-seq tools so I could be wrong here. Transcript reconstruction combines reads from the same transcript and should be robust to partial transcripts. In addition, Nanopore direct RNA-seq produces a lot of partial transcripts. I tend to believe many transcript reconstruction tools can deal with them.

As to what tool to use, StringTie is a popular one and has been shown to be better than Flair. IsoQuant seems promising. Its developers argued it is more accurate than StringTie, Bambu, Talon and Flair. Squanti only works for iso-seq.
> 2 votes

# Answer
You could have a go using LAST. I've been developing and using a protocol for mapping cDNA reads that is based around LAST-mapping since 2017. The most recent RY model addition in LAST seems to clean up the majority of the random gene counting noise issues that we were having.

See the gene counting part of the protocol here. If you're working from direct/native RNA rather than cDNA (or don't care about strand-specific mapping), then the strand correction in previous steps is not needed:

https://dx.doi.org/10.17504/protocols.io.5qpvonn2bl4o/v17
> 1 votes

---
Tags: rna-seq, assembly, long-reads, isoform
---


title: motif searching in a multi-fasta file

text:
# Question
Title: motif searching in a multi-fasta file
I am having trouble trying to search for motifs in a multi fasta file. I have used two techniques one gives me the name of the sequence where the motif is found but doesn't give me the the motif and its position. The other method does not return anything. Please see the code below, any assistance will be highly appreciated.

This is the protein sequence: (Please note I have only put in the first 3 sequences, the actual file I am working with has 100 sequences)

```
>P1
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQXXTXKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

>P2
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQTVTKKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

>P3
MSDNGPQNQRNAPRITFGGPSDSTGSNQNGERSGARSKQRRPQGLPNNTASWFTALTQHG
KEDLKFPRGQGVPINTNSSPDDQIGYYRRATRRIRGGDGKMKDLSPRWYFYYLGTGPEAG
LPYGANKDGIIWVATEGALNTPKDHIGTRNPANNAAIVLQLPQGTTLPKGFYAEGSRGGS
QASSRSSSRSRNSSRNSTPGSSRGISPARMAGNGGDAALALLLLDRLNQLESKMSGKGQQ
QQGQTVTKKSAAEASKKPRQKRTATKAYNVTQAFGRRGPEQTQGNFGDQELIRQGTDYKH
WPQIAQFAPSASAFFGMSRIGMEVTPSGTWLTYTGAIKLDDKDPNFKDQVILLNKHIDAY
KTFPPTEPKKDKKKKADETQALPQRQKKQQTVTLLPAADLDDFSKQLQQSMSSADSTQA*

```

The first one is:

```
infile=open("sequence.fasta",'r')
out=open("Result.csv",'w')
pattern=re.compile(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]")
for line in infile:
  line = line.strip("\n")
  if line.startswith('>'):
     name=line
else:
   s=re.finditer(pattern,line)
   print('%s:%s' %(name,s))
   out.write('%s:\t%s\n' %(name,s))

```

This one above returns, the output below. It gives me the name of the sequence which is P100, but it does not give me (1) The motif found (2) the position of the motif. And the Results.csv file that it generates is blank.

```
#output
>P100:<callable_iterator object at 0x000001D0611B7AF0

```

I then tried a different technique:

```
infile=open("sequence.fasta",'r')
open=open("result.csv",'w')
pattern=re.compile(r"(P[A-Z]{2}P")
for line in infile:
  line=line.strip("\n")
  if line.startswith('>'):
     name=line
else:
   s = re.finditer(pattern,line)
   for match_obj in s:
       print(match_obj)
       print(match_obj.group())
  out.write('%s:\t%s\t%s\n' %(name, match_obj,match_obj.group()))

```

The one above does not return anything at all. the Results.csv file is also blank. I'm still new in this python, and some of the techniques I used here I found on bioinformatics stackexchange.

The output I desire, is where I can get (1) Name of sequence, (2) Motif found and (3) position of motif (where it starts and ends) as shown in the example below in the context of every sequence in the multi fasta file

```
>P100
motif:
position:

```

Any format is fine as long as I can get something similar to the format above. Could I also be assistet with the indention as that cofuses me.

# Answer
First of all, you can't do this this way. Fasta sequences can, and indeed yours do, span multiple lines so in order to find a pattern, you need to join those lines into one to make sure you don't miss any that are found across more than one line.

Next, as a general rule, it is a bad idea to hardcode output file names. That means the user cannot control the output, and that any existing files will be overwritten, so just print to standard output and redirect to a file:

```
foo.py > out_file

```

You also don't want to hardcode the name of the input file since that just makes the script harder to use. Instead, read the input file name as an argument. With all this in mind, here is a working version of your script:

```
#!/usr/bin/env python3
import re
import sys

def print_matches(seq_name, matches):
     if matches:
          for match in matches:
               print(name,match.start(), match.end(),match.group(),sep='\t' )

pattern = re.compile(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]")
seq = ""

try:
     with open(sys.argv[1],'r') as infile:
          for line in infile:
               line = line.strip("\n")
               if line.startswith('>'):
                    name=line.lstrip('>')
                    if seq != "":
                         matches = re.finditer(pattern, seq)
                         print_matches(name, matches)
                         seq = ""
               else:
                    seq += line
          else:
               matches = re.finditer(pattern, seq)
               print_matches(name, matches)

except IOError as e:
     print("Failed to open file %s with error %s" % (sys.argv[1], e))
     sys.exit(1)

```

Running that on your example gives:

```
$ foo.py sequence.fasta
P2  44  53  LPNNTASWF
P2  122 131 YGANKDGII
P2  209 218 MAGNGGDAA
P2  219 228 ALLLLDRLN
P2  304 313 AQFAPSASA
P2  391 400 VTLLPAADL
P2  410 419 MSSADSTQA
P3  44  53  LPNNTASWF
P3  122 131 YGANKDGII
P3  209 218 MAGNGGDAA
P3  219 228 ALLLLDRLN
P3  304 313 AQFAPSASA
P3  391 400 VTLLPAADL
P3  410 419 MSSADSTQA
P3  44  53  LPNNTASWF
P3  122 131 YGANKDGII
P3  209 218 MAGNGGDAA
P3  219 228 ALLLLDRLN
P3  304 313 AQFAPSASA
P3  391 400 VTLLPAADL
P3  410 419 MSSADSTQA

```

The output is 4 tab-separated columns per match with the sequence name, the match's start position, its end position and the actual sequence matched.
> 3 votes

# Answer
Read the sequences with Biopython, reverse translate the sequences and the pattern peptides, construct pattern regexes, search (`find_matches()`) for the patterns for each sequence in a loop using DNA Chisel.

*Disclaimer: I'm the maintainer of DNA Chisel*
> 1 votes

---
Tags: proteins, python, motifs, multi-fasta
---


# Question
Title: bedtools coverage - Report the depth at each position in each A feature
I am using bedtools coverage to compute the sequencing depth at every positions of a chromosome but it didn't work as I expected. Instead it reported 0 coverage at every positions.

This is how I did it :

I made a `chr1.bed` bedfile :

```
chr start   end
1   1   158534110

```

I ran this command :

```
coverageBed -a -d chr1.bed -b aln.bam -s > coverage.perposition_chr1.tsv

```

However I don't find any position in the output that is actually covered :

```
chr start end 1-pos depth
1   1   158534110   1   0
1   1   158534110   2   0
1   1   158534110   3   0
1   1   158534110   4   0
1   1   158534110   5   0
1   1   158534110   6   0
1   1   158534110   7   0
1   1   158534110   8   0
(...)
1   1   158534110   158534109   0

```

I must do something wrong because I already computed depth at specific pattern positions using the same `aln.bam` and it was working fine...

**Suggestions and tries :**

* Swapping `-a` & `-d` parameters produces the same results
* Adding a strand column to the bedfile still gives 0 coverage.

`chr1.bed`

```
1       1       158534110       .       0       +
1       1       158534110       .       0       -

```

`output`

1       1       158534110       .       0       +       1       0
1       1       158534110       .       0       +       2       0
(...)
1       1       158534110       .       0       -       158534108       0
1       1       158534110       .       0       -       158534109       0
```
</code></pre>

# Answer
> 2 votes
If you prefer to have a more concise report than a per-base one, bedtools genomecov could be a better choice.

Its `-bga` option allows you to report consecutive bases that have the same coverage as a range in a single line.

In your case, you could achieve this (unfortunately they only offer the strandedness option for a specified strand) by running the following command:

```
bedtools genomecov -bga -ibam aln.bam -g chr1.bed -strand + > coverage.perposition_chr1.tsv

```

Output will look like:

```
chr1  554304  554309  5
chr1  554309  554313  6
chr1  554313  554314  1
chr1  554315  554316  6
chr1  554316  554317  5

```

---
Tags: sequence-alignment, bedtools, coverage
---


# Question
Title: multi-sequence alignment of samples with multiple contigs each
I want to perform a multi-sequence alignment on 12 samples that clustered based on cgMLST. Ultimate goal is to find out whether they differ by the presence of a few genes.

I have performed multi-sequence alignment in the past using MAFFT, the main difference being that, back then, I was aligning protein sequences, and now I want to align multiple **denovo assembled genomes, with multiple contigs each**

I first had a look at MAFFT, the program that I used in the past for protein multi-sequence alignment, though I noticed that it considers each contig in a denovo assembly as independent sequence, which is not what I want

Then I had a look at Mauve, but was discouraged by installation trouble-shooting. Also I read on their webpage:

> When an individual file contains several sequence entries, they will be concatenated and the whole concatenated sequence will be aligned to the sequences in the other files. This behavior allows multi-chromosomal genomes to be aligned by including all chromosomes in a single sequence file. Similarly, incomplete genomes that consist of several sequence contigs can be aligned, though beware that incorrectly ordered contigs will appear as genome rearrangements in the Mauve viewer.

So I guess concatenating all contigs in a single fasta file is not an option, as artefacts (genome rearrangements) might appear

It turns out that MAFFT has also an experimental feature to multi-align fastq reads, which sounds attractive to me. Therefore I would like to ask if anyone can recommend this new/ experimental feature. Also if this option is available on the command-line tool, as the mafft help is quite limited (it does not specify the format of input and output files):

```
  MAFFT v7.505 (2022/Apr/10)
  https://mafft.cbrc.jp/alignment/software/
  MBE 30:772-780 (2013), NAR 30:3059-3066 (2002)
------------------------------------------------------------------------------
High speed:
  % mafft in > out
  % mafft --retree 1 in > out (fast)

High accuracy (for <~200 sequences x <~2,000 aa/nt):
  % mafft --maxiterate 1000 --localpair  in > out (% linsi in > out is also ok)
  % mafft --maxiterate 1000 --genafpair  in > out (% einsi in > out)
  % mafft --maxiterate 1000 --globalpair in > out (% ginsi in > out)

If unsure which option to use:
  % mafft --auto in > out

--op # :         Gap opening penalty, default: 1.53
--ep # :         Offset (works like gap extension penalty), default: 0.0
--maxiterate # : Maximum number of iterative refinement, default: 0
--clustalout :   Output: clustal format, default: fasta
--reorder :      Outorder: aligned, default: input order
--quiet :        Do not report progress
--thread # :     Number of threads (if unsure, --thread -1)
--dash :         Add structural information (Rozewicki et al, submitted)

```

Any other recommendation on another strategy to align multiple denovo assembled genomes with multiple contigs each are welcome

# Answer
MUMmer4 is a versatile alignment tool for DNA and protein sequences. It supports one reference genome and up to 32 query genomes. MUMmer4 will align every contig in each genome to the reference genome.
> 2 votes

# Answer
If you have a reference genome (or are willing to designate one of your de novo assemblies a reference), you may find QUAST helpful.

QUAST will perform an alignment of all genomes against the reference genome and use that as the basis for comparisons. I believe that under the hood it is running minimap2 for this purpose.

It includes a large-scale alignment viewer ("Icarus") to investigate possible differences within contigs, computes high-level statistics on genome completeness/etc, and has other handy features.

You can see the most recent QUAST paper or the manual for more information.
> 2 votes

# Answer
in the end, I just loaded the `.bam` files into artemis and, by inspecting the depth (heat-map), I could check which samples had the genes I was looking for:
> 2 votes

# Answer
It won't work on more than 3 genomes at a time, but SynMap3D from the Comparative Genomics webserver (https://genomevolution.org/coge/) should do the trick. I'm mentioning it here because the advantage, for me (in contrast to, say, Mauve), was user-friendliness.

Small note: Be mindful of what you name your uploaded samples. You can set your sequence not to be publicly available, but if you give it a new name, as opposed to a species name already in the database, that will be visible.
> 1 votes

---
Tags: sequence-alignment, fasta, sam, multiple-sequence-alignment
---


# Question
Title: finding overlapping motifs to increase length of motif
I am able to find matching motif in my sequence, and I would like to now find overlapping motifs. Basically, after matching my motif, I want to find the 6 amino acids after it. This is the code below that I used to find the motif:

```
import Bio
import regex

from Bio import SeqIO
input_file = 'sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     result=regex.finditer(r"[YFWLIMVA]..[LMALVN]..[AGSTCD].[LAIVNFYMW]",sequence)
     for x in result:
     print(name, x.start(), x.end(), x.group())

```

The above code works perfectly becasue it give me the sequence id, positions and the motif. The output is below

```
P1  33 41 VTLLPAADL

```

Right now, what I want to do is to also get the 6 amino acids after matching this motif, such that I get an output like the one below where the length will now be 15,.

```
P1 33 47 VTLLPAADLLMAIID

```

The code that I have tried to get the 6 amino acids after my match is below.

```
import Bio
import regex

from Bio import SeqIO
input_file = 'sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     result=regex.finditer(r"[YFWLIMVA]..[LMALVN]..[AGSTCD].[LAIVNFYMW]",sequence)
     for x in result:
     print(name, x.start(), x.end() + 6, x.group())

```

This the output it gives me:

```
#It does not extend my motif by 6 amino acids, after getting the match, it only extends the
#end position by 6. From 41 to 47.
P1  33 47 VTLLPAADL

#My desired output is this which include the overlapping LMAIID, the LMAIID is after the
#match.
P1   33 47 VTLLPAADLLMAIID

```

I also tried the code below, but it returns an error.

```
import Bio
import regex

from Bio import SeqIO
input_file = 'sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     result=regex.finditer(r"[YFWLIMVA]..[LMALVN]..[AGSTCD].[LAIVNFYMW]",sequence)
     for x in result:
     print(name, x.start(), x.end() + 6, x.group() +6)

```

# Answer
This is how it managed to give the desired output.

```
import Bio
import re

from Bio import SeqIO
input_file = 'sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     matches=re.finditer(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]",sequence)
     for x in matches:
     print(name, x.start(), x.end()+6,sequence[x.start():x.end()+6])

#Output
P1 33 47 VTLLPAADLLMAIID

```
> 1 votes

# Answer
Just add 6 extra aminoacids to your motif:

```
matches=re.finditer(r"[YFWLIMVA]..[LMALVN].[AGSTCD]..[LAIVNFYMW]......",sequence)

```

That's all you need, now the existing code will give you the answer you want.
> 0 votes

---
Tags: proteins, python, motifs, multi-fasta
---


# Question
Title: How to identify mutations in a viral genome
I have a fasta file with multiple sequences comprising 38 sequences. The length of the sequences are around 11000 bp.

How can i get changes in the genomes based in a reference genome? (aa subtitutions or deletions)

# Answer
> 1 votes
Python3 via the `difflib` library/package. `AlignIO` can be leveraged (Biopython) to construct the `align` list. To be honest if your coding isn't strong you could just cut and paste the alignment into this code. Nucleotides not amino acids ... same idea just easier to type. It will work for aa and gaps. "Ref" (output) refers to the first alignment position only.

```
import difflib

align=['AGCTAGATGATATGA', 'AGCTAGATGATATGA', 'AGCTAGGTGATATGA', 'AGCTAGATGATATGA', 'AGCTAGATGATATGA', 'AGCTAGAAGATATGA', 'AGCTAGAGGATATGA', 'AGCTAGATCATATGT']
ref= 'AGCTAGATGATATGA'
for query in align:
    print('{} (ref) => {} (query)'.format(ref,query))
    for loc,pos in enumerate(difflib.ndiff(ref, query)):
        if pos[0]==' ':
            continue
        elif pos[0]=='-':
            print(u'Ref "{}"'.format(pos[-1]))
        elif pos[0]=='+':
            print(u'Query has "{}" at position {}'.format(pos[-1],loc))
    print ('')

```

Output,

```
AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)

AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)

AGCTAGATGATATGA (ref) => AGCTAGGTGATATGA (query)
Ref "A"
Query has "G" at position 7

AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)

AGCTAGATGATATGA (ref) => AGCTAGATGATATGA (query)

AGCTAGATGATATGA (ref) => AGCTAGAAGATATGA (query)
Ref "T"
Query has "A" at position 8

AGCTAGATGATATGA (ref) => AGCTAGAGGATATGA (query)
Ref "T"
Query has "G" at position 8

AGCTAGATGATATGA (ref) => AGCTAGATCATATGA (query)
Ref "G"
Query has "C" at position 9
Query has "T" at position 16

```

Sequences must be aligned, i.e. ref sequence within the global alignment and separated. Hope that makes sense. It will work with multiple mutations between ref and query as the last example points out (two mutations per example query was just too much to type).

A reproducible example in future please.

---
Tags: phylogenetics, genome, phylogeny, reference-genome, mutations
---


# Question
Title: What options should I explore to improve the output of the protein reconstruction?
We developed a neural network-based protein reconstruction tool to reconstruct the main chain from only CA atoms.

1. we generated data from some selected PDBs from the RCSB website to train an NN model.
2. we then use those data to train the NN model.
3. we select some test PDBs, strip all atoms except CA atoms, and save them in files.
4. we pass those CA-only PDBs through the NN model, obtain a reconstructed main chain, and save them in files.
5. we compare original and reconstructed PDB files and calculate CRMSD values.

If we obtain a CRMSD value of 0.3559 and want to improve it what options should we explore?

* Using Keras.
* Not CNN nor RNN. It's a straightforward MLP (4 layers in total). No feedback, no convolution.
* training set is actually large: 1398438 rows and 102 columns in the data file.

# Answer
The first step is to assess the accuracy, recall and precision of the model then look at confusion matrix. Thereafter consider ROC AUC and if all those are fine ... it's biological. If the model is not 'accurate' where accuracy \>0.85 (minimum) its about training the model and that is a different question.

My guess (*edit: which from the comments was wrong*) is:

* the sample size is insufficient for a presumed deep learning model.
* this will generate problems in the accuracy of the model.

The solution to low sample is called 'data augmentation', a lot of biological knowledge is needed to successfully augment a data set. In this case I suspect the level of biological understanding is strong, but the understanding of deep learning is less advanced.

In summary, low accuracy is a feature of small sample sizes within a deep learning model. To give some idea deep learning needs 500 000 sample size to be reasonable.

---

**Update from the comments**, the data set is large:

* a sample of 1398438 rows and 102 columns (data points).

This is a good sample for deep learning and low sample size is not the cause.

The model is `keras` via a 4-layered MLP, so thats fine, but not clear whether its CNN (probably) or RNN.

The likely issue is under-fitting and this is diagnosed via the above statistics for example accuracy and ROC AUC and solution is to enhance the MLP. If there is over-fitting occurring thats a different question.
> 1 votes

---
Tags: protein-structure, machine-learning, public-dataset, deep-learning
---


# Question
Title: How to resolve Error in fix.by(by.x, x) : 'by' must specify a uniquely valid column while using merge function in R
I have two csv formatted data frames that I want to merge. First data frame contains 9 columns and the second contains three columns. Based on column names: start, end, function I am trying to merge the two data frames. Following is what I tried:

```
df1<-read.csv(file = "LuC_A_TE_gff.csv")
df2<-read.csv(file = "LuC_A_TE_1kb_range_ids.csv")
merged.data<-merge(df1, df2, by.x=c("start","end","function"), by.y=c("start","end","function"))

```

> Error in fix.by(by.x, x) : 'by' must specify a uniquely valid column.

```
df1

```

Output

```
chr,RM,sim,start,end,stat,strand,id,function
Chr01A,RepeatMasker,similarity,1,1122,4438,+,.,Target=rnd-4_family-4692 2 652;Class=LTR;PercDiv=4.9;PercDel=9.9;PercIns=5.4;
Chr01A,RepeatMasker,similarity,1123,2001,6668,+,.,Target=25:218931..219866_LTR 171 936;Class=LTR/Gypsy;PercDiv=6.5;PercDel=6.5;PercIns=0.7;
Chr01A,RepeatMasker,similarity,1998,2421,3164,+,.,Target=rnd-4_family-4692 102 559;Class=LTR;PercDiv=6.1;PercDel=8.0;PercIns=0.0;
Chr01A,RepeatMasker,similarity,2002,2119,1182,+,.,Target=rnd-5_family-5429 198 324;Class=DNA/Tc1-Mariner;PercDiv=5.1;PercDel=5.7;PercIns=0.0;
Chr01A,RepeatMasker,similarity,2284,3113,6265,+,.,Target=25:218931..219866_LTR 32 908;Class=LTR/Gypsy;PercDiv=4.6;PercDel=6.4;PercIns=0.7;
Chr01A,RepeatMasker,similarity,3003,3103,1238,+,.,Target=rnd-5_family-5429 159 202;Class=DNA/Tc1-Mariner;PercDiv=9.1;PercDel=5.5;PercIns=1.6;

```

.

```
df2

```

output

```
start,end,function
10000060,10000222,Target=rnd-6_family-1538 599 686;Class=DNA/CMC-EnSpm;PercDiv=24.8;PercDel=1.1;PercIns=3.2;
10000190,10000200,Target=33:34356884..34364107_INT-int 202 270;Class=LTR;PercDiv=20.5;PercDel=5.9;PercIns=1.6;
10000201,10000302,Target=rnd-6_family-3392 448 555;Class=DNA;PercDiv=10.9;PercDel=6.9;PercIns=0.9;
10000223,10000356,Target=rnd-6_family-1538 196 343;Class=DNA/CMC-EnSpm;PercDiv=20.3;PercDel=11.2;PercIns=0.7;
10000312,10000423,Target=fa.05_5_22233 5 115;Class=DNA;PercDiv=3.6;PercDel=0.0;PercIns=0.9;
10001256,10001332,Target=rnd-4_family-2291 200 273;Class=DNA;PercDiv=17.6;PercDel=0.0;PercIns=4.0;
10001590,10001909,Target=rnd-4_family-274 376 722;Class=DNA/CMC-EnSpm;PercDiv=12.3;PercDel=9.4;PercIns=0.9;
10001650,10001757,Target=rnd-6_family-7385 84 183;Class=DNA;PercDiv=16.3;PercDel=0.9;PercIns=9.0;
10002071,10002421,Target=rnd-3_family-11 116 430;Class=DNA/hAT-Ac;PercDiv=4.2;PercDel=0.0;PercIns=11.4;
10002933,10003171,Target=rnd-5_family-572 47 262;Class=DNA;PercDiv=17.4;PercDel=4.6;PercIns=15.7;
10003024,10003085,Target=Chr23:11562195..11564516 1491 1552;Class=RC/Helitron;PercDiv=4.8;PercDel=0.0;PercIns=0.0;
10003030,10003207,Target=rnd-4_family-825 65 227;Class=DNA;PercDiv=12.9;PercDel=7.3;PercIns=17.2;
10003132,10003142,Target=rnd-5_family-1284 903 950;Class=DNA/Kolobok;PercDiv=16.7;PercDel=7.7;PercIns=2.4;
10003143,10003941,Target=fa.24_3_1586 1 753;Class=DNA;PercDiv=17.3;PercDel=4.8;PercIns=11.2;
10003209,10003451,Target=rnd-5_family-502 491 715;Class=DNA/IS3EU;PercDiv=19.6;PercDel=7.4;PercIns=16.0;
10003271,10003496,Target=rnd-6_family-3401 557 818;Class=DNA/IS3EU;PercDiv=22.4;PercDel=19.5;PercIns=3.0;
10003414,10003567,Target=rnd-1_family-1158 5 166;Class=Unknown;PercDiv=25.1;PercDel=9.7;PercIns=4.3;
10003494,10003560,Target=rnd-4_family-825 2 63;Class=DNA;PercDiv=13.7;PercDel=4.5;PercIns=12.9;
10003631,10003854,Target=rnd-1_family-378 1 406;Class=DNA/P;PercDiv=4.2;PercDel=0.0;PercIns=0.1;
10004081,10004200,Target=fa.00_5_3944 125 245;Class=DNA;PercDiv=10.8;PercDel=0.8;PercIns=0.0;
10004147,10004282,Target=fa.13_3_23430 1409 1547;Class=DNA/hAT;PercDiv=6.7;PercDel=2.9;PercIns=0.7;
10004172,10004199,Target=fa.02_2_13330 434 524;Class=DNA;PercDiv=16.6;PercDel=3.1;PercIns=7.6;
10004200,10004427,Target=rnd-1_family-377 2 229;Class=DNA;PercDiv=11.3;PercDel=3.1;PercIns=3.1;
10004279,10004320,Target=fa.13_3_23430 882 923;Class=DNA/hAT;PercDiv=6.8;PercDel=0.0;PercIns=0.0;
10004321,10004656,Target=fa.15_1_6197 348 728;Class=DNA/Tc1-Mariner;PercDiv=13.1;PercDel=13.4;PercIns=0.0;
10004428,10004480,Target=fa.02_2_13330 525 709;Class=DNA;PercDiv=16.6;PercDel=3.1;PercIns=7.6;
10004482,10004804,Target=rnd-1_family-77 1 331;Class=SINE/tRNA-Core-RTE;PercDiv=7.5;PercDel=4.0;PercIns=1.5;
10004515,10004566,Target=rnd-6_family-2042 304 370;Class=DNA/TcMar-Tc1;PercDiv=7.7;PercDel=28.9;PercIns=0.0;
10004658,10005425,Target=rnd-5_family-3762 370 1163;Class=DNA/hAT-Ac;PercDiv=7.8;PercDel=3.8;PercIns=0.4;

```

# Answer
When I loaded these files into R, I got the following result, which matches your error (except that it's been translated into German):

```
> df1 <- read.csv("df1.csv")
> df2 <- read.csv("df2.csv")
> merge(df1, df2, by=c("start","end","function"))
Fehler in fix.by(by.x, x) : 'by' muss gültige Spalte(n) spezifizieren

```

This error is obscure, and it took me a while to get my head around what was going on. The problem is that 'function' is a reserved word in R, so it renamed the `function` column to add a period \[`.`\] at the end of it:

```
> colnames(df1)
[1] "chr"       "RM"        "sim"       "start"     "end"       "stat"      "strand"
[8] "id"        "function."
> colnames(df2)
[1] "start"     "end"       "function."

```

So R was looking for a column name called "function", didn't find one, and got confused.

I think the easiest fix to this is to use the new names that R has given the tables:

```
> merge(df1, df2, by=c("start","end","function."))
[1] start     end       function. chr       RM        sim       stat      strand    id
<0 Zeilen> (oder row.names mit Länge 0)

```

\[note: this has no output, because the examples you provided have no common start/end/function\]

Alternatively, let R work out the common fields itself, as its default for `merge` is to intersect the column names (which is what you want in this case):

```
> merge(df1, df2)
[1] start     end       function. chr       RM        sim       stat      strand   id
<0 Zeilen> (oder row.names mit Länge 0)

```

Which, for completeness, is equivalent to this:

```
> merge(df1, df2, by=intersect(names(df1), names(df2)))
[1] start     end       function. chr       RM        sim       stat      strand    id
<0 Zeilen> (oder row.names mit Länge 0)

```
> 1 votes

---
Tags: r, merge
---


# Question
Title: finding the frequency of a motif-matching peptide
*This question has also been asked on Biostars*

I want to find the frequency of each motif-matching peptide. The sequences are:

```
> Sequence_1
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL

>Sequence_2
MAEVRKFTKRLSKPGTAAELRQSVSEAVRGSVVLEKAKLVEPLDYENVITQRKTQIYSDP
LRDLLMFPMEDISISVIGRQRRTVQSTVPEDAEKRAQSLFVKECIKTYSTDWHVVNYKYE
DFSGDFRMLPCKSLRPEKIPNHVFEIDEDCEKDEDSSSLCSQKGGVIKQGWLHKANVNST
ITVTMKVFKRRYFYLTQLPDGSYILNSYKDEKNSKESKGCIYLDACIDVVQCPKMRRHAF
ELKMLDKYSHYLAAETEQEMEEWLIMLKKIIQINTDSLVQEKKDTVEAIQEEETSSQGKA
ENIMASLERSMHPELMKYGRETEQLNKLSRGDGRQNLFSFDSEVQRLDFSGIEPDVKPFE
EKCNKRFMVNCHDLTFNILGHIGDNAKGPPTNVEPFFINLALFDVKNNCKISADFHVDLN
PPSVREMLWGTSTQLSNDGNAKGFSPESLIHGIAESQLCYIKQGIFSVTNPHPEIFLVVR

>Sequence_3
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

> Sequence_4
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL

>Sequence_5
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

```

The fasta file actually has more sequences, I just chose the top 5, to put on here.

The code used is below:

```
import Bio
import regex

from Bio import SeqIO
input_file = 'Sequences.fasta'
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
     name, sequence = fasta.id, str(fasta.seq)
     matches=regex.finditer(r"(P[A-Z]{2}P..)",sequence)
     for m in matches:
        print(name, m.start(), m.end(), m.group())

```

This is giving me my desired output, which is below:

```
Sequence1 73 79 PFLPGK
Sequence1 281 287 PLPPKL
Sequence1 288 294 PSSPAH
Sequence4 73 79 PFLPGK
Sequence4 281 287 PLPPKL
Sequence4 288 294 PSSPAH

```

Now I want to:

1. Calculate the number of sequences in that file. I know in linux environment I would have to use `grep -c ">" Sequences.fasta` to get the total number of sequences in the fasta file, but how do I do it in python?.
2. Once I know how many sequences are in file in total, I want to count how many of those sequences have the motif-matching peptides above. For example how many sequences have PFLPGK at positions 73 79? In the example above PFLPGK its in sequence 1 & 4, so the **frequency** is 2, meaning that two sequences have PFLPGK. How do I find that frequency for each of the peptides using python.

Many thanks for your assistance.

# Answer
> 0 votes
For the first part is easy,

```
from Bio import SeqIO
seqs = SeqIO.to_dict(SeqIO.parse("myfile.fa", "fasta"))
print (len(seqs))

```

The second part is solvable Python's amazing slicing, this is the concept ....

```
from Bio import SeqIO
loci = {'PFLPGK':[1,10], 'PFLPGK':[50:60] }
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
mydict = {}
for key, values in loci.items():
    myrange = '-'.join(values)
    mydict = {key:myrange} # you need to initiate the dictionary without defaultdict
    for record in fasta_sequences:
        if record.seq[values[0]-1:values[1]-1] == key:
             mydict[key][myrange].append(record.id)

```

To print it:

```
key, value in mydict.items():
    print("No. of" + key + '=' + len(mydict[key][value]))

```

**Note 1**, the second answer is concept only.

**Note 2** there can be an issue between the `list` location and the sequence location ... that is **super** important otherwise it will not work. I've added `-1` from both locations of the amino acid and the location on the `list`, *but this needs checking*. This is because `list`s start at `0` whilst amino acids start at `1`

---

**Update** *Advice on Python* requested by OP

1. Answer to Q1 has been debugged, tested and works.
2. Use Q2 as the basis ... just debug it (see point 4.), it will be a bit fiddly because moving from `lists` to aa strings.
3. Use an IDE, I use Spyder: don't get too caught up in "the best IDE" they just make life easier.
4. Book: Mark Lutz "Learning Python 5th Edition" by O'Reilly section on `lists` (slicing) and `dictionaries` *****.
5. Online resources `GeekforGeeks` on `lists` (slicing) and `dictionaries`.
6. BioPython docs on SeqIO search for `record.id` and `record.seq` ... you need to understand the SeqIO object.

***** This is physically a huge book so an electronic version is strongly preferred just for logistics. This goes beyond the full modern implementation of `OOP` in bioinformatics (which is not very modern in truth and can become unreadable).

---
Tags: phylogenetics, proteins, python, motifs, homework
---


# Question
Title: How can I find the chromosomal location of a list of genes?
I have a list of genes nearly 20000: gene name (column1) and coordinates (columns 2 and 3)

```
head genes
    Ma0_73Ns    2359    2923
    Ma0_73Ns    4727    4868

```

I want to find the chromosomal location of each gene, on which chromosome each gene located? my study organism is *Helianthus annuus* (sunflower)

# Answer
Helianthus annuus is an Ensembl plants genome. You should be able to use the bioMart there to pull out genomic co-ordinates (chromosomes etc) for genes. However, searching for the gene names you posted does not return any results. The same is true of the two websites posted above. Are they real gene names or just examples you made up for the post?
> 3 votes

# Answer
Ok so I found two websites that could help. The one I found first is more intuitive to go through, but it looks like it could be outdated as the second website I found seems more recent so it probably has a more complete genome to use as the reference. (Sorry, I haven't worked personally with sunflowers and don't know for sure where you would go for the "official" reference genome, or if there even is an "official" one, but hopefully one of these websites will help.)

1) This is a website I found that looks like it is a website for the sunflower genome and includes gene annotation information. (The recent publications list only goes up till 2015 though so it may not be the most recent data)

https://www.sunflowergenome.org/annotations-data/

If you scroll down to the bottom of this page there is a text file detailing the descriptions of how the files are formatted.

For example, the linked to file gives the following:

> * Ha412v1r1\_CDS\_v1.0.fasta.gz
>
> + Description: Nucleotide CDS for each gene
>     + Header format: \>TranscriptID GeneID Chromosome\_start\_end\_strand
>     + Example: \>Ha1\_00044280-RA Ha1\_00044280 Ha1\_94590914-94591860\_-

If you download this fasta file you could then create a python script that will read in the file and extract the headers. If each header is a string variable called "mytext", you could then extract the TranscriptID, GeneID, and chromosomal info and using mytext.split(" "). You could then split on "\_" to further parse the chromosomal location info. You could save the extracted information in a dictionary using the gene name as a key and parse through the original file you mentioned containing the list of genes, in order to get the chromosomal location information you desire. The one caveat is that this is assuming that your gene names are the same as the geneIDs the above website uses.

I did not exhaustively examine the list of files and perhaps there is another one that would be a better one to parse through for your purpose, but the linked to text file gives more details on the files and their formats and presumably they could be parsed in a similar manner.

2) Double checking my first answer I found this article:

> Badouin H, Gouzy J, Grassa CJ, Murat F, Staton SE, Cottret L, Lelandais-Brière C, Owens GL, Carrère S, Mayjonade B, Legrand L. The sunflower genome provides insights into oil metabolism, flowering and Asterid evolution. Nature. 2017 Jun;546(7656):148.

Which links to a website as follows:

> "The integrative web interface Heliagene provides visualization, querying tools for data mining and network exploration for the community (https://www.heliagene.org)."

So it looks like this website also has reference genomes for sunflowers, although a lot of the annotation information seems to be locked unless you have a login. Like I mentioned previously, this website is a little less intuitive to go through but you might be able to find something helpful here if you think the other website isn't official or recent enough.
> 2 votes

# Answer
If you can add a column of RefSeq identifiers for each gene to your table, first download a whole-genome GFF file for your organism:

```
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/4232/101/GCF_002127325.2_HanXRQr2.0-SUNRISE/GCF_002127325.2_HanXRQr2.0-SUNRISE_genomic.gff.gz

```

and then open the file in Python using a package like gffutils and extract all the rows which contain a RefSeq identifier located in table (the RefSeq identifier is the first column in the gff file linked in the wget command above). For all the rows you extract, which chromosome it belongs to will be stated in the Attributes column, which you can also extract using gffutils.
> 1 votes

---
Tags: gene, sequence-annotation
---


title: How to wget a gff from NCBI database?

text:
# Question
Title: How to wget a gff from NCBI database?
Say I want to use wget to download a gff file for the X chromosome of release 5.57 of the Drosophila melanogaster genome. The NCBI page with the gff file for just this is here.

However, I do not see any way to directly wget the gff from this page. What I can do is click "Send to" at the top, select File, and change the Format dropdown to GFF3, and finally click 'Create File', and the gff file downloads onto my computer.

But then I lose the ability to use the -a flag to log the download details, I can't download the gff using wget from a remote server directly but instead I need to download it locally and then upload it, etc. Is there any way to use wget to download gff files from NCBI programmatically instead of having to download it manually from the website?

# Answer
> 3 votes
my suggestion would be to go via the assembly page. For your example link, I would

1. Click the link to the BioProject: PRJNA164
2. On the right side, under "Related information" choose the assembly
3. On the assembly page, under "Access the data" on the right side, I would choose "FTP directory for RefSeq assembly"

Once in the RefSeq folder, copy the link to the `*_genomic.gff.gz` and use wget.

```
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/215/GCF_000001215.4_Release_6_plus_ISO1_MT/GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.gff.gz

```

---

**Edit**: The initial question requested an older versions, which can be found in RefSeq's archive of prior releases. For this specific Drosophila assembly chromosome X, the version 5 can be found at https://ftp.ncbi.nlm.nih.gov/genomes/archive/old\_refseq/Drosophila\_melanogaster/RELEASE\_5/CHR\_X/NC\_004354.gff

# Answer
> 1 votes
I believe I figured out a way to do this for my particular example (although in general cases Adrian's answer works). To get release 5.57 (an earlier release) of chromosome X of Drosophila melanogaster, I believe you can go to the Datasets Genome page of D. melanogaster, scroll down to the Revision History and select Release 5, and at the top of the page, NCBI lets you click on a 'datasets' or 'curl' box to copy the exact command line argument needed to download a GFF (and other files) of that whole genome. Then, I believe it should be possible to open the gff in Python using a package like gffutils, and filter only for chromosome X rows, and then write those rows to a new gff file.

---
Tags: phylogenetics, bash, data-download, gff3
---


# Question
Title: Understanding regularization step in Meta-cell pre processing
I'm reading through this paper: Baran, Y., Bercovich, A., Sebe-Pedros, A., Lubling, Y., Giladi, A., Chomsky, E., ... & Tanay, A. (2019). MetaCell: analysis of single-cell RNA-seq data using K-nn graph partitions. Genome biology, 20(1), 1-19.

Maybe I am not familiar with this notation but I don't understand the regularization procedure as follows:

1. What is happening in these max functions? What does the 0 indicate? They have some terms $x$ and then in the max function they write $max( x, 0)$.
2. I assume $K$ is some regularization parameter, like a cost. Is $s\_{ij} * s\_{ji}$ the rank product? I.e. the geometric mean?

# Answer
> 2 votes
Without looking in too much detail:

* We may suppose that $S$ is a matrix with dimensions indexed as $(i, j)$. So $s\_{ij}$ is the cell of $S$ indexed to the $i$th row and $j$th column, and $s\_{ji}$ is the cell vice versa. Therefore, symmetrizing $S$ is a matter of ensuring that $s\_{ij} == s\_{ji}$, which you can achieve by multiplying the (ranked?) entries. Symmetric matrices are nice for doing certain kinds of math, and yield a handy representation of undirected graphs, which the authors appear to be using for kNN.
* $max(x, 0)$ is just saying that you want either $x$ if $x$ is non-negative, or $0$ if it is negative. Thus, you make the whole matrix non-negative, which makes some handy matrix math possible. You can see how this works computationally by inspecting e.g. the Python `max()` function.

---
Tags: scrnaseq, statistics, graphs
---


# Question
Title: Why do we need to find minimum energy in a protein chain?
The above research paper is about a software tool for reconstructing a protein's main chain model only from its alpha-carbon backbone (aka, C-alpha trace).

This tool provides two options: conversion (1) with energy minimization, and (2) without energy minimization.

The above paper says that:

> ...
> The goal of energy Minimization is to find a set of coordinates representing the minimum energy conformation for the given structure.
> ...
> In molecular modeling we tend to are particularly curious about minimum points on the energy surface. Minimum energy arrangements of the atoms correspond to stable states of the system; any movement off from a minimum provides a configuration with a better energy.
> ...

I haven't understood the quote from the second paper.

Why do we need to know the minimum energy conformation of a protein chain? I.e., why do we want to do *energy minimization*?

# Answer
> 6 votes
After translation, proteins fold to achieve a specific 3D structure which allows them to perform their specific functions. The structure of a protein computationally inferred through energy minimization often corresponds quite well to that specific, real-life folded 3D structure.

If you are interested in the exact reason why proteins tend to fold to a conformation of minimal energy (even though this is an oversimplification of what's really happening) you would have to dive a bit deeper into the thermodynamics of proteins: https://en.wikipedia.org/wiki/Protein\_dynamics

**TL;DR: Energy minimization is often used to get an idea of the protein's real-life conformation (i.e. structure) through computational means.**

---
Tags: proteins, protein-structure, protein-protein-interaction
---


# Question
Title: Multiple values for memory resources shown when Snakemake is run
I am running a Snakemake (V7.14.0) pipeline on the cluster and provided resources in a resource config file like:

```
threads:
    tool1: 4
mem_gb:
    tool1: 12

```

Point to it in my rules like:

```
resources:
        mem_gb=config['mem_gb']['tool1']
threads:
        config['threads']['tool1']

```

And the snakemake cluster command like:

```
snakemake --configfiles \
    config/config_resources.yml \
    config/config_parameters.yml \
    --latency-wait 60 \
    --use-conda \
    --use-singularity \
    --printshellcmds \
    --jobname MyPipe{name}_{jobid} \
    --jobs 300 \
    --cluster \
    "bsub -W 200 -n {threads} -M {resources.mem_gb}G \
    -R 'rusage[mem={resources.mem_gb}G]' \
    -o run_logs/{name}_%J.out \
    -e run_logs/{name}_%J.err"

```

What I do not understand when I run the pipeline is that it shows this in the terminal:

> threads: 4
>
> resources: mem\_mb=19442, disk\_mb=19442, tmpdir=/tmp, mem\_gb=12

So next to the one I provided (mem\_gb), there are: mem\_mb, and disk\_mb. They are clearly not the same value. In the cluster job it looks as though the MEMLIMIT is indeed the one I give, but the MEM usage is different (higher). Is it because Snakemake does not understand the mem\_gb and I specifically need to use that word?

Thanks in advance

# Answer
> 1 votes
You have defined `mem_gb` but not `mem_mb` in your resource configuration file. You can set them to same value by adding something like `mem_mb=config['mem_gb']['tool1']`. Additional information in their documentation.

---
Tags: snakemake, hpc
---


title: file path of BiocManager:install()

text:
# Question
Title: file path of BiocManager:install()
I was trying to format my GWAS summary statistics with the `MungeSumstats` R package and would need to install `SNPlocs.Hsapiens.dbSNP155.GRCh37` and `BSgenome.Hsapiens.1000genomes.hs37d5` as reference panel.
I typed `BiocManager::install("SNPlocs.Hsapiens.dbSNP144.GRCh37") BiocManager::install("BSgenome.Hsapiens.1000genomes.hs37d5")` but failed due to Internet connection problems. So I planned to download the files from the web browser and copy them to the R library. I was wondering where I should put the files. If they were downloaded from R, where would they be saved?
I would really appreciate any answers, thanks!

# Answer
**First, I would try again to download via the bioconductor command. If your internet is intermittent it may just take a few tries.**

If that doesn't work, I would then strongly recommend using `BioCManager` to install from a local source wherever you've downloaded the files. If you try to copy the files around there are a lot of ways that things can go wrong, and there are a lot of ways that it still might not work even if you put it in the right place, if you don't use the package manager as intended.
> 3 votes

---
Tags: r, bioconductor, reference-genome
---


# Question
Title: How can I extract the Protein distance map data in Julia?
I calculated the distance between C-alphas using Julia, but there is no way to extract distance map data from Julia. I don't even know the data file form, png or csv or etc. I don't want to extract it as an image, I just want all raw data.

# Answer
> 1 votes
for someone who have same question as mine... this helps you to extract distance map data file.

```
julia> using DelimitedFiles

julia> writedlm("calphadistances.out",dm.data,"")

```

---
Tags: proteins, protein-protein-interaction
---


# Question
Title: Bcftools 1.16 able to add F_MISSING tag?
I tried adding the F\_MISSING tag using bcftools 1.16. When I run this command:

```
bcftools +fill-tags input.vcf.gz -- -t 'F_MISSING' | bcftools view -i 'INFO/F_MISSING<0.25' -Oz -o output.vcf.gz

```

I get the following error:

```
Error parsing "--tags F_MISSING": the tag "F_MISSING" is not supported

```

This command runs fine using bcftools 1.15. However, version 1.15 gives complications with other packages I use in my snakefile. Do you maybe know alternatives for how to add F\_MISSING using bcftls 1.16 (if possible)?

When I type `bcftools +fill-tags --version`:

bcftools 1.9 using htslib 1.9

plugin at 1.9 using htslib 1.9

##SOLUTION##

Indeed the issue was that I was not installing the most recent version of Conda. I solved it by changing the .condarc file to **solelely** include the following lines:

```
channels:
    - conda-forge
    - bioconda
    - defaults

```

The order is crucial as well.

# Answer
The issue was that not installing the most recent version of Conda. It was solved by changing the .condarc file to **solely** include the following lines:

```
channels:
    - conda-forge
    - bioconda
    - defaults

```

The order is crucial as well.

*Converted to answer from question.*
> 0 votes

---
Tags: bcftools
---


# Question
Title: Can the information regarding atomic bonds be interpreted only from atom lines in a PDB file?
The following are the first 14 atoms from the protein **1CRN**:

```
ATOM      1  N   THR A   1      17.047  14.099   3.625  1.00 13.79           N
ATOM      2  CA  THR A   1      16.967  12.784   4.338  1.00 10.80           C
ATOM      3  C   THR A   1      15.685  12.755   5.133  1.00  9.19           C
ATOM      4  O   THR A   1      15.268  13.825   5.594  1.00  9.85           O
ATOM      5  CB  THR A   1      18.170  12.703   5.337  1.00 13.02           C
ATOM      6  OG1 THR A   1      19.334  12.829   4.463  1.00 15.06           O
ATOM      7  CG2 THR A   1      18.150  11.546   6.304  1.00 14.23           C

ATOM      8  N   THR A   2      15.115  11.555   5.265  1.00  7.81           N
ATOM      9  CA  THR A   2      13.856  11.469   6.066  1.00  8.31           C
ATOM     10  C   THR A   2      14.164  10.785   7.379  1.00  5.80           C
ATOM     11  O   THR A   2      14.993   9.862   7.443  1.00  6.94           O
ATOM     12  CB  THR A   2      12.732  10.711   5.261  1.00 10.32           C
ATOM     13  OG1 THR A   2      13.308   9.439   4.926  1.00 12.81           O
ATOM     14  CG2 THR A   2      12.484  11.442   3.895  1.00 11.90           C

```

These fourteen atoms are members of two **THR** residues.

Can the information regarding atomic bonds be interpreted only from these fourteen lines?

For example,

* the first atom, N, has three bonds attached to CA.
* the second atom, CA, has four bonds; three are attached to N, and one is to C.
* the third atom, C, has four bonds; one is attached to CA, one to oxygen, and two are empty.
* the fourth atom, O, has two bonds; one is attached to the third carbon, and one is to the fifth atom beta-carbon.
* the fifth atom, beta-carbon, has four bonds; one is attached to the fourth atom, two are attached to the first branch (OG1), and one is attached to the second branch (CG2).

Or do I have to use my background knowledge of amino acids to find bonds?

**Note:** I am not talking about hydrogen bonds. I am talking about atomic bonds only in the atoms in the residue.

# Answer
> 2 votes
The official list of bonds in each residue in the PDB is in the Chemical Component Dictionary.

You could also use interatomic distances to infer connectivity, but that's less reliable.

# Answer
> 1 votes
You can use topology files to extract bond information. Here is an example at CHARMM36

---
Tags: proteins
---


# Question
Title: Why am I not being able to detect all the bonds in a residue?
Let us consider the first seven atoms of the protein **`1CRN.pdb`**:

```
ATOM      1  N   THR A   1      17.047  14.099   3.625  1.00 13.79           N
ATOM      2  CA  THR A   1      16.967  12.784   4.338  1.00 10.80           C
ATOM      3  C   THR A   1      15.685  12.755   5.133  1.00  9.19           C
ATOM      4  O   THR A   1      15.268  13.825   5.594  1.00  9.85           O
ATOM      5  CB  THR A   1      18.170  12.703   5.337  1.00 13.02           C
ATOM      6  OG1 THR A   1      19.334  12.829   4.463  1.00 15.06           O
ATOM      7  CG2 THR A   1      18.150  11.546   6.304  1.00 14.23           C

```

This represents an amino acid called ***Threonine***.

The general bond diagram for Threonine is as follows:

So, it has at least six bonds (I am actually not sure!) which are also evident from the seven `ATOM` lines.

Each atom has its radius. In general, if the distance between two atoms is less than the sum of the radii, there is a bond between these two atoms.

Based on the above information, I tried to identify the bond information, and I could detect only two:

```
     atom1           atom2        distance       atm1 Radi      atom2 Radi        sum radi        is bond?
         N              CA            1.50            0.65            0.70            1.35           False
         N               C            2.44            0.65            0.70            1.35           False
         N               O            2.67            0.65            0.60            1.25           False
         N              CB            2.48            0.65            0.70            1.35           False
         N             OG1            2.75            0.65            0.60            1.25           False
         N             CG2            3.86            0.65            0.60            1.25           False
        CA               N            1.50            0.70            0.65            1.35           False
        CA               C            1.51            0.70            0.70            1.40           False
        CA               O            2.36            0.70            0.60            1.30           False
        CA              CB            1.57            0.70            0.70            1.40           False
        CA             OG1            2.37            0.70            0.60            1.30           False
        CA             CG2            2.61            0.70            0.60            1.30           False
         C               N            2.44            0.70            0.65            1.35           False
         C              CA            1.51            0.70            0.70            1.40           False
         C               O            1.24            0.70            0.60            1.30            True
         C              CB            2.49            0.70            0.70            1.40           False
         C             OG1            3.71            0.70            0.60            1.30           False
         C             CG2            2.98            0.70            0.60            1.30           False
         O               N            2.67            0.60            0.65            1.25           False
         O              CA            2.36            0.60            0.70            1.30           False
         O               C            1.24            0.60            0.70            1.30            True
         O              CB            3.12            0.60            0.70            1.30           False
         O             OG1            4.34            0.60            0.60            1.20           False
         O             CG2            3.74            0.60            0.60            1.20           False
        CB               N            2.48            0.70            0.65            1.35           False
        CB              CA            1.57            0.70            0.70            1.40           False
        CB               C            2.49            0.70            0.70            1.40           False
        CB               O            3.12            0.70            0.60            1.30           False
        CB             OG1            1.46            0.70            0.60            1.30           False
        CB             CG2            1.51            0.70            0.60            1.30           False
       OG1               N            2.75            0.60            0.65            1.25           False
       OG1              CA            2.37            0.60            0.70            1.30           False
       OG1               C            3.71            0.60            0.70            1.30           False
       OG1               O            4.34            0.60            0.60            1.20           False
       OG1              CB            1.46            0.60            0.70            1.30           False
       OG1             CG2            2.54            0.60            0.60            1.20           False
       CG2               N            3.86            0.60            0.65            1.25           False
       CG2              CA            2.61            0.60            0.70            1.30           False
       CG2               C            2.98            0.60            0.70            1.30           False
       CG2               O            3.74            0.60            0.60            1.20           False
       CG2              CB            1.51            0.60            0.70            1.30           False
       CG2             OG1            2.54            0.60            0.60            1.20           False
2
Done!

```

How can I detect all bonds?

EDIT: As per pippo1980's request:

```
public class ADTResidue: Dictionary<int, ADTAtom>, IEquatable<ADTResidue>
{
    public string ProteinName { get; set; }
    public string ChainName { get; set; }
    public int ResidueNo { get; set; }
    public Point3d Coordinate { get; private set; }
    public string Code { get; set; }
    public string Full { get; set; }
    public Color ColorRGB { get; set; }

    public ADTResidue(string code, string full, Color colorRgb)
    {
        Coordinate = new Point3d();

        Code = code;
        Full = full;
        ColorRGB = colorRgb;
    }

    private readonly Dictionary<int, CovalentBond> BondDictionary;

    public List<CovalentBond> GetBonds()
    {
        return new List<CovalentBond>(BondDictionary.Values);
    }

    public ADTResidue()
    {
        BondDictionary = new Dictionary<int, CovalentBond>();
    }

    public void CalculateBonds()
    {
        foreach (ADTAtom atom1 in this.Values)
        {
            foreach (ADTAtom atom2 in this.Values)
            {
                if (atom2.Radius <= 0)
                {
                    throw new Exception("zero radius in " + atom2.Name);
                }

                if (atom1 != atom2)
                {
                    double distance = Point3d.Distance(atom1.Coordinate, atom2.Coordinate);

                    bool exists = ABond.IsBondExists(atom1.Coordinate, atom2.Coordinate, atom1.Radius, atom2.Radius);

                    if (exists)
                    {
                        int hashCode = atom1.GetHashCode(atom2);
                        if (!BondDictionary.ContainsKey(hashCode))
                        {
                            BondDictionary.Add(hashCode, new CovalentBond(atom1, atom2));
                        }
                    }
                }
            }
        }

        int count = BondDictionary.Count;

        Console.WriteLine(count);
    }

    #region equality comparison
    // IEquatable<> implementation...
    #endregion
}

public abstract class ABond: Pair<ADTAtom, ADTAtom>
{
    public ABond():base()
    {
    }

    public ABond(ADTAtom first, ADTAtom second) : base(first, second)
    {
    }

    public static bool IsBondExists(Point3d first, Point3d second, double radius1, double radius2)
    {
        double distance = Point3d.Distance(first, second);

        return (distance <= (radius1 + radius2));
    }
}

public class CovalentBond : ABond
{
    public CovalentBond():base()
    {
    }

    public CovalentBond(ADTAtom first, ADTAtom second) : base(first, second)
    {
    }
}

```

# Answer
Your method doesn't work for two possible reasons;

1. The radius information may be incorrect.
2. Your molecule requires some minimization. This option helps to optimize distance between atoms.

You can find some information about the minimization at the following links

1. Wikipedia
2. Gromacs
> 3 votes

---
Tags: proteins, protein-structure, pdb, 3d-structure
---


# Question
Title: Can I perform Expression Analysis on this 10x dataset?
10x Genomics: 20k Human PBMCs is the dataset.

**Description of the dataset:**

> Inputs/Libraries Human peripheral blood mononuclear cells (PBMCs) of a healthy male donor aged 30-35 were obtained by 10x Genomics from AllCells.

> Gene expression and V(D)J libraries were generated from ~33,000 cells (18,470 cells recovered) as described in the Chromium Next GEM Single Cell 5' HT Reagent Kits v2 User Guide (CG000421) using the Chromium and sequenced on an Illumina NovaSeq 6000 to a read depth of approximately 25,000 mean reads per cell for Gene Expression, 15,000 mean reads per cell for TCR Amplified, and 25,000 mean reads per cell for BCR Amplified libraries.

> Paired-end, dual indexing

> Read 1: 26 cycles (16 bp barcode, 10 bp UMI)
> i5 index: 10 cycles (sample index)
> i7 index: 10 cycles (sample index)
> Read 2: 90 cycles (transcript)

This is the RPubs Link for the Single Cell Analysis work I've performed on it.

Can someone confirm that I can do both Single-Cell RNA Seq Analysis *AND* Differential Gene Expression Analysis?

I know that Expression Analysis comes from having a control and a treatment group. I honestly don't know anything about bioinformatics but I'm diving into the fire to get a solid understanding once I've coded everything.

I followed this video and want to follow this other video.

Thank you guys.

# Answer
> 0 votes
Typically what this PMBC RNA-seq data set would be used for is for a control against individuals who would exhibit a pronounced cellular immune response. For example, following a COVID19 infection the PMBCs would be have lots of T-cell CD4+ and CD8+ and are essential for controlling the infection. Thus it is a really good approach to understand the reasons for different clinical outcomes, such chronic infection and severe infections.

I would suspect a cancer patient would show a pronounced shift in cellular immune response, but I'm not certain. Anyway, by comparing a control data set against a COVID19 patient a better understanding of the cellular immune response can be obtained and a differential gene expression analysis would be an important inroad to assess the T-cell up-regulation. DEG would show very pronounced differences between a healthy individual and a patient.

Having a hypothesis about the expected results from DEG is important, because interpreting DEG output is otherwise difficult.

With regards your question, a DEG analysis IMO would be a good analysis and is easy to perform, however there only appears to be a single sample. A DEG analysis in context is control X vs patient Y. The sample you have appears to be a control sample so you would need a second sample. Performing DEG is easy, but do stress that interpreting DEG isn't easy because everything is relative. Thus what is sought usually in PMBC comparisons healthy vs patient is evidence of T-cell CD8+/CD4+ up regulation, any indication about which possibly which HLA type is active (I dunno whether RNAseq will do this - good if it can) and information on the signalling cascade. Basically as much information on the cellular immune response that can be gleaned via RNAseq without performing laborious and difficult immunological experiments requiring specialist apparatus.

There are investigators here with more experience than I in eukaryotic RNAseq analysis, so there input would be helpful.

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


title: Creating a PDB from a figure in paper

text:
# Question
Title: Creating a PDB from a figure in paper
I want to create a PDB file of a drug that is given in a paper. Can somebody share a tutorial or similar resources for a beginner-level user?  Thank you.

# Answer
You can use Image2Smiles tool. It produces smi format, then you can convert smi to pdb using openbabel.

If you just have a few figures, I suggest using PubChem Sketcher.
> 2 votes

# Answer
or you could taket it from :

https://www.rcsb.org/structure/7WMV pdb file

```
HETATM 4865  C13 1YI A 701     106.585 114.218 112.698  1.00 74.13           C
HETATM 4866  C15 1YI A 701     106.050 112.633 114.547  1.00102.13           C
HETATM 4867  C17 1YI A 701     104.024 111.140 110.766  1.00 89.07           C
HETATM 4868  C20 1YI A 701     103.626 108.779 110.527  1.00 83.00           C
HETATM 4869  C21 1YI A 701     104.626 108.507 111.437  1.00 84.86           C
HETATM 4870  C22 1YI A 701     102.234 110.466 109.205  1.00 85.32           C
HETATM 4871  C24 1YI A 701     100.911 110.073 109.846  1.00 66.83           C
HETATM 4872  C26 1YI A 701      92.545 108.407 114.139  1.00 95.68           C
HETATM 4873  C28 1YI A 701     100.141 109.053 109.319  1.00 63.43           C
HETATM 4874  C11 1YI A 701     105.694 113.343 111.806  1.00 75.91           C
HETATM 4875  C12 1YI A 701     105.840 111.929 112.354  1.00 91.76           C
HETATM 4876  C14 1YI A 701     105.990 114.087 114.099  1.00 85.43           C
HETATM 4877  C16 1YI A 701     105.029 110.853 111.663  1.00 99.03           C
HETATM 4878  C18 1YI A 701     105.319 109.548 112.006  1.00100.47           C
HETATM 4879  C19 1YI A 701     103.333 110.096 110.188  1.00 81.96           C
HETATM 4880  C23 1YI A 701     104.801 110.649 116.200  1.00 88.38           C
HETATM 4881  C25 1YI A 701     102.892 107.590 109.921  1.00 99.79           C
HETATM 4882  C27 1YI A 701     100.517 110.745 110.979  1.00 88.80           C
HETATM 4883  C29 1YI A 701      98.556 109.383 111.080  1.00 77.66           C
HETATM 4884  C30 1YI A 701      97.245 109.006 111.756  1.00 83.68           C
HETATM 4885  C31 1YI A 701      97.108 107.503 111.935  1.00 85.27           C
HETATM 4886  C32 1YI A 701      99.328 110.410 111.602  1.00 90.75           C
HETATM 4887  C33 1YI A 701      98.954 108.713 109.937  1.00 76.18           C
HETATM 4888  C34 1YI A 701      96.107 107.127 113.039  1.00 98.03           C
HETATM 4889  C35 1YI A 701      94.780 107.913 113.042  1.00 92.89           C
HETATM 4890  C36 1YI A 701      91.696 107.883 115.298  1.00113.05           C
HETATM 4891  C37 1YI A 701      92.889 109.868 114.425  1.00 87.90           C
HETATM 4892  C38 1YI A 701      91.711 108.277 112.861  1.00 77.85           C
HETATM 4893  N8  1YI A 701      93.804 107.673 114.094  1.00 93.06           N
HETATM 4894  N9  1YI A 701      90.968 106.633 115.149  1.00106.00           N
HETATM 4895  O2  1YI A 701     105.306 111.885 113.648  1.00 99.43           O
HETATM 4896  O3  1YI A 701     106.050 113.444 110.457  1.00 80.66           O
HETATM 4897  O4  1YI A 701     106.560 115.552 112.284  1.00 86.11           O
HETATM 4898  O5  1YI A 701     106.634 114.878 115.050  1.00 99.48           O
HETATM 4899  O6  1YI A 701      94.512 108.703 112.208  1.00 84.41           O
HETATM 4900  O7  1YI A 701      91.625 108.487 116.316  1.00101.67           O
HETATM 4901  S1  1YI A 701     105.267 112.404 116.160  1.00121.62           S

```

picture :

obtained with : https://litemol.org/Viewer/

OOOPSSSSSSSSSS

realized that ligand in pdb its not the full ligand of your picture,

need to get it from Ligand page on RCSB\_PDB : https://www.rcsb.org/ligand/1YI

download links highlighted in pics :

download `/1YI_ideal.sdf` convert via openbabel :

http://www.cheminfo.org/Chemistry/Cheminformatics/FormatConverter/index.html

result

```
WARNING: damped steplength 0.0184262 to 0.00725371
WARNING: damped steplength 0.00508674 to 0.00100138
COMPND    1YI
AUTHOR    GENERATED BY OPEN BABEL 2.3.2
HETATM    1  C   LIG     1       3.537  -2.101   2.235  1.00  0.00           C
HETATM    2  C   LIG     1       2.812  -4.346   1.266  1.00  0.00           C
HETATM    3  C   LIG     1       2.262  -2.902  -2.712  1.00  0.00           C
HETATM    4  C   LIG     1       0.447  -1.602  -3.663  1.00  0.00           C
HETATM    5  C   LIG     1       0.405  -0.869  -2.474  1.00  0.00           C
HETATM    6  C   LIG     1       1.522  -3.533  -5.065  1.00  0.00           C
HETATM    7  C   LIG     1       0.305  -3.796  -5.950  1.00  0.00           C
HETATM    8  C   LIG     1      -5.920  -1.833 -12.368  1.00  0.00           C
HETATM    9  C   LIG     1      -0.581  -4.848  -5.684  1.00  0.00           C
HETATM   10  C   LIG     1       2.908  -1.583   0.938  1.00  0.00           C
HETATM   11  C   LIG     1       3.139  -2.484  -0.313  1.00  0.00           C
HETATM   12  C   LIG     1       3.714  -3.615   2.273  1.00  0.00           C
HETATM   13  C   LIG     1       2.202  -2.172  -1.496  1.00  0.00           C
HETATM   14  C   LIG     1       1.247  -1.156  -1.414  1.00  0.00           C
HETATM   15  C   LIG     1       1.366  -2.666  -3.802  1.00  0.00           C
HETATM   16  C   LIG     1       1.091  -5.111   3.310  1.00  0.00           C
HETATM   17  C   LIG     1      -0.531  -1.162  -4.707  1.00  0.00           C
HETATM   18  C   LIG     1      -0.009  -2.939  -7.020  1.00  0.00           C
HETATM   19  C   LIG     1      -2.131  -4.074  -7.418  1.00  0.00           C
HETATM   20  C   LIG     1      -3.479  -4.118  -8.115  1.00  0.00           C
HETATM   21  C   LIG     1      -3.492  -4.837  -9.473  1.00  0.00           C
HETATM   22  C   LIG     1      -1.224  -3.051  -7.715  1.00  0.00           C
HETATM   23  C   LIG     1      -1.778  -4.991  -6.415  1.00  0.00           C
HETATM   24  C   LIG     1      -4.883  -4.859 -10.129  1.00  0.00           C
HETATM   25  C   LIG     1      -5.503  -3.491 -10.484  1.00  0.00           C
HETATM   26  C   LIG     1      -7.501  -1.753 -12.406  1.00  0.00           C
HETATM   27  C   LIG     1      -5.427  -1.719 -13.834  1.00  0.00           C
HETATM   28  C   LIG     1      -5.340  -0.600 -11.644  1.00  0.00           C
HETATM   29  N   LIG     1      -5.408  -3.127 -11.821  1.00  0.00           N
HETATM   30  N   LIG     1      -8.129  -2.821 -13.038  1.00  0.00           N
HETATM   31  O   LIG     1       3.072  -3.897  -0.073  1.00  0.00           O
HETATM   32  O   LIG     1       3.421  -0.244   0.755  1.00  0.00           O
HETATM   33  O   LIG     1       2.698  -1.683   3.334  1.00  0.00           O
HETATM   34  O   LIG     1       5.094  -3.922   1.972  1.00  0.00           O
HETATM   35  O   LIG     1      -6.061  -2.821  -9.650  1.00  0.00           O
HETATM   36  O   LIG     1      -8.120  -0.753 -12.076  1.00  0.00           O
HETATM   37  S   LIG     1       1.026  -4.388   1.646  1.00  0.00           S
HETATM   38  C   LIG     1      -9.774  -4.156 -16.922  1.00  0.00           C
HETATM   39  C   LIG     1      -9.839  -1.879 -16.189  1.00  0.00           C
HETATM   40  N   LIG     1      -9.383  -3.231 -15.850  1.00  0.00           N
HETATM   41  C   LIG     1      -9.964  -3.665 -14.556  1.00  0.00           C
HETATM   42  C   LIG     1      -9.568  -2.827 -13.310  1.00  0.00           C
CONECT    1   10   12   33
CONECT    2   31   12   37
CONECT    3   15   13
CONECT    4   17   15    5
CONECT    5    4   14
CONECT    6    7   15
CONECT    7    6    9   18
CONECT    8   28   29   27   26
CONECT    8
CONECT    9    7   23
CONECT   10   11    1   32
CONECT   11   13   10   31
CONECT   12    1    2   34
CONECT   13    3   14   11
CONECT   14    5   13
CONECT   15    6    4    3
CONECT   16   37
CONECT   17    4
CONECT   18    7   22
CONECT   19   23   22   20
CONECT   20   19   21
CONECT   21   20   24
CONECT   22   18   19
CONECT   23    9   19
CONECT   24   21   25
CONECT   25   35   24   29
CONECT   26    8   30   36
CONECT   27    8
CONECT   28    8
CONECT   29   25    8
CONECT   30   26   42
CONECT   31   11    2
CONECT   32   10
CONECT   33    1
CONECT   34   12
CONECT   35   25
CONECT   36   26
CONECT   37    2   16
CONECT   38   40
CONECT   39   40
CONECT   40   38   39   41
CONECT   41   40   42
CONECT   42   41   30
MASTER        0    0    0    0    0    0    0    0   42    0   42    0
END

```

picture here :

more about the compound at:

PDB Chemical Component 1YI

**UPDATE**

Last link doesnt work anymore :

> Ligand Expo Ligand Expo is down from January 1-31, 2023. Requests to this service will return a 410 http error code (“Gone”). Please help us better understand which Ligand Expo features are utilized and should be integrated with RCSB.org by taking this survey.

> Users can access CCD data in mmCIF format at https://files.wwpdb.org/pub/pdb/refdata/chem\_comp/. For other data download during server downtime, please contact chenghua.shao@rcsb.org.

Can use : https://www.ebi.ac.uk/pdbe-srv/pdbechem/chemicalCompound/show/1YI
> 2 votes

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


# Question
Title: Why do Illumina 850k/EPIC arrays ignore CpGs which are "GC" in the forward strand?
CpGs are symmetrical, in that a CG sequence on the forward strand is hybridized to a GC --- and both dinucleotides on each opposing strand are CpGs dinucleotides which can be methylated. Conversely, CpGs can be GC on the forward strand but CG on the reverse strand.

```
FORWARD -> 5'--CG--3'  [OR]  5'--GC--3' <- FORWARD
REVERSE -> 3'--GC--5'        3'--CG--5' <- REVERSE

```

The assignment of "forward" and "reverse" strandedness is more or less arbitrary.

Given the above, why does it seem like the Illumina 850k (aka EPIC) array only profiles methylation from CpGs which are CG in the forward strand, while ignoring CpGs which are GC in the forward strand? I would also love to hear if my premises are wrong.

Here is some R code illustrating the issue at hand:

suppressPackageStartupMessages({
  library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
  library(tidyverse)})
data(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)

IlluminaHumanMethylationEPICanno.ilm10b4.hg19 %>%
  getAnnotation() %>%
  as_tibble() %>%
  count(forward_seq=str_extract(Forward_Sequence, "\\[[ATCG]{2}\\]"))

# Results:
# # A tibble: 3 × 2
#   forward_seq      n
#   <chr>        <int>
# 1 [CA]          2922
# 2 [CG]        862927
# 3 [CT]            10
```

# Answer
My question was kindly answered by Biostars user Papyrus, who pointed out that a GC in the 5' to 3' direction is a GpC dinucleotide, not a CpG dinucleotide. GpC dinucleotides are then biologically and chemically distinct from CpG dinucleotides, due to the directionality of the DNA polymer.
> 1 votes

---
Tags: genome, database, sequence-annotation, genomics, public-databases
---


title: Matches, mismatches and indels

text:
# Question
Title: Matches, mismatches and indels
I have a bam-file of reads and reference genome. How can I located matches, mismatches and indels in three different files with Phred score in each data of the files?

# Answer
Try one of these:

GATK: https://gatk.broadinstitute.org/hc/en-us/articles/360036194592-Getting-started-with-GATK4

freebayes: https://bioinformaticsworkbook.org/dataAnalysis/VariantCalling/freebayes-dnaseq-workflow.html#gsc.tab=0
> 4 votes

---
Tags: sequence-alignment, indel, phred-scores
---


# Question
Title: hisat2 --rna-strandness option and downstream htseq-count analysis
I've got some doubts on the hisat2 `--rna-strandness` option and its output for downstream analysis. Is it expected to see a difference in alignment and counting results given the default usage of `--rna-strandness` in hisat2 followed by `htseq-count -s reverse` on a strand-specific assay?

Please see below.

I understand that the `--rna-strandness` option produces an XS tag to indicate where a transcript is from (on the + or - strand) for downstream transcriptome assembly analysis. I have a paired-end stranded sequencing library that was aligned to the genome using hisat2 without specifying the `--rna-strandness` (in other words, the default unstranded was the usage). Following this, the reads were assigned to genes using htseq-count and this time `-s reverse` was specified given the strand-specific sequencing assay type.

Since `--rna-strandness` is for transcriptome assembly using the XS tags generated and htseq does not use XS tags for counting, I presume there should be no practical impact from the above. Could you also shed light on this? in case I may have been overlooking other facts of the usages of the tools.

To help verify the above, I re-aligned and counted the reads from 2 samples by switching on `--rna-strandness RF` in hisat2. I attach the alignment and count features info. below for assessment.

Overall alignment rate of Sample 1:

```
94.52% (--rna-strandness RF)
94.12% (--rna-strandness unstranded)

```

Overall alignment rate of Sample 2:

```
94.57% (--rna-strandness RF)
94.15% (--rna-strandness unstranded)

```

Feature counts of Sample 1 (following --rna-strandness RF + -s reverse):

```
__no_feature        6327294
__ambiguous     2954776
__too_low_aQual     3784481
__not_aligned       688856
__alignment_not_unique      4858182

```

Feature counts of Sample 1 (following --rna-strandness unstranded + -s reverse):

```
__no_feature        6291151
__ambiguous     2911298
__too_low_aQual     4075017
__not_aligned       754400
__alignment_not_unique      16136045

```

Feature counts of Sample 2 (following --rna-strandness RF + -s reverse):

```
__no_feature        5417882
__ambiguous     1708510
__too_low_aQual     3532352
__not_aligned       564596
__alignment_not_unique      2859501

```

Feature counts of Sample 2 (following --rna-strandness unstranded + -s reverse):

```
__no_feature        5359434
__ambiguous     1676091
__too_low_aQual     3813344
__not_aligned       623122
__alignment_not_unique      2891792

```

These results look comparable to me across pipelines.

Thanks Guan

# Answer
If you reran the command with the correct settings, just leave it at that. (It is not at all clear to me that strandedness rf is correct)

If you want people to tell you if you ran the commands right, you need to put down what commands you used.
> 1 votes

# Answer
Sample 2 looks comparable but Sample 1 shows a drastic reduction in unique alignments.

Your strandedness param is going to dictate how your alignments are done in the sense that discordant pairs will either be detected or not.

I would identify the strandedness of your library and use the appropriate respective param.
> 1 votes

---
Tags: rna-seq
---


# Question
Title: How to get full headers with tabs as separators with segmasker?
I try to pass a *.fasta* file with headers where tabs are used as separators as input file to *NCBI BLAST+ segmasker*:
`>A0A059VAR9 IPR022879 V-type ATP synthase regulatory subunit B/beta`

Unfortunately, the *segmasker* output's headers are shortened to the first tab:
`>A0A059VAR9`

Is there a way to get whole headers in the output file (without replacing all the tabs in the input file with white spaces)?

# Answer
> 1 votes
I used a local version of `segmasker` *****. One way to solve this is,

```
perl -p -i -e 's/ //g' bioinfo.fa
perl -p -e 's/\t/^/g' bioinfo.fa | segmasker | perl -p -e 's/^/\t/g'  > bioinfo2.sm

```

*Edit: I changed `,` for `^` because the OPs code already has `,` and will insert an extra tab (\t). Same change in the line below*

Here the input file you presented in the link is called `bioinfo.fa`. If you need the pipe explaining let me know, but it's fairly straight forward.

If you want to preserve space characters, i.e. ' ' or "\s" - because in the code above I delete them - then (not tested):

```
perl -p -e 's/\t/^/g' bioinfo.fa | perl -p -e 's/ /_/g' | segmasker | perl -p -e 's/_/ /g' | perl -p -e 's/^/\t/g'  > bioinfo2.sm

```

The issue the OP encounter is not unique to tabs (\t) but also occurs with spaces (\s), at least with `makebasedb`.

**input**

```
>A0A000 IPR004839   Aminotransferase, class I/classII   PF00155 41  381
GVSAIDGMNTLVRGQKLPIFSSSGQPHSELAMQIARQASVPEEEEGGDDEEGSEFAVIFGAMGITQEEANEFMEDFERTGALERSVVFMNLADDPAVERTVTPRMVLTTAEYLAFEKDYHVLVILTDMTNYCEALREIGAAREEVPGRRGYPG
>A0A000 IPR010961   Tetrapyrrole biosynthesis, 5-aminolevulinic acid synthase   TIGR01821   12  391
...

```

**output**

```
>A0A000 IPR004839       Aminotransferase        class_I/classII        PF00155 41      381
41 - 53
>A0A000 IPR010961       Tetrapyrrole_biosynthesis       _5aminolevulinic_acid_synthase  TIGR01821       12      391
41 - 53
>A0A000 IPR015421       Pyridoxal phosphate dependent transferase        major domainG3DSA:3.40.640.10  48      288
106 - 126
...

```

Please don't ask about solving this via an upstream `makeblastdb` because that gets complicated. The rate limiting step, i.e. why it is ignoring tabs is because of `makeblastdb`. Thats a consequence of the way its written ... its a developer issue (see **blurb**)

---

***** Local `segmasker` installation ...

1. Install `conda` via `Anaconda` or `miniconda`
2. `conda create -n blast3 python=3.11` (creating an environment ... yep 3.11 is out (wow!, if you don't understand why **WOW!** ask in chat)
3. `conda activate blast3` (this is the conda environment, you can call it anything)
4. `conda install -c bioconda blast`

`perl | segmasker | perl` will now work. The Python installation is just good practice.

---

**blurb about NCBI's blast** The issue is `makeblastdb`, what happens is that long headers are a complete pain in algorithms because they use loads of RAM and you've got to carry them through reiterative loops (convergence etc ...). The way to solve it is via recoding (thats how I've done it). Essentially, you have a mini database forcing header recoding which activates at the start of the algorithm and reassembles the original header at the end of the algorithm. The user never sees this step, but it's really important for efficiency. Thus within the reiteration the header isn't consuming RAM, because it's an internal code, beyond base 10 to squash the overhead even further. The thing is header can become a lot bigger than the sequence data and thats not cool when you've loads and loads of data, e.g. a blast database. NCBI just haven't bothered, like if you can't see it why do it attitude. What they've said is that 50 characters is the limitation of the header for `makeblastdb`... what NCBI is looking for is a unique id at the start of the header it then truncates everything beyond that, i.e. tabs, spaces ... anything it can grab. NCBI have carried this through into `segmasker` it's not needed, but `makeblastdb` and `segmasker` would normally be used together - thats why. In the OPs case this results in duplicate headers, so the file format falls over. This isn't the OPs fault.

---
Tags: fasta, phylogenetics, segmasker
---


title: replacement has 1 row, data has 0

text:
# Question
Title: replacement has 1 row, data has 0
I have tried to use UpsetR to visualize the input file which can be found here

```
> library("UpSetR")
> orthogroups_df<- read.table("orthogroups.GeneCount.tsv",  header=T, stringsAsFactors = F)
> #All species
> selected_species <- colnames(orthogroups_df)[2:(ncol(orthogroups_df) -1)]
> selected_species
 [1] "Atha" "Cann" "NQLD" "Natt" "Ngla" "Nlab" "Nsyl" "Ntab" "Ntom" "Slyc" "Stub" "Vvin"
> head(orthogroups_df)
  Orthogroup Atha Cann NQLD Natt Ngla Nlab Nsyl Ntab Ntom Slyc Stub Vvin Total
1  OG0000000    0    0  965    0    0    3    0    0    0    0    0    0   968
2  OG0000001    0    1    3    0    0  448    0    0    0    0    0    0   452
3  OG0000002    0    1  313    0    0  120    1    0    1    0    0    0   436
4  OG0000003    0   93   15   21   46   16   33   63   36   25   39   26   413
5  OG0000004    1   42    2   34  109    6    8  154   11    9    4    0   380
6  OG0000005    0    2   61    1   34   44   91   70   43   20    1    0   367
> ncol(orthogroups_df)
[1] 14
> orthogroups_df[orthogroups_df > 0] <- 1
> upset(orthogroups_df,
+       nsets = ncol(orthogroups_df),
+       sets = rev(c(selected_species)),
+       queries = list(list(query = intersects, params = list("NQLD", "Nlab", "Nsyl"), color = "#238c45", active = T),
+                 list(query = intersects, params = list("NQLD", "Nlab"), color = "#ffd977", active = T)))

Error in `$<-.data.frame`(`*tmp*`, "freq", value = 45L) :
  replacement has 1 row, data has 0

```

How is it possible to fix the above error?

# Answer
I suspect the error is caused by a 'lack' of intersection between "NQLD", "Nlab", and "Nsyl"? If you remove "Nlab" from the first query the command appears to complete as expected, e.g.

```
library(UpSetR)
orthogroups_df <- read.table("~/Desktop/orthogroups.GeneCount.tsv",
                            header=T, stringsAsFactors = F)
selected_species <- colnames(orthogroups_df)[2:(ncol(orthogroups_df) -1)]
selected_species
#>  [1] "Athaliana_447_Araport11.protein_primaryTranscriptOnly"
#>  [2] "Cann"
#>  [3] "NQLD"
#>  [4] "Natt"
#>  [5] "Ngla"
#>  [6] "Nlab"
#>  [7] "Nsyl"
#>  [8] "Ntab"
#>  [9] "Ntom"
#> [10] "Slyc"
#> [11] "Stub"
#> [12] "Vvin"
head(orthogroups_df)
#>   Orthogroup Athaliana_447_Araport11.protein_primaryTranscriptOnly Cann NQLD
#> 1  OG0000000                                                     0    0  965
#> 2  OG0000001                                                     0    1    3
#> 3  OG0000002                                                     0    1  313
#> 4  OG0000003                                                     0   93   15
#> 5  OG0000004                                                     1   42    2
#> 6  OG0000005                                                     0    2   61
#>   Natt Ngla Nlab Nsyl Ntab Ntom Slyc Stub Vvin Total
#> 1    0    0    3    0    0    0    0    0    0   968
#> 2    0    0  448    0    0    0    0    0    0   452
#> 3    0    0  120    1    0    1    0    0    0   436
#> 4   21   46   16   33   63   36   25   39   26   413
#> 5   34  109    6    8  154   11    9    4    0   380
#> 6    1   34   44   91   70   43   20    1    0   367
ncol(orthogroups_df)
#> [1] 14
orthogroups_df[orthogroups_df > 0] <- 1

upset(orthogroups_df,
      nsets = ncol(orthogroups_df),
      sets = rev(c(selected_species)),
      queries = list(list(query = intersects,
                          params = list("NQLD", "Nsyl"),
                          color = "#238c45", active = T),
                     list(query = intersects,
                          params = list("NQLD", "Nlab"),
                          color = "#ffd977", active = T)))

```

Created on 2022-11-28 with reprex v2.0.2

You can see in the plot that "NQLD", "Nlab", and "Nsyl" don't appear to have any values in common, i.e. the 'dots' corresponding to those variables aren't joined. So, to fix the error you need to select variables that intersect in your queries.
> 2 votes

---
Tags: r, ggplot2, dataframe, upsetr
---


# Question
Title: How to convert Sars-CoV-2 ORF1ab codon positions to ORF1b positions?
Some websites mention amino acid mutations as ORF1ab, most however use ORF1a/b instead.

So how can I translate from one to the other? What do I need to subtract from say `ORF1ab:4588I` to get the corresponding `N` in `ORF1b:<N>I`?

# Answer
`ORF1ab:4588I` corresponds to `ORF1b:187I` so one needs to subtract `4401` to go from `1ab` to `1b` and add `4401` to go the other way.
> 2 votes

# Answer
There is now an Open Source online tool that does the conversion, Codon2Nucleotide by Theo Sanderson.
> 2 votes

---
Tags: sars-cov-2, covid-19, codon, nomenclature, orf
---


# Question
Title: Keep first 20 nucleotides of sequences using R
I have an Excel sheet from the CRISPR library where I have sequences of gRNAs (with 30 nucleotides) in a column and I only need to keep the first 20 nucleotides for those gRNAs and delete the rest nucleotides. The column has 1000 rows that have 30 nucleotide gRNA sequences. How can I delete extra nucleotides using R

# Answer
> 2 votes
You will need to "read" your data. There are packages that can read Excel files directly but I would save the data as csv first and then use `read.csv()`. The following would create a new column with your desired output:

```
data_frame$new_col <- substring(data_frame$col_with_seqs, 1, 20)

```

Alternatively you can just use Excel's `LEFT()` since your data is already in Excel:

```
=LEFT(A1, 20)

```

---
Tags: r, scrnaseq, bioconductor, crispr
---


# Question
Title: How to solve Nextflow error: "Trace file already exists"?
When trying to run epi2me-labs/wf-artic, I get the following error:

```
❯ nextflow run epi2me-labs/wf-artic \
 -r v0.3.18 --fastq  ~/Downloads/barcode95.fastq.gz \
 --scheme_version Midnight-ONT/V3

N E X T F L O W  ~  version 22.10.1
Launching `https://github.com/epi2me-labs/wf-artic` [ecstatic_marconi] DSL2 - revision: a60a1e1e73 [v0.3.18]
Trace file already exists: /Users/corneliusromer/Downloads/BJGMC_ONT_Run71/output/execution/trace.txt

```

It seems that nextflow doesn't like to rerun after it failed. Is there a convenient way to force nextflow to rerun, in a similar way to `snakemake --force-rerun`? Or do I need to manually delete the trace file every time? That works, but it's tedious.

# Answer
Rotation of report files (trace, execution, timeline etc) was removed recently (since version 22.10.0 in `a762ed5`). To overwrite the existing trace file you will need to add the following to your Nextflow config:

```
trace.overwrite = true

```

See: https://github.com/nextflow-io/nextflow/issues/3317
> 5 votes

---
Tags: assembly, sars-cov-2, nextflow, workflow-management
---


# Question
Title: What are the conditions to perform Gene Expression Analysis on data you've performed RNA Sequence Analysis on?
I want to perform Differential Gene Expression Analysis. I recently completed RNA Sequence Analysis using the Seraut Tutorial up to *this* point: Finding Differentially Expressed Features (Cluster Biomarkers).

My project can be found at *this* Rpubs Link which contains the link to the dataset on the 10x genomics website.

**My question is... is the following table derived from performing the steps found in the first link I provided?**

**If so, once I derive that table ( ID, Name, L2FC, p-value, etc.), does that mean I can perform Differential Gene Expression Analysis?**

I am a complete newbie and am teaching myself all of this. Thanks, y'all.

**UPDATE:** Yes, it did provide me that table. With this table, am I able to perform gene expression?

# Answer
The table pictured here *is* the result of a differential gene expression (DGE) test. In the table, each row is a gene (labeled with the gene ID and name) and each column (besides the first two) are the results of DGE tests. Each cluster has two sub-columns: L2FC (log2 fold-change) which tells you the average relative expression changes between cells of the cluster vs. others; and p-value, which is how significant the change is.

You can make this table yourself by performing DGE tests for all clusters. You've already done this yourself in your Rpubs script! To make it easier though, so you don't have to manually test each cluster, use the FindAllMarkers function. It will compare each cluster to all other clusters and report the L2FCs and p-values.

However, I recommend using the p\_val\_adj (adjusted p-value) instead of the p-value. This is due to a problem called multiple testing: you are performing many, many tests (each gene in each cluster is a separate statistical test), and we expect 5% of differentially expressed genes (DEGs) to be DE by random chance. Over something like 100,000 tests, that means you'd expect 5,000 of those to be false positives (which is huge). So p-value is adjusted to account for multiple tests.

Once you filter the DEG table (p-adj\<0.05), you can do other downstream analyses to describe what has changed between the conditions. You already have a good start in your Rpub script with the violin plots. Cluster 6 is clearly distinct in terms of those genes' expression levels. You could try to figure out what's special with that cluster. Good luck!
> 4 votes

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