# Question
Title: Why ORF13 and ORF14 of Bat Sars coronavirus Rp3 have no corrispondence in Sars2?
I'am comparing bats coronavirus with Sars2, i.e. COVID-19. I have one question...

* one Bat Sars coronaviurs identify in GenBank as DQ071615.1 has the ORF13(28122..28415) and ORF14(28572..28784) proteins
* Sars2 identify in GenBank as MN908947\_3, has range 28284..28577 not labeled as ORF13 and the same 28734..28955 is not labeled as ORF14?

The ORF13(28122..28415) of Bat coronavirus and PSEUDO\_ORF13(28284..28577) of sars2 share 72% homology and ORF14(28572..28784) of Bat coronavirus and PSEUDO\_ORF14(28734..28955) of sars2 share 75% homology

---

This seems one publication for ORF14 and ORF13 in Bat Viruses and relations in Sars2 here

The author say that ORF14 in Sars2 is all contained in N protein... And ORF13 is know in that above link as ORF9B and in that written is said ORF9B is implied in immunitary innate suppression...

# Answer
> 1 votes
The virus Rp3 stands for Rhi pea isolated in China in 2004.

The trees is here and groups with the Bulgarian sarbecoviruses and is closely related to the African sarbecoviruses. Thus although it is a Chinese sarbecovirus is not with the "Chinease originating" clade as is SARS and Covid-19.

Thats why the more mysterious ORFs ORF13 and 14 are missing, because over notable genetic distances these genes along with ORF10 are present/absent. This is a common feature of the betacoronaviruses.

In summmary this virus is notably distance from the major human pathogenic SARS/covid-19 viruses. Thus the absence of these auxilary genes isn't really surprising.

Its the very last taxa on the SARS-related viruses clade (monophyly) and is in a different group to SARS, which is the Beta-1 clade on this most excellent paper.

---

To answer the comment. They are outgroups, basically ancestral to COVID-19 and SARS-1 clade. So they lost these genes whilst the ingroup including COVID-109 and SARS-1 retained ORF13, 14 ... but I think lost ORF10 from memory

---
Tags: covid-19, sars-cov-2, coronavirus
---


# Question
Title: How to print alternative allele from VCF with htslib c++
I would like to print the alternative allele from a vcf file which i load in with htslib. The description of some functions: https://github.com/samtools/htslib/blob/develop/htslib/vcf.h

I have used bcf1\_t objects by creating a header and reading the information.

```
htsFile *test_vcf = NULL;

  // creates header
  bcf_hdr_t *test_header = NULL;

  // initialize and allocate bcf1_t object
  bcf1_t *test_record = bcf_init();

  test_vcf = vcf_open("Testfile.vcf", "r");

  // returning a bcf_hdr_t struct
  test_header = bcf_hdr_read(test_vcf);

  while(bcf_read(test_vcf, test_header, test_record) == 0){
    std::cout << "chr" << bcf_hdr_id2name(test_header, test_record->rid) << ":" <<test_record->pos+1 << std::endl;
  }
  bcf_hdr_destroy(test_header);
  bcf_destroy(test_record);
  bcf_close(test_vcf);

```

But im not able to extract the reference or alternative allele information using this structure bcf1\_t.

As far as i can read its perhaps possible with the structure and function:

```
bcf_dec_t
bcf_unpack()

```

But i am not able to get it to work since all the other functions seems limited to the bcf1\_t struct.

Can anyone help with extracting the REF and ALT sequence? Thanks for your time

# Answer
> 1 votes
see the code in https://github.com/samtools/htslib/blob/develop/vcf.c#L3104

```
 bcf_unpack((bcf1_t*)v, BCF_UN_ALL);// if you need genotypes too, might use another constant to avoid decoding all the structure
 (...)
 kputc('\t', s); // ALT
 if (v->n_allele > 1) {
        for (i = 1; i < v->n_allele; ++i) {
            if (i > 1) kputc(',', s);
            kputs(v->d.allele[i], s);
        }
    }

```

---
Tags: vcf, c++, htslib
---


# Question
Title: How to replace sequence identifiers in a fasta with OTU IDs from another file?
I'm pretty new to Unix and bioinformatics and having a hard time accomplishing the following. I have one FASTA file with sequences and headers, and one OTU ID mapping file (.txt) with OTU IDs and their sequences. The sequence header in the .fasta file partially matches the sequence listed in the OTU ID file. What I want to do is replace the current headers in the FASTA file with just the OTU IDs from the .txt file. The files look like something like this:

OTU\_ID.txt

```
denovo123 U1019.1010.2020_20384
denovo453 U1472.2992.2020_20443
denovo763 U1928.1898.2020_20343

```

.FASTA

```
>U1019.1010.2020_20384 0000-AHDJ-38347_223
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA
>U1928.1898.2020_20343 0000-AHDJ-24247_242
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA

```

and my desired output is that the .fasta file would look like this:

```
>denovo123
GAAATACCGAGAGAGAAATACCGAGAGAGAAAGAAATACCGAGAGAGAAATACCGAGAAATACCGAGAGAGAAATACCGAGAGA
>denovo763
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA

```

Is it possible to do this in the command line?

# Answer
> 1 votes
You can do it with bash quiet easy. First create a backup of your fasta file with

```
cp otu.fasta otu_backup.fasta

```

Then you can run:

```
while IFS=" " read -r otu fastaId rest;
do
sed -E -i "s/>$fastaId .+/>$otu/" otu.fasta; done < ids.txt;
done < OTU_ID.txt

```

What does the code?

With `IFS=" "` you declare a whitespace as field seperator and by reading each line of your `OTU_ID.txt`, you have the two variables `otu` and `fastaId` which hold the otu name and the id of the Fasta Header.

Within the `while` loop, sed is looking each loop for the corresponding id within a RegEx and substitute it with the corresponding otu name. By using the `-i` option, you overwrite the existing file. This is why you need to create a backup.

# Answer
> 1 votes
You can put you fasta file in a single line per sequence:

```
perl -pe 'if(/\>/){s/\n/\t/; s/>//; s/\n//}' file.fasta > file.tsv

```

This line will create the next output:

```
U1019.1010.2020_20384 0000-AHDJ-38347_223   GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA
U1928.1898.2020_20343 0000-AHDJ-24247_242   GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA

```

Then with awk you can match columns separated with tab:

```
awk -F' ' 'FNR==NR{a[$2]=$1;next} ($1 in a) {print ">"a[$1]"\n"$3}' otu_id.txt file.tsv

```

You will get this output:

```
>denovo123
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA
>denovo763
AAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA

```

You can add **`$1`** to check if you are getting the right id.

```
awk -F' ' 'FNR==NR{a[$2]=$1;next} ($1 in a) {print ">"a[$1], $1, "\n"$3}' otu_id.txt file.tsv

```

If everything is alright you can redirect to a new file:

```
awk -F' ' 'FNR==NR{a[$2]=$1;next} ($1 in a) {print ">"a[$1]"\n"$3}' otu_id.txt file.tsv > file_ids.fasta

```

How it works, well `awk` implicitly loops through each file, one line at a time.

`FNR==NR{a[$2]=$1;next}`

`NR` is the number of lines that awk has read so far and `FNR` is the number of lines that `awk` has read from the current file. Thus, if `FNR==NR`, we are still reading the first-named file: otu\_id.txt. For every line in file.tsv, we assign `a[$2]=$1`.

Here, `a` is an associative array, and `a[$2]=$1` means saving file.tsv first column, denoted `$1`, as `a` value in array a using file.tsv second column, `$2`, as the key.

`next` tells `awk` to skip the rest of the commands and start over with the next line.

`($1 in a) {print ">"a[$1]"\n"$3}`

This means that we are reading the second file: otu\_id.tsv. If match the first field of the line in file.tsv, as determined by the contents of array `a`, it will print out a line with the values of fields `a[$1]` and `$3` separated with a new line character.

---
Tags: fasta, shell, identifiers, grep
---


# Question
Title: Recommended coverage depth for ChIP Seq (H3K9me3 )
According to ENCODE recommendations

> For broad-peak histone experiments, each replicate should have 45 million usable fragments

According to my understanding, a usable fragment means a uniquely mapped read / pair of reads.

However

> H3K9me3 is an exception as it is enriched in repetitive regions of the genome. ... Tissues and primary cells should have 45 million total mapped reads per replicate.

I fail to grasp Encode's language.

* Does "total mapped reads" mean uniquely mapped reads?
* If so, why is 45 million "an exception" ? It is the same as the other broad peaks.

Any ideas where else can I look for an explanation?

# Answer
> 1 votes
The ENCODE Data Coordination Center, that developed and maintains the ENCODE portal, have implemented scripts validating adherence to the standards the ENCODE Consortium agreed upon when evaluating ChIP-seq analysis results. For any ChIP-seq experiment, except for H3K9me3, the adherence to the read depth standard is assessed looking on the number of "usable fragemnts" (https://www.encodeproject.org/data-standards/terms/#read-depth) (see answer by James Hawley). For H3K9me3 the read depth is assessed not by looking on usable fragments, but by looking on the total number of mapped reads (if it is a PE sequencing run, the number of mapped read pairs). It is an exception, because the number of mapped reads, especially for H3K9me3 is substantially higher than the number of usable fragments (where the duplicated and low mapping quality reads are filtered out). The decision to apply different read depth standard in case of H3K9me3 came up due to the logic explained here (https://www.encodeproject.org/chip-seq/histone/), briefly: because the H3K9me3 mark is enriched in repetitive regions of the genome, the number of usable fragments (post-filtering) is not expected to be as high as other marks (45M), despite high-quality libraries and sufficient sequencing depth. Instead of marking these experiments as failing read depth standard, the script checks the total number of mapped reads (or read pairs for PE sequencing runs) and makes sure it is above the read depth threshold (45M).

# Answer
> 1 votes
In the initial description for "usable fragments", ENCODE has the following (emphasis mine):

> Usable fragment – A fragment is defined as the sequencing output **corresponding to one location in the genome**. If single-ended sequencing is performed, one read is considered a fragment. If paired-ended sequencing is performed, one pair of reads is considered a fragment. **Fragments are considered usable if they pass the various filters in the ChIP-seq uniform processing pipelines**. Used to evaluate ChIP-seq data.

Aligning reads to repetitive regions is difficult because they don't give unique alignments. Reads from these regions will often have poor alignment scores and may not pass alignment QC checks. Reads may also be counted as multi-mapped, and often methods call for uniquely mapped reads, which again, may not be true for a large number of reads from H3K9me3 ChIP-seq.

So for this particular histone mark, ENCODE loosens the definition for what you should count towards a minimum sequencing depth threshold.

---
Tags: chip-seq
---


# Question
Title: How would we test for significant changes in expression when many of the normalised counts are clustered together?
I've been handed RNA-seq data with a lot of associated covariates. This data has been put through the DESeq2 package and as a result normalised the data. One of the transcripts of interest still has non-normal distribution for two genotypes. This is because many of the normalised counts are now clustered around a value of 4.

This is making it tricky to perform statistical tests because in order to perform parametric tests (in this case I want to perform a ANCOVA or multiple linear regression to include the covariates) the data must be normally distributed.

I also don't want to use a non-parametric test because I can't include all the covariates.

Should I just perform parametric tests ignoring this violation or do I have to just perform non-parametric on all my data for a fair comparrison? Or is there another alternative option I am missing?

Thanks in advance

# Answer
TLDR- just do the analysis in DESeq2 (it can handle multiple covariates)

If I've interpreted your question correctly there seem to be a few confusions here.

First off, for a multiple linear regression aka ANCOVA, one of the assumptions are the *residuals* are normally distributed - however this is almost never an issue in practice. The distribution of the data is not so important, you make the assumption that the data has come from a population which is continuous and normally distributed that's why you've chosen to do an ANCOVA. But, RNAseq produces count data - not continuous and bounded at zero. You need to be looking at using a Poisson GLM, or negative binomial GLM as done by DESeq2 (lookup overdispersion)

The normalisation done by DESeq is to account for different library sizes,it doesn't *transform* the data to make it conform to a normal distribution. You'd need to do some form of data transformation to do that, but don't it's unnecessary.

I'd suggest reading the DESeq vignettes and following their standard pipeline.

Hope that is of some help.
> 1 votes

# Answer
Your assessment is perfectly correct. The easiest approach is to perform a transformation to normalise your data and then you are good to go. DeSeq2 has a standard transformation (I think it is called deseq2) and you can check your results via a Q-Q plot and you are looking for an approximately linear relationship.

There are a multitude of transformations, but the transformation proceedure has been heavily investigated and standardised in RNAseq

---

If normalisation by deseq2 has failed (Q-Q) you simply examine an alternative transformation. However, do not transformation a transformation. If you cannot normalise the data then you are stuck with non-parametric stats in my personal opinion.

However, mixed models can work if you trap the underlying distribution and in a regression will result in a very low residual.
> 0 votes

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


title: FASTA and PDB: How to specify chain?

text:
# Question
Title: FASTA and PDB: How to specify chain?
For proteins that have multiple chains (e.g., 1EMS), is there an easy way to specify which chain I want to use for blastp?

I cannot imagine that I am the first person to have this problem, but so far my only solution has been to write my own script that gets the chain I want from the .pdb file and convert that to a .fasta file.

# Answer
WEB VERSION OF BLASTP

When you run a blastp query using a fasta file with multiple chains you can choose which chain's results to examine using a drop-down box that appears **after** the search has run.

TERMINAL BLASTP

The results for each chain will be listed separately. For ease of searching you can specify the output to be a csv of your desired format (See blastp -help and look up the -outfmt option for more information).

Note that for some proteins, such as the 1EMS you chose for an example, the multiple chains are identical. When you get the fasta file for that specific protein it will only have one entry. However if you really need a fasta file formated so each chain has its own sequence you can get the fasta file for all PDB entries and extract the specific protein\_chain you want.
> 0 votes

# Answer
The NCBI BLAST can be used with the PDB DB (which NCBI has). The PDB codes are stored as 4 letter codes underscore chain, *e.g.* `1GFL_B`. The catch is segment identifiers —but generally they are the same peptide so shouldn't be an issue.

You can search specifically the PDB DB in NCBI (not the RCSB PDB) by setting the database to `PDB`.
> 3 votes

---
Tags: fasta, pdb, blastp
---


title: Rejected publication

text:
# Question
Title: Rejected publication
A publication that we have been working on for over 6 months on COVID-19 has been rejected by 2 journals. The objective of our research was to identify effective compounds against the COVID-19 virus by using Molecular Docking studies performed with AutoDock Vina. Which Journal should we submit our manuscript to? Following are the Rejection emails I received from them. What do you think is the reason for the rejection and what should I do to resolve the issue to get my manuscript successfully published.

***Rejection Notification from International Journal of Molecular Sciences*** "We are writing to inform you that we will not be able to process your submission further. Submissions sent for peer-review are selected on the basis of discipline, novelty and general significance, in addition to the usual criteria for publication in scholarly journals. Therefore, our decision is not necessarily a reflection of the quality of your work."

***Rejection Notification from International Journal of Biological Macromolecules*** "Thank you for submitting your manuscript to International Journal of Biological Macromolecules. I have reviewed your manuscript and I must regretfully inform you that your work presents information that is specific to a selective audience in its nature and interest and is more appropriate for a more specialized journal. The focus of your present work, a computational study, is unfortunately outside the scope of the journal, therefore, unsuitable for publication in IJBIOMAC. For IJBIOMAC publication, the computational analysis should be used only as a supporting tool to non-computational analysis. We appreciate you submitting your manuscript to International Journal of Biological Macromolecules and thank you for giving us the opportunity to consider your work."

Thank you for your time.

# Answer
> 2 votes
I cannot comment on your specific manuscript since I obviously cannot read it but in general the question is whether a study does add an interesting novelty.

The problem with pure in silico studies is that there is essentially no confirmation towards the in vivo relevance. Drug companies simulate and also test many thousands of compounds in a high-throughput fashion on a daily basis so it is often questionable to what extend individual pure in silico studies without in vivo confirmation do actually add to the bigger picture.

That *might* be the reason your work was rejected but as said this is a rather general comment and I can neither judge content nor quality of your work so please do not feel offended and read it with a grain of salt.

---
Tags: covid-19, publishing
---


# Question
Title: Assessing PyMol sequence alignment object
I use cealign in PyMol for structure alignment. Instead of visualization, I want to return the alignment object to my python script for further analysis. Is there a function to return the object without I/O to a clustalw file?

edit:

By alignment object, I meant the sequence alignment in cealign.

# Answer
In the application, unfortunately the only option is doing

* `cealign ,, object =`
* `Menu > Export Alignment...` on your cealign  object
* and converting the Clustal format using one of the millions of tools online

## get\_fastastr

As you know `cealign` adds gaps, but `cmd.get_fastastr('')` strips them. There is no option to keep the gaps —annoyingly especially with missing loop info! The help has this to say:

```
 get_fastastr(selection='all', state=-1, quiet=1, key='', _self=<module 'pymol.cmd' from '/Users/matteo/miniconda3/lib/python3.7/site-packages/pymol/cmd.py'>)

```

## Correction and dissection

> So I had a quick check to see if there is an alternative route to get the sequence string with gaps. These are my observations using the Python3 module, sorry if they aren't polished up.

In the comments I said something incorrect. `cmd.cealign` returns a dictionary with the 4x4 rototranslational matrix —no need for an object. Using Python3 and two TIM barrels as an example:

```
import pymol2

with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1TIM', 'tim')
    pymol.cmd.remove('not chain A')
    pymol.cmd.fetch('5AHE', 'hisA')
    x = pymol.cmd.cealign('tim', 'hisA')
    print(x)
# {'alignment_length': 184, 'RMSD': 3.8800425577361755, 'rotation_matrix': [...]}

```

While the object does not contain any useful info from PyMol *directly*\*:

```
    ...
    pymol.cmd.cealign('tim', 'hisA', object ='rototrans')
    pymol.cmd.get_session('rototrans')['names'][0][5]
    # useful info is generally in ['names'][0][5]

```

(\*) The code is actually wrapped C++ code, so the object `rototrans` is in the C++ layer and can only be access by session or possibly by some functions in `cmd.exporting`. In fact, it can be exported to an aln string (ClustalW).

```
    aln = pymol.cmd.exporting.get_alnstr('rototrans')

```

In the PyMOL application, doing the `cealign` with an object (`cmd.cealign('tim', 'hisA', object ='rototrans')`, allows the Menu\>Export Alignment...\>'rototrans' which is saved as a ClustalW file, as you mention. Namely:

```
CLUSTAL

tim          APRKFFVGGN----WK-------MNGKRKSLGELIHTLDGAKLSADTEVVCGAPSI------Y
hisA         MIIPALDLIDGTVVRLHQGDYARQRDYGN---DPLPRLQDYAAQGAGVLHLVDLTGAKDPAKR
                           ..       ......   .....*..................      .

tim          LDFARQK---LDAKIGVAAQNCYKVPKGAFTGEISPAMIKDI-GAAWVILGHSERRHVFGESD
hisA         QIPLIKTLVAGVNVPVQVGGGVRTEE--------DVAALLKAGVARVVIGS---------TAV
             .......   ................        ..*..... .*..**..         ...

tim          ELIGQKVAHALAEG---LGVIACIGEKLDEREAGIT-----------EKVVFQETKAIADNVK
hisA         KSPDV-VKGWFERFGAQALVLALDV---RIDEHGTKQVAVSGWQENSGVSLEQLVETYLPVGL
             ..... *.......   ..*.*...   ...*.*..           .....*..........

tim          DWSKVVLAYEPVWAIGTGKTATP------QQAQEVHEKLRGWLKTHVSDAVAVQSRIIYGGSV
hisA         ----KHVLCTDI----------SRDGTLAGSNVSLYEEVCARYP---------QIAFQSSGGI
                 ........          .      .......*.......         *......*..

tim          TGGNCKELASQHDVDGFLVGGASLKPEFVDIINAKH------------
hisA         GDIDDIAALRGTGVRGVIVGRALLEGKFTVKEAIQ?WQNVKGHHHHHH
             .............*.*..**.*.*...*..

```

However, there is no method in `cmd.exporting` that gives out just the sequence with gaps.

Whereas the `cmd.exporting` has the following... none applicable.

```
 print(dir(pymol.cmd.exporting))
['DEFAULT_ERROR', 'DEFAULT_SUCCESS', 'QuietException', 'SDF', 'SDFRec', 'Shortcut', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__spec__', '_cmd', '_feedback', '_get_dump_str', '_get_mtl_obj', '_raising', '_resn_to_aa', '_session_convert_legacy', '_unit2px', 'absolute_import', 'assign_atom_types', 'cPickle', 'cache', 'cache_action_dict', 'cache_action_sc', 'cmd', 'copy', 'copy_image', 'export_coords', 'fb_mask', 'fb_module', 'get_alnstr', 'get_bytes', 'get_ccp4str', 'get_cifstr', 'get_fastastr', 'get_maestr', 'get_mol2str', 'get_pdbstr', 'get_pqrstr', 'get_psestr', 'get_sdfstr', 'get_session', 'get_str', 'get_xyzstr', 'io', 'is_dict', 'is_error', 'is_list', 'is_ok', 'is_tuple', 'loadable', 'lock', 'multifilenamegen', 'multifilesave', 'multisave', 'os', 'png', 'print_function', 'pymol', 're', 'save', 'savefunctions', 'selector', 'sys', 'thread', 'traceback', 'unlock']

```
> 0 votes

---
Tags: sequence-alignment, pymol
---


# Question
Title: How can I locate duplicated regions in a sequence?
I am facing an issue when trying to align short reads against a region in human chr5. The two Sensory Motor Neuron genes, (SMN1 and SMN2) are almost 100% identical and this causes the aligner to fail to align reads correctly since each read matches perfectly to two separate regions of the reference, leading to very low alignment scores.

If I BLAT the SMN1 gene's full genomic sequence against the human genome (hg19) on UCSC, I get (only showing the 1st two lines):

```
QUERY       SCORE START   END QSIZE IDENTITY  CHROM  STRAND  START       END   SPAN
--------------------------------------------------------------------------------------------------------------
NG_008691.1 35072     1 35072 35072   100.0%  chr5   +    70215768  70250839  35072
NG_008691.1 34279     1 35072 35072    99.8%  chr5   +    69340350  69375383  35034

```

And on the browser:

To get around this, I want to mask the SMN2 gene. However, the region of similarity extends beyond the annotated boundaries of the gene, so I would like to identify the full extent of this duplication.

So, given a large nucleotide sequence in fasta format, how can I compare it to itself and identify identical or near identical regions?

# Answer
If you know approximately where are the duplicated regions, you can cut them out from the genome and align the two block using tools like MUMmer or Mauve. Although Mauve is not scalable solution (You can not scan genomes for duplicates with it) it's very handy for aligning small(ish) blocks. It can also align multiple sequences if you would like to check some tri/multiplications.

I downloaded the sequence of ch5 you linked and extracted big regions around the genes (300bp blocks). The first block is from 70100000 - 70400000 (with SMN1 sequence at 115768 - 150839) and the second from 69200000 - 69500000 (with SMN2 sequence at 140350 - 175383).

I used progressive mauve to align the two (click on File -\> align with ProgressiveMauve -\> add Sequence ... -\> Align).

I found that the duplicated region is almost 200k long and maps between 44690 - 225520 of the first block and 19760 - 200940 of the second block, which we still need translate to the original genome coordinates, but I believe I can leave this to you :-).

P.S. you can also somehow upload annotations to Mauve, but I was lazy, so I "hand drawn" the arrow where I knew the genes are.
> 2 votes

# Answer
GeneBlocks is a Python library for comparing DNA sequences. It can be used to find common blocks (identical regions) in a group of DNA sequences. Using a modified version of your example (I duplicated a small region):

```
import geneblocks

sequence = geneblocks.load_record('sequence.gb', name='seq')
sequence2 = geneblocks.load_record('sequence.gb', name='seq2')
blocks = geneblocks.CommonBlocks.from_sequences([sequence, sequence2])

blocks.plot_common_blocks()

```

Note that it wasn't designed to plot only one sequence, so as a workaround the sequence is given twice, with different names. As an alternative, you can use one sequence and plot with DNA Features Viewer:

```
import dna_features_viewer

blocks = geneblocks.CommonBlocks.from_sequences([sequence])
seq_to_plot = blocks.sequences_with_annotated_blocks()["seq"]
graphic_record = dna_features_viewer.BiopythonTranslator().translate_record(seq_to_plot)
ax, _ = graphic_record.plot(figure_width=10, strand_in_label_threshold=7)

```

See the documentation for more details. A web interface is also provided at EGF CUBA: find common blocks.

*Disclaimer: I'm the current maintainer of GeneBlocks*
> 2 votes

---
Tags: fasta, genome, sequence-analysis, sequence-homology, duplications
---


# Question
Title: What is the best way to drop alignment artifacts in sequence objects in biopython?
I have a bunch of sequences created with `Roary` which I want to preprocess. When I read the sequences (genomics data) I want to translate them to protein sequences using:

```
import pandas as pd
from Bio.Seq import Seq
from Bio import SeqIO

def read_fasta_file(fn):
    with open(fn,'r') as file:
        fasta_sequences = SeqIO.parse(file, 'fasta')
        sequences = pd.DataFrame([{
            'filename': basename(fn),
            'id': rec.id,
            'seq': str(rec.seq.translate())} for rec in fasta_sequences])
    return sequences

```

The problem is that `.translate()` fails with:

> TranslationError: Codon 'T--' is invalid

Is there a build-in function to ignore the dashes?

# Answer
> 2 votes
I think the `.ungap()` is what I was looking for.

```
def read_fasta_file(fn):
    with open(fn,'r') as file:
        fasta_sequences = SeqIO.parse(file, 'fasta')
        sequences = pd.DataFrame([{
            'filename': basename(fn),
            'id': rec.id,
            'seq': str(rec.seq.ungap().translate())} for rec in fasta_sequences])
    return sequences

```

---
Tags: python, seqio
---


title: Convert a fasta file to Newick format

text:
# Question
Title: Convert a fasta file to Newick format
My problem is (should be) pretty simple.

I have a series of 10 sequences which I wish to convert into a Newick format file - what would be good software to do this? I've tried using Jalview and am getting nowhere.

Any help appreciated.

# Answer
> 3 votes
The best would be Biopython using AlignIO, given they are all in a single file otherise SeqIO. The alignment approach,

```
from Bio import AlignIO
align = AlignIO.read("/home/path/myseq.fa", "fasta")
print(align)
with open("/home/path/example.nex", "w") as handle:
    SeqIO.write(align, handle, "nexus")

```

(not debugged)

MEGAX or ClustalX which are nice GUI's that will also work. If you've not got Biopython installed then ClustalX is probably the easiest.

# Answer
> 1 votes
@M\__'s answer led me to clustal.org and some searching led me to this page (Christoph's).

I tried compiling Clustalx (GUI - problem with qmake) which didn't work, but Clustalw (compiled to exe clustalw2) - the command line tool - worked a treat.

So, I had my `clustalw2` executable and my 25 *Drosophila* alcohol dehydrogenase protein sequences and I ran the following commands (from Christoph's page):

```
./clustalw2 -align -type=protein -infile=droso.fasta

```

This produces a droso.aln file which looks like this (only a few lines shown):

```
D_paulistorum               ETRVAEKLLEHPTQTTLACAQNFVKAIELNKNGAIWKLDLGTLEPIEWTK
D_insularis                 ETRVAEKLLEHPTQTTLACAQNFVKAIELNKNGAIWKLDLGTLEPIEWTK
D_equinoxialis              ETRVAEKLLEHPTQTTLACAQNFVKAIELNKNGAIWKLDLGTLEPIEWTK
D_willistoni                ETRVAEKLLEHPTQTTLACAQNFVKAIELNKNGAIWKLDLGTLEPIEWTK
                            *. *:**** ** *.:  *.:*********:*****:** .***.* *::

```

Then on Christoph's page, there is a section:

* Generate a Newick tree file

So, I followed the instructions and ran:

```
./clustalw2 -tree -infile=droso.aln

```

and obtained a file droso.ph which looks like a Newick format file to me (some lines):

```
(
(
(
D_madeirensis:0.00000,
D_subobscura:0.00787)
:0.00257,
D_guanche:0.00530)
:0.00336,
(
(
D_pseudoobscura:0.00000,
D_pseudoobscura_pseudo:0.00000)
:0.00262,
D_persimilis:0.00131)
:0.01239)
:0.01132,
D_ambigua:0.01821)
:0.00839,

```

So, it looks like I have a path to doing what I wanted in the first place - getting from .fasta to Newick format.

Then, using the programs jalview and figtree, I was able to draw trees and the like!

# Answer
> 1 votes
I think you can do this in Jalview. In your alignment window, open the **Calculate** menu and select **Calculate Tree or PCA...**. This will open a **Choose Calculation** window which you can use to select what kind of tree to build. When you click on the **Calculate** button you will get another window with the tree in it. In that window you can open the **File** menu and then select **Save as...****Newick Format**. Hopefully this will do what you need.

---
Tags: phylogenetics, fasta, phylogeny, format-conversion
---


# Question
Title: Is there a tool for synonymous recoding of DNA sequences?
I've got a DNA sequence that I'd like to make synonymous mutations throughout, thereby preserving the amino acid sequence. Does anyone know of a tool to achieve this

# Answer
There seems to be a bioconductor package named `synmut` which does just that. Apparently, it can take codon usage into account when generating synonymous mutations.

Copy pasting an example that matches your use case from the official doc:

```
library("SynMut")
filepath.fasta <- system.file("extdata", "example.fasta", package = "SynMut")
filepath.csv <- system.file("extdata", "target_regions.csv", package = "SynMut")
region <- read.csv(filepath.csv)
rgd.seq <- input_seq(filepath.fasta, region)
mut.seq <- codon_random(rgd.seq, keep = TRUE)
Biostrings::writeXStringSet(get_dna(mut.seq), "rgd.fasta")

```
> 1 votes

# Answer
DNA Chisel can introduce synonymous mutations into your sequence. Annotate the coding sequence with the @cds constraint and the ~change objective.

Additionally, you may also optimize the codon usage to match your species with ~match\_codon\_usage or one of the other options.

A web interface is also provided at EGF CUBA (Collection of Useful Biological Apps): Sculpt a Sequence

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

---
Tags: r, bioconductor, mutations, synonymous-mutations
---


# Question
Title: Question of Padj value (very high, infinite) got from dds
After I got dds from my count matrix by DESeq(), I use results() method to get Padj and log2 fold changes for volcano plotting, then I found some gene's Padj values are infinite and the P-value is very close to 0; this makes those dots on my volcano plotting located very top to Y-aix and also confused me:

1. is this normal or not, such a p-value (very close to 0) and -logPadj is infinite;
2. how should I treat those data: are they very different expressed between samples and control or just abnormal data which should be taken off?

Here follows the table:

table 1 is dds results, -log10(padj) show inf; table 2 is count matrix from featureCount table 3 is size factor of samples

Please help me !

# Answer
> 0 votes
1. The fold changes for those are huge (on the order of 32x for the largest), so it's unsurprising for the adjusted p-values to be absurdly small.
2. There's no reason to treat them differently.

---
Tags: rna-seq, deseq2
---


# Question
Title: Metagenomics: Identifying most common sequences
I am working on a project and used the following command:

```
vsearch --derep_fulllength filtered_merged.fa -sizeout -relabel Uniq -output dereplicated_filtered_merged.fa

```

and got the following output:

```
87373926 nt in 203453 seqs, min 310, max 480, avg 352
Sorting 100%
10981 unique sequences, avg cluster 2.0, median 1, max 1287
Writing output file 100%

```

The output had provided me with the data that 10981 unique sequences have been identified. But I cant seem to identify how many reads of the most common sequence were present in the input data.

Any suggestions will be kindly appreciated!

# Answer
According to the VSEARCH docs, since you have specified `--sizeout` your abundances have been written into the FASTA headers:

> --sizeout
>
> Take into account the abundance annotations present in the input fasta file (search for the pattern ’\[\>;\]size=integer\[;\]’ in sequence headers). That option is active by default when rereplicating.
>
> **Add abundance annotations to the output fasta file (add the pattern ’;size=integer;’ to sequence headers).** If --sizein is specified, each unique sequence receives a new abun- dance value corresponding to its total abundance (sum of the abundances of its occur- rences). If --sizein is not specified, input abundances are set to 1, and each unique sequence receives a new abundance value corresponding to its number of occurrences in the input file.
> 1 votes

---
Tags: fasta, clustering, metagenome
---


# Question
Title: How do you visualize minimap2's PAF output format?
As described by the developer of minimap2:

> PAF is a text format describing the approximate mapping positions between two set of sequences.

He also suggest to use either Ribbon to visualize PAF files, or to use `minidot` from `miniasm`:

> ```
> minimap2 -c -x??? asm1.fa asm2.fa > out.paf
> minidot out.paf > out.eps && epstopdf out.eps
>
> ```

What does `-x???` stand for? And how do you visualize PAF files? Could you think of a way to use IGV to visualize PAF files?

# Answer
> 5 votes
As described in the minimap2 manual, `-x STR` is used to specify a preset option like `-x map-pb` for `PacBio/Oxford Nanopore read to reference mapping (-Hk19)`.

In general, I use the `-a` flag to output alignments in sam format that can then be sorted and converted to bam format, and viewed in any genome browser. The support for PAF format when it comes to visualization appears to be limited. A few options that I have found, but not tested are:

# Answer
> 4 votes
Depending on what level of information you need, it is not difficult to just homebrew a plot in your plotting package of choice, as PAF has relevant coordinate information. Unless you are interested in some more complex features of a plotter such as IGV, I think that this is a good flexible solution.

For example, this R code makes a dot plot of the synteny between the reference and the query:

```
library(ggplot2)
paf = read.table("out.paf", sep="\t")
ggplot(data=paf, aes(x=V3, xend=V4, y=V8, yend=V9, color=factor(V1))) + geom_segment() + labs(x="Reference coordinate", y="Query coordinate")

```

The plot looks like this:

It is effectively a rather minimalist version of what pafr (as suggested by vkkodali) or minidot is doing.

(h/t Zev K.)

---
Tags: minimap2
---


title: Create identicon of a DNA sequence

text:
# Question
Title: Create identicon of a DNA sequence
How can I create an identicon (a visual representation) of a DNA sequence?

# Answer
Sequenticon is a Python library for generating identicons for DNA sequences. For instance, the sequence `ATGGTGCA` gets converted into the following icon:

A web interface is also provided at EGF CUBA (Collection of Useful Biological Apps): Render Sequenticons

*Disclaimer: I'm the current maintainer of Sequenticon*
> 5 votes

---
Tags: visualization, egf
---


# Question
Title: How to calculate RNA copies in qPCR using R?
I have these Ct data from qPCR for a series of samples A-H, I want to quantify the RNA copies of my samples using R:

```
Sample_ID   Replicate1_Ct Replicate2_Ct
A           22.92         24.21
B           24.29         22.16
C           24.29         22.16
D           27.91         23.76
E           21.44         22.81
F           16.50         18.31
G           17.79         16.56
H           17.34         19.09

```

In addition, I have measured a standard (STD) curve of copy number of the RNA in question, that I want to use as a reference:

```
STD_Ct      STD_copies      Log10(STD_copies)
33          4.00E+03        3
28          4.00E+04        4
24          4.00E+05        5
20          4.00E+06        6
15          4.00E+07        7
12          4.00E+08        8

```

It seems

This the original data form:

# Answer
> 1 votes
Now that we've straightened out the data, it seems like the simplest solution is to just linearly interpolate the copy number from your standard curve.

When I read your standard curve (`std`) and your sample data (`ct`) into R, I see that the log10(Ct) fits the copy number fairly well, whereas the linear is a bit harder:

So we can fit a model to predict copy number based on Ct on the STD data:

```
> fit = lm(std$Log10.STD_copies. ~ std$STD_Ct)
> fit

Call:
lm(formula = std$Log10.STD_copies. ~ std$STD_Ct)

Coefficients:
(Intercept)   std$STD_Ct
    10.6847      -0.2357

# we plug those numbers into the linear equation,
# taking the mean of the two replicates for a best estimate.

> rowMeans(ct[,2:3]) * -.2357 + 10.6847
[1] 5.130429 5.210567 5.210567 4.595390 5.469837 6.582341
[7] 6.636553 6.391424

```

So those are log10(copy number) estimates, of a sort. You can get much more elaborate, of course.

---
Tags: r, statistics, qpcr
---


title: Creating a tab delimited column

text:
# Question
Title: Creating a tab delimited column
I have a blast file produced. I executed a `blast(x)` command outputting both `"qeseqid"` and `"sseqid"`:

```
QRv313_NP342_d0_h2_l9    YN13213
QRv313_NP9080_d0_h1_l1   YN5345
QRv313_NP123_d0_h1_l7    YN756
QRv313_NP123_d0_h1_l113  YN9768
QRv313_NP654_d0_h2_l6    YN432
QRv313_NP8_d0_h1_l1      YN3242
QRv313_NP756_d0_h1_l2    YN85686

```

I have written a command in `nano` within `command-line` to obtain the following desired output:

```
NP342    YN13213
NP9080   YN5345
NP123    YN756
NP123    YN9768
NP654    YN432
NP8_d0   YN3242
NP756    YN85686

```

I have written a `nano` script to provide me a tab delimited column of my query and subject id. I am just having trouble moving forward from here. I am unsure as to how I would modify my script to provide me with my desired output.

```
import sys
file_object = open(sys.argv[1])

for my_data in file_object:

  list =  my_data.split("\t")

  print (list [0], list [1])

```

Is there a way t alter my command so I can receive the desired output?

Any suggestions would be kindly appreciated!

# Answer
> 3 votes
as far as I understand your question, you want to reformat your output file that consists of two columns:

```
QRv313_NP342_d0_h2_l9    YN13213

```

to extract the second "field" from your first column, where fields are delimited by `_`:

```
NP342    YN13213

```

In general, your code is close. You could just split again by `_`:

```
import sys
file_object = open(sys.argv[1])

for my_data in file_object:
  list =  my_data.strip().split("\t")
  print (list[0].split("_")[1], list[1], sep="\t")

```

Now when I teach students, I usually tell them to always expect the unexpected, so you might also want to add some minimal error checking, so here's a version that checks some minimal assumptions:

```

from sys import argv
from os import path

usage = f"""
usage: {argv[0]} <blast_output>
"""

def main(in_file):
    with open(argv[1], "r") as in_handle:
        for lno, line in enumerate(in_handle):
            line = line.strip()
            cols = line.split("\t")
            if len(cols) != 2:
                raise ValueError(
                    f"Error in line {lno}: Expected two columns,"
                    " found {len(cols)}. Offending line: {cols}"
                )
            query = cols[0].split("_")
            if len(query) != 5:
                raise ValueError(
                    f"Error in line {lno}: Expected five fields"
                    "in query string, found {len(query)}. Offending line {query}"
                )
            print(query[1], cols[1], sep="\t")

if __name__ == "__main__":
    if len(argv) != 2:
        print(usage)
        raise ValueError("Illegal number of arguments")
    if not path.exists(argv[1]):
        raise FileNotFoundError(f"File {argv[1]} could not be found. Check spelling")

    main(argv[1])

```

And this is what it produces:

```
>python3 ex.py in.txt
NP342   YN13213
NP9080  YN5345
NP123   YN756
NP123   YN9768
NP654   YN432
NP8     YN3242
NP756   YN85686

```

---
Tags: python, file-formats
---


# Question
Title: What is a sensible number of gene/observations to explain PCA variance?
I am working with a set of RNASeq dataset. I have about 4000 observations (genes) on 20 samples and plotting a PCA I found the clustering doesn't vary much when I use different number of genes, but the % variance of PC1 and PC2 does. I have ~1300 DEG. But how do you select a sensible cutoff for the number of genes/observations to explain PCA variance? Would it be correct to select the number of genes that are DEG which account for \>50% of the variance between PC1 and 2?

Few examples:

* All genes: PC1(34%), PC2(14%)
* Top1000: PC1(45%), PC2(17%)
* Top500: PC1(49%), PC2(19%)
* Top50: PC1(55%), PC2(24%)
* Top 5: PC1(75%), PC2(24%)

# Answer
> 3 votes
The general goal of PCA in RNA-seq can be stated as, "I'd like a low-dimension representation of my data to allow easy assessment of the gross structure of my samples, specifically for assessing missing batch effects, samples swaps, etc." Ideally, this low-dimension representation can fit into a plot or two. In other words, we'd like the first 2 (maybe 3) PCs to explain a reasonable amount of variation. A few things to note from this alone:

* The number of differentially expressed genes were not mentioned. They have no relevance here.
* There was no mention made of the treatment groups clustering together. Seeing this is **NOT** a goal of PCA.

As you have noticed, the more input rows (genes, assuming you haven't already transposed it) the less variation explained by the first 2 principal components (PCs). On the flip side, the more rows you input the better PCA is describing the structure of your entire dataset.

Question: Do I really care about the gross structure of the entire dataset?

Answer: No, I just want to see problems that I might need to account for.

So you really don't have to use all of the data in the PCA, just the bits that will allow you to see any issues. In general that'll be the few hundred most-variable genes (after some reasonable transformation so you're not just looking at the most highly expressed genes!). **The exact number shouldn't really matter very much, something in the range of 300-500 most variable genes should suffice.** A common practice is to just randomly pick a nice round number somewhere in or around this range. This is not a number that needs optimization unless perhaps you have a very large number of samples and therefore PC1 and PC2 only account for a small (say 10-20%) percentage of the variation, in which case you're better off looking at either a number of difference principal components or using a different dimension reduction technique, such as UMAP.

---
Tags: rna-seq, deseq2, pca
---


# Question
Title: Is is possible to predict ncRNAs from sequence and homology alone?
I'm working with a set of homologous genes (let's call it gene A) from several bacterial species.

I know (from previously published research) that in gene B (a close paralogue of gene A), there is a non-coding, antisense RNA. This antisense RNA regulates gene B by forming a duplex with the gene B mRNA.

I have noticed that, sometimes, an antisense RNA would be predicted for gene A in the same location in the IMG/MER database. This is not predicted for all of the available gene A sequences, only some of them. IMG/MER uses INFERNAL to predict RNAs, I think.

My question is - is there any way to assess how real this is without wet lab data? As it is, it could be that there's no actual antisense RNA in gene A, and this is a spurious prediction due to homology. Or it could be that the antisense RNA is present in all of the gene A sequences, and the RNA prediction software simply missed some of them. Is there a computational way of testing which of these options (if any) is more likely? Thank you for your time!

# Answer
> 1 votes
Is there any conserved secondary structure in this RNA? One could align all known instances of this RNA and see if there is any significant covariation according to R-scape. That could serve as a computational signal that the RNA is real (although wet lab experiments would still be needed).

Another question is whether this RNA is already in Rfam? If so, then the job of making an alignment of known instances is easier.

---
Tags: rna, software-recommendation, non-coding
---


# Question
Title: Detecting whether genome is genetically engineered?
Is there a methodology for detecting whether a genome has been tampered with by a biologist? For instance, say there has been an ebola outbreak such as in "The Andromeda Strain". Is there a way a researcher can tell by just looking at the ebola genetic data and determine it is a genetically engineered strain?

Some ideas I have:

1. organism is abnormally close to unrelated organisms, i.e. some of the genetic code from the other organisms has been inserted into the first
2. subsequences are very different from the rest of the genome. If there is a consistent rate of mutation applied to the entire genome, then a subsequence that is abnormally noise free could be inserted.
3. combination of 1 and 2: finding clear instances of unrelated organism genomes inserted into the genome being examined.

Additionally, are there any publicly accessible reference genomes that are known to be genetically engineered, that I could use as examples to test such ideas?

UPDATE: the wikipedia article on detecting genetically modified organism is not quite what I am looking for, since it claims current methods cannot detect unknown genetic modification.

https://en.wikipedia.org/wiki/Detection\_of\_genetically\_modified\_organisms

> Currently, it is highly unlikely that the presence of unexpected or even unknown GMOs will be detected, since either the DNA sequence of the transgene or its product, the protein, must be known for detection.

# Answer
> 2 votes
The definition of a GMO for example UK GMO regulatory bodies (HSE) is strictly the artificial introduction of foreign genetic material using targeted genetic engineering, such as splicing, genome integration etc... Attenuation for example by serial passage the basis of the yellow fever vaccine, does not count as a GMO.

The simple answer is look at phylogenetic trees between the GMO strain against the total distribution of isolates. If large regions have been spliced the 'silent mutation rate' looks very different,.because the GMO parent is from a very different environment to the recipient. The switch in replication rate affects the numbers of 'silent mutation', which is fast in RNA viruses.

If the GMO is a point mutation, ie limited amino acid change then I've actually done this with dengue virus and it is fairly easy because the novel mutation make the GMO 'stick out'.of the phylogeny against all natural mutation patterns. At this point the outlier is tracked back to Genbank where the relavant publication defines the engineered virus.

If you are inferring a conspiracy theory that is a different process and a different discussion, particularly if it is Ebolavirus which is a Category 4 pathogen. The investigation is again within the phylogeny, but in this case the answer is always the isolate is a natural zoonose from a sylvatic (wild) transmission cycle which has spilled over into humans.

# Answer
> 0 votes
If you need to engineer a trait, then it has not evolved, so it probably does not confer a selectional advantage to the organism.

So a fingerprint of an engineered organism is that the genetic difference between it and its relatives does not have a selectional explanation.

---
Tags: phylogenetics, phylogeny, covid-19, genetics, sars-cov-2
---


# Question
Title: Picking number of samples for clustering
I have expression data for 100 samples and clinical data for 50 among them. I want to carryout k-means clustering using expression data and then do a survival analysis (Kaplan-Meier analysis) using the labels obtained from k-means.

Can I used the expression data of all 100 samples to carryout clustering or should I use only 50 samples for clustering?

# Answer
> 1 votes
In case there is no systematic difference between the samples with or without the clinical data (i.e. one set is microarray and the other is RNA-seq), I would use all of the samples. And if the cluster labels obtained with "clustering 50 samples" and those obtained with "clustering 100 samples" does not differ too much, then doing either of two would not have much effect anyway. If you would like to compare cluster labels in that way and assuming you would want to test a number of `k`s, it would be better to compare each combination of `k`s and number of samples.

---
Tags: clustering
---


# Question
Title: Why don't the "child" peptide masses sum to the parent in this De novo sequencing example?
In that example, 129+216+327 does not equal the parent's mass of 492. My understanding is that during MS, the protein is broken down into peptide masses of varying length of approx. 5 to 20 amino acids. They also lose ions such as water or ammonia like in the example., the experimental spectra I thought would show the masses of these "children" that can be summed up to the parent (if you also add all the ammonia and water ions that are broken off). But that is not the case. Where am I wrong in my understanding?

# Answer
> 1 votes
In a peptide tandem mass spectrometry experiment, each fragment (child) ion is typically the result of breaking a single bond in the peptide backbone yield b and y ions. Because each ion fragment ion is the result of one fragmentation event, we don't expect the fragment ion masses to sum to the precursor mass.

If your peptide was fragmented well, you should get a ladder of b and y ions where the differences in masses between ions in the b and y ion series correspond to the masses of amino acids (+/- any modifications and neutral losses, such as the NH3 and H2O mentioned above).

A great resource for understanding these concepts is this tutorial:
https://www.ionsource.com/tutorial/DeNovo/introduction.htm

---
Tags: mass-spectrometry
---


# Question
Title: how to train a gene dataset with a nearest shrunken centroid classifier?
I have a data file named "geneexp.csv".

the data contains information about gene expression of three different cell types (CD4 and CD8, CD19) I want to classify cells by performing the nearest shrunken centroid classification of training data in which the threshold is chosen by cross-validation. I split the data (70% train and 30% test).

```
data = read.csv("geneexp.csv")

splitData <- function(data, trainRate) {
  n <- dim(data)[1]
  idxs <- sample(1:n, floor(trainRate*n))
  train <- data[idxs,]
  test <- data[-idxs,]
  return (list(train = train, test = test))
}

split <- splitData(data, .7)
train <- split$train
test  <- split$test

```

then with the use of `pamr` package I tried to buid the following model and plot :

```
y <- train[[ncol(train)]]
x <- t(train[,-ncol(train)])
mydata <- list(
  x = x,
  y = as.factor(as.factor(y)),
  geneid = as.character(1:nrow(x)),
  genenames = rownames(x)
)

# Training and cross-validating threshold
model <- pamr.train(mydata)
cvmodel <- pamr.cv(model, mydata)
pamr.plotcv(cvmodel)

```

but I can't make it work. I get the following error:

> Error in `contrasts<-`(`*tmp*`, value = contr.funs\[1 + isOF\[nn\]\]) :
> contrasts can be applied only to factors with 2 or more levels

I have already transfer y to the factors. Can you help me? How can I fix it?

# Answer
> 0 votes
It's a bit different from usual R, but if you check the help:

> data: The input data. A list with components: x- an expression genes in the rows, samples in the columns)

So in your case, you need to transpose the matrix:

```
mydata <- list(
  x = t(x),
  y = as.factor(y),
  geneid = as.character(1:nrow(x)),
  genenames = rownames(x)
)

```

As I don't have your amazing data, I can only use iris below:

```
pamr.train(list(x=t(iris[,1:4]),y=iris[,5]))

pamr.train(data = list(x = t(iris[, 1:4]), y = iris[, 5]))
   threshold nonzero errors
1   0.000    4       6
2   0.841    4       7
3   1.682    4       10
4   2.523    4       11
5   3.364    4       13
6   4.205    4       18
7   5.046    3       22
8   5.887    3       23

```

---
Tags: r, gene, genome, machine-learning, data-preprocessing
---


title: chipseeker annotation issue

text:
# Question
Title: chipseeker annotation issue
```
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
library(clusterProfiler)
library(org.Hs.eg.db)

peaks <- readPeakFile("MACS_folder/pol2_chip/Pol_2pausing/pol2_diff_peak.bed")

print(peaks)

#peakAnno <- annotatePeak("merge/blast.bed", tssRegion=c(-3000, 3000),
#                        TxDb=txdb, annoDb="org.Hs.eg.db")

peakAnno <- annotatePeak(peaks, tssRegion=c(-1000, 1000),
                         TxDb=txdb, annoDb="org.Hs.eg.db")

```

So this is basically what im doing to annotate the peaks.

Now the issue comes with annotation for example I see TLR4 annotated to this

```
118082896   118083202   307 *   shp1_R1_peak_679    Exon (ENST00000413759.1/ENST00000413759.1, exon 1 of 1) 9   117708197   117708979   783 1   7099    ENST00000490685.1   374699  ENSG00000136869 TLR4

```

Now The peak is 679 as you can see from the IGV plot as well as annotated file. But when i see in IGV browser TLR4 is far away from the peak where i get the annotation.

So now is it an issue with the code im using or due to the inherent nature of the data.

Any help or suggestion would be really appreciated.

# Answer
I had a similar case just recently.

The `tssRegion` argument ChIPseeker is used to define what is considered a promoter, which is frequently the question at hand. However, ChIPseeker will not perform any kind of filtering on the annotated peak file.

In the results there is a column `annotation` that contains information what kind of genomic region the peak is in. If it is within your `tssRegion` window it will be labeled as a promoter.
In your case the annotation is `1st Exon = Exon (ENST00000413759.1/ENST00000413759.1, exon 1 of 1) 9`. This tells you that the peak is not part of the promoter region you defined and located within the first exon of the indicated transcript. The peak is really far away from TLR4 (374699 bp) but I think all other sequences nearby are pseudogenes and this is why TLR4 is selected as the closest gene.

You can get an overview over the annotations using `plotAnnoPie(peakAnno)`

tl;dr: you need to filter the results of the ChIPseeker annotation to get the genomic regions you want.
> 2 votes

---
Tags: ngs, chip-seq
---


title: shiny code with seurat object

text:
# Question
Title: shiny code with seurat object
I am trying to create a shiny board for my scRNA seq data. I have saved my seurat objects adn I use following code as below:

```
library(shiny)
library(shinythemes)
library(Seurat)
library(uwot)

set_1 <- readRDS("/temp2/data/mydataB10.rds")
set_2 <- readRDS("/temp2/data/mydataB5.rds")

ui <- fluidPage(

    # Application title
    titlePanel("scRNAseq data browser"),

    fluidRow(
        column(4,

               sidebarPanel(
                   radioButtons("dataset", label = h3("Dataset"),
                                choices = list("sample1" = "set_1", "sample2" = "set_2"),
                                selected = "set_1"),
                   helpText("Enter gene names in CAPS.  Gene names must be exact."),
                   textInput("gene2", label = "Gene Name", value = "GAPDH")

               )
        ),
        column(8,
               column(6,
                      mainPanel(
                          plotOutput("dimPlot1", width = "150%"),
                          plotOutput("genePlot2", width = "150%")
                      )
               )
        )
    )
)

server <- function(input, output, session) {

    datasetInput <- reactive({
        df <- (input$dataset)
            })

    output$dimPlot1 <- renderPlot({
        DimPlot(datasetInput(), reduction = "umap")
    })

    output$genePlot2 <- renderPlot({
        FeaturePlot(datasetInput(), features = (input$gene2), reduction = "umap")
    })

}

shinyApp(ui, server)

```

That error that I get on the section where image needs to be displayed is "Subscript out of bound".

Following dataset can be downloaded and used as an example `curl -Lo pbmc_10k_v3.rds https://www.dropbox.com/s/3f3p5nxrn5b3y4y/pbmc_10k_v3.rds?dl=1`

Can anyone help in fixing this shiny script so that it works.

Many thanks.

# Answer
> 3 votes
Not easy to reproduce your error but here is my guess. This is a string, you can check it:

```
input$dataset
class(input$dataset)

```

This is a simple solution. If you want to call out the dataset depending on input, use a list:

```
set_1 <- readRDS("/temp2/data/mydataB10.rds")
set_2 <- readRDS("/temp2/data/mydataB5.rds")
data_list = list(set_1=set_1,set_2=set_2)

```

Then in the server :

```
datasetInput <- reactive({
        df <- data_list[[input$dataset]]
            })

```

---
Tags: r, scrnaseq, seurat
---


# Question
Title: Snakemake wildcard constraints reject too much
I have the following two rules:

```
rule extract:
    input: "data/{dir}/{file}.gz"
    output: "data/{dir}/{file,^.*(?<!\.gz)$}"
    conda: "config/conda-extract-env.yml"
    shell: "cd 'data/{wildcards.dir}'; gunzip -k {wildcards.file}.gz"

rule extract_dot:
    input: "data/{dir}/wtdbg2.3.dot.gz"
    output: "data/{dir}/wtdbg2.3.dot"
    conda: "config/conda-extract-env.yml"
    shell: "cd 'data/{wildcards.dir}'; gunzip -k wtdbg2.3.dot.gz"

```

The second one is actually a hack because the first does not work on files named `wtdbg2.3.dot`. I have also tried giving the constraint using the `wildcard_constraints` member of the rule, and prefixed the string with an `r` or not.

Whatever I do, as soon as I uncomment the `extract_dot` rule, snakemake complains about the missing input file `data/C.elegans/wtdbg2.3.dot` in another rule. Note that there is a rule with output `data/C.elegans/wtdbg2.3.dot.gz`, and when changing this rule to output `data/C.elegans/wtdbg2.3.dot` instead, everything works as well.

So I conclude that the problem is in the regex. I have tested the regex on regex101.com, and it accepts `wtdbg2.3.dot`. But snakemake seems not to. Does anyone know if there is anything special about the snakemake regexes? Should they not be normal python regexes?

# Answer
> 2 votes
> Does anyone know if there is anything special about the snakemake regexes? Should they not be normal python regexes?

No, nothing special. The Snakemake documentation just links to the Python re module documentation.

I think this is something strange in how the anchors `^` and `$` behave when matching wildcards, and not anything else in your particular regular expression. Check out the pattern-matching with this dummy Snakefile:

```
rule all:
    input: ["data/subdir/example", "data/example", "example"]

rule anchored_example_of_subdir:
    """Can't add either ^ or $ or this will fail"""
    output: "data/{subdir,subdir}/example"

rule anchored_example_in_subdir:
    """Can't add ^ or this wil fail"""
    output: "data/{file,example$}"

rule anchored_example:
    """This works with both ^ and $"""
output: "{file,^example$}"

```

(There are no inputs or actions so it's just a matter of what filenames match what rules, and the regular expressions do nothing but match exact strings.) This "works" as is (`snakemake -nrp` shows what you expect) but if you add in anchors as described in the docstrings above the upper two rules no longer match the expected filenames. It's as if (and maybe is?) under the hood Snakemake is trying to apply the expression to the entire path.

Can you just drop your `^` and `$`? The snakemake documentation implies that it's using the equivalent of `re.match` rather than `re.search` already, with wording like "Here we could restrict the wildcard `dataset` to consist of digits only using `\d+` as the corresponding regular expression" so I don't think it's necessary anyway.

---

Edit: I think I found the issue in the snakemake source. It munges your constraints together into one larger regular expression for the entire file path, so `^` and `$` get embedded *inside* the resulting expression. This looks like something that should at least be clarified in the documentation.

---
Tags: snakemake
---


# Question
Title: Virus genome: protein inside another protein
The genome sequence for the GenBank virus JX869059.2 (a human betacoronavirus) contains, among others, two proteins:

* N protein (range 28566..29807) and
* Orf8b (range 28762..29100).

The interval for Orf8b is contained in the interval for N.

How are these proteins translated? Is the orf8 protein created with post-translational processing after N is translated, or does the ribosome translate the two proteins separately?

# Answer
> 3 votes
Note that the protein sequence of ORF8b is not contained within the protein sequence of protein N. That means that they are translated from different frames (though they are in the same orientation). Indeed, if we just subtract the start coordinates we do not find a multiple of 3.

Thus, they cannot be results of the same translational event. This means that the smaller protein could result from translational slippage (as suggested by Matteo Ferla), "leaky scanning" of the translation machines for start codons, or it could simply result from a different ribosome assembly site, recruited via the viral mechanism of an internal ribosome entry site, which allows a ribosome to be recruited in a mRNA cap-independent fashion (though this may be less likely).

For some mechanistic details of how one such ORF is translated, see this paper.

It looks like such things are sometimes called "subgenomic mRNAs". For a paper on such nested genes in coronavirus, including I believe Orf8b, see here.

# Answer
> 2 votes
Essentially @MaximilianPress hit the nail on the head. sgRNA is a major feature of SARS-CoV-2.

This is where the 5' UTR is spliced at numerous locations in the genome, such that although it is a single stranded positive sense RNA genome it behaves like a segmented viral genome because it is translated as a series of sub-genomic RNA molecules. So influenza is a classic example of a segmented genome. SARS-CoV-2 is more slightly complicated because it is using frameshifts to achieve to achieve the different sgRNAs.

However it isn't a particularly strong example of frameshifts because virus such as HBV (hep B) are pass masters, where two different proteins are encoded within a single locus.

My personal view is that slippage isn't operating here, unlike say HIV, but that it is a function of the sgRNA construction, thus N protein occurs on one sgeRNA and ORF8b occurs on a second sgRNA and sg signal is within the N protein.

If SARS-CoV-2 is your thing then stuff like this is a big player in interpreting bioinformatic output.

---
Tags: ngs, covid-19, sars-cov-2, coronavirus
---


# Question
Title: How do I rename output folder created using WDL's glob function, with cromwell execution engine?
I made a wdl (genomics pipeline language by the Broad Institute) task for fastqc that outputs the zip files, of which there will be as many as reads that the sample I'm working on at that moment has. We create the output using glob:

```
  output {
    Array[File] fastqcArray = glob("*.zip")
  }

```

I would like to give the folder a name that is representative of what it contains, for example `fastqc`, rather than glob-xxxxxxx, so that it's clear what it is in my task's output folder. Is this possible, and if so, how can I achieve it?

# Answer
Ok, thanks to Ruben Vorderman on the openwdl github page, I found a solution that does not involve renaming the glob folder (which I believe is not possible?). The only way is to write out the names of the files and pass them explicitly to the Array\[File\] output. This is how I did it (simplified), in case it helps anyone else in the future:

```
task fastqc_task {
  File forwardReads
  File reverseReads
  String fwdName = sub(basename(forwardReads),".f.*q.*$","")
  String revName = sub(basename(reverseReads),".f.*q.*$","")
  command {
    mkdir -p ./qc/fastqc_zip
    fastqc ${forwardReads} ${reverseReads} --outdir ./qc/fastqc_zip
  }
  output {
    Array[File] fastqcArray = ["./qc/fastqc_zip/${fwdName}_fastqc.zip", "./qc/fastqc_zip/${revName}_fastqc.zip"]
  }
}

```
> 2 votes

---
Tags: wdl
---


title: Mapping protein refseq to Gene ID

text:
# Question
Title: Mapping protein refseq to Gene ID
I have a protein refseq (eg, NP\_000029). How can I get the corresponding gene ID (ag, APC) from NCBI using an R package?

# Answer
> 2 votes
You can do this with the reutils package, which provides an API to NCBI's E-utilities. Here's an example for your specific question:

```
install.packages("reutils")
library(reutils)

# Get universal identifier
uid <- esearch("NP_000029", db = "gene")
# Fetch summary
sm <- esummary(uid, db = "gene")
# Extract specific tag
sm$xmlValue("//Name")

```

---
Tags: gene, genome, phylogenetics, proteomics, refseq
---


# Question
Title: How do I change the identity of a sample from spatial 'Visium' data preprocessing in Seurat v3?
I'm wanting to create a merged object in Seurat using 2 10x Visium 'slices'. However, when I create an object, Seurat assigns an identity "SeuratProject" to the objects (by default I'm guessing; see attached plot). This will be problematic when I attempt to merge the objects as there will be unintended duplicates of the objects.

Here's what I'm attempting to use to create the 10x Spatial Object:

```
brain1 <- Load10X_Spatial(data.dir = "~/xxx", filename = "filtered_feature_bc_matrix.h5", assay = "Spatial", slice = "slice1", filter.matrix = TRUE, to.upper = FALSE)

```

I'd need the ability to assign unique ident to each of these (and any subsequent) slices/data.

Thank you!!

# Answer
> 1 votes
You can re-assign `orig.ident` after creating the object. For example:

```
brain1 <-Load10X_Spatial(data.dir = "~/xxx", filename = "filtered_feature_bc_matrix.h5", assay = "Spatial", slice = "slice1", filter.matrix = TRUE, to.upper = FALSE)
brain1$orig.ident <- "Visium"

```

---
Tags: seurat, single-cell
---


# Question
Title: Summing up one column entries for each sample
I have a dataframe with 27 samples and for each of them I have different number of values

I want to sum up the values for each sample, leaving each sample with a single value

```
> dput(d[1:10,c(1,8)])
structure(list(SampleID = c("LP6008460-DNA_A04", "LP6008460-DNA_A04",
"s30", "s30", "s67", "s15", "s67", "s67", "s67", "s30"), X.1 = c(8.17407e-05,
1.40067e-06, 3.86017e-05, 9.50247e-05, 7.08333e-07, 3.19133e-05,
3.1602e-05, 5.00513e-05, 3.29613e-05, 1.003e-06)), row.names = c(NA,
10L), class = "data.frame")
>

```

I started to do that manually but that was a pain and inaccurate

# Answer
> 2 votes
It's worth learning how to use `dplyr` for this sort of thing:

```
dat %>%
    dplyr::group_by(SampleID) %>%
    dplyr::summarise(sample_mean = mean(X.1))

```

```
# A tibble: 4 x 2
  SampleID          sample_mean
  <chr>                   <dbl>
1 LP6008460-DNA_A04   0.0000416
2 s15                 0.0000319
3 s30                 0.0000449
4 s67                 0.0000288

```

---
Tags: r
---


title: GC normalization ATAC seq data

text:
# Question
Title: GC normalization ATAC seq data
Is it really required to do GC normalization ATAC seq data ?

One of the paper where they have ATAC seq data did the normaisation of GC bias after peak calling. When i looked for the library and the reason why GC bias i found this

```
Could it be true that genes with higher GC content are higher expressed? It has been suggested that genes that are either extremely high or extremely low expressed are under some form of selection leading to “extreme” GC content. What CQN does, is making the effect of GC content comparable across samples, and we show in [1] that this leads to improved inference. It also flattens the effect of GC content on gene expression, but we believe this is better than having the effect of GC content depend on the sample.
https://www.bioconductor.org/packages/release/bioc/vignettes/cqn/inst/doc/cqn.pdf

```

Now that scenario is for RNA seq would it be logical to do the same for ATAC seq?

# Answer
The only reason to normalize for GC content in RNA-seq is if it differs notably between samples/groups. If that's not the case and you aren't trying to compare genes withing samples then you have no reason to try to account for GC content. The same goes for ATAC-seq, though there the danger with trying to normalize for GC content is that you then mask changes in global accessibility.

In short, normalize for things that should cause problems with your analysis and nothing else.
> 2 votes

---
Tags: atac-seq
---


# Question
Title: What is the best way to map/align my reads on a given genome?
I am frequently using a ballgown package for my rnaseq analysis, but recently I have had a new task to have my reads mapped on two different genomes to understand the level of alignment between the two, therefore I wonder what is the best tool to do this? I have single-end reads. What I need is to have a percentage of uniquely mapped fragments, mapped to multiple loci and, say, unmapped. Can `bowtie` do the job? Is there a standard explanation procedure explained elsewhere? Thanks.

# Answer
This is a pretty straightforward analysis, assuming that you want to measure something like "how many reads align to genome 1, and how many reads align to genome 2?".

You should be able to then use samtools to compute statistics on your alignments.
> 0 votes

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


# Question
Title: How to find Genotype, from list of mutations:lineage list?
I'm very new to python and I'm getting along nicely (I'd like to believe); however, I must be missing something here. I'm looking to read each and every file and compare each line against one another.

The aim of this is to see if any mutations are shared in any of the strains - to determine the genotype.

The inner-workings of the files are single spaced headers, as those you would find on GISAID.

`os.chdir("C:)`

with open('A D614G.txt', 'r') as file2:
with open('A E780Q.txt', 'r') as file3:
with open('A G476S.txt', 'r') as file4:
with open('A L18F.txt', 'r') as file5:
with open('A N439K.txt', 'r') as file6:
with open('A S477N.txt', 'r') as file7:
with open('A T478I.txt', 'r') as file8:
with open('A V483A.txt', 'r') as file9:

files = [file1, file2, file3,file4, file5, file6, file7, file8, file9]

it = itertools.permutations((files), len(files))

for x in it:
    print(x)

os.chdir("C:)
with open('tester.txt', 'wt') as file_out:
    for x in it:
        file_out.write()```

1 I recieve the following error on a very long repeated stretch of text

```<_io.TextIOWrapper name='A G476S.txt' mode='r' encoding='cp1252'>```

</code></pre>

# Answer
```
import re
import os
import pandas as pd
import logging
import sys
import csv
import itertools

mutations = ['A222V', 'D614G', 'E484Q', 'E780Q', 'G476S', 'L18F', 'N439K',
             'S477', 'S477N', 'T478I', 'V483A']

combinations = []
for M in range(1, len(mutations)+1):
    for subset in itertools.combinations(mutations, M):
        combinations.append(subset)

        combinations = ['_'.join(sorted(x)) for x in combinations]
        combinations = [x.split('_') for x in list(set(combinations))]

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

combination_labels = []
combination_counts = []

for lineage in lineages:

    df = pd.read_csv('Results/' + lineage).dropna()

    for combination in combinations:

            combination_df = df[list(df)]

            for mutation in combination:
                combination_df = combination_df[combination_df[mutation] == 1]

        #print(combination_df.shape[0])

            combination_labels.append('_'.join(combination))
            combination_counts.append(combination_df.shape[0])

            out_df = pd.DataFrame({'combination':combination_labels,
                               'count':combination_counts})
            out_df['percentage'] = (out_df['count'] / df.shape[0]) * 100
            out_df = out_df.sort_values('percentage', ascending = False)

            out_df.to_csv('Results_2/' + lineage.replace(".csv", "") + '_3.csv',
                     header = True,
                     index = False)

```
> 0 votes

---
Tags: genomics, covid-19, sars-cov-2, python
---


# Question
Title: How to get % similarity between strains and mutation files
I’m very new to python, and having some difficulty getting hang of some more complicated things

I have multiple files which look like so:

> hCoV-19/Singapore/4/2020|EPI\_ISL\_410535|2020-02-03

> hCoV-19/USA/WA13-UW9/2020|EPI\_ISL\_413601|2020-03-02

> hCoV-19/USA/WA-UW142/2020|EPI\_ISL\_416680|2020-03-11

**Please be aware that the lines above are meant to be one file**

I want to extract the EPI\_ISL\_000000 for an easy comparison among files.

Could someone please advise on:

1. A programme to extract this data into new files (There’s many lines in each file - 1000+)
2. A programme to then give a % comparison between two or more files - comparing all lines in one file against all lines in a second+ file

# Answer
```
   left_lineagelist = [x.split('_')[-1].split('|')[0]
                          for x in left_lineagelist]
        right_lineagelist = set([x.split('_')[-1].split('|')[0]
                          for x in right_lineagelist])

```

Allows for extraction of 6 digit EPI, provided the file has had sequences removed prior; as such:

```
for line in lines:
    if line[0] == '>':
        print(line[1:])

```
> 0 votes

---
Tags: phylogenetics, genomics, covid-19, sars-cov-2
---


# Question
Title: How to count the kmer occurrence in FASTA file considering overlapping and reverse complement?
I am using count\_overlap() for counting the kmers from Biopython .

Does it take the reverse complement of kmer into account?

I need to count reverse complement too

# Answer
There's a lot detail missing here, but at a guess could you try,

```
from Bio import SeqIO
from Bio.Seq import Seq
for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    reversec = record.seq.reverse_complement())
    revcount = reversec.count_overlap('AT')
    print (revcount + count_overlap('AT'))

```

Essentially `count_overlap` only counts the forward strand if this is the Biopython method. Not debugged the above.
> 1 votes

# Answer
Here's a C++ kmer counter I wrote for my own purposes. It can count canonical and reverse complement kmers from FASTA files up to length 31. Run `kmer-counter --help` or read the README for count options.

For longer kmers, you'd want to look at an approximation tool like Jellyfish.
> 1 votes

# Answer
I would NOT recommend using biopython for a task like this. First, it's slow, but also it's really not necessary. There are several other existing tools optimised for exactly this purpose.

My favourite is KMC. It's very fast, and by default, it reports canonical kmers (Parameter `-b` would turn this off). So, assuming you have KMC installed, this is how you get all the kmers from a range of occurances (1 - 10000).

```
# just making a directory for temprary files
mkdir tmp

# building a kmer database with k = 21,
# 16 cores and
# up to 64GB of memory.
# Counting occurances between 1 to 10000,
# `-fa` specifies that the input is fasta, you can also input other formats
# the database name will be kmer_counts
# & saving temporary data to tmp directory;
kmc -k21 -t16 -m64 -ci1 -cs10000 -fa my_seq.fasta kmer_counts tmp # run kmc

# extract kmers from the database in a text file
kmc_dump -ci1 -cx10000 kmer_counts kmer_k21.dump

```
> 0 votes

---
Tags: python, covid-19, sars-cov-2, k-mer, seqio
---


title: Receptor-ligand database

text:
# Question
Title: Receptor-ligand database
How can I find ligands of receptors (or *vice versa*)?

For example, I want to find the receptors of IL-8 (CXCL8).

# Answer
CellPhoneDB is a publicly available repository of *curated* receptors, ligands and their interactions.

The interactions can be searched or downloaded.
> 3 votes

# Answer
Have you tried STRING?

It's a publicly available protein-protein interaction database, so it won't be a complete picture.

Ingenuity Pathways Analysis (IPA) is really nice but is a paid service. Ask your institution if they have a license for one.
> 2 votes

# Answer
CellTalkDB is another option for a searchable ligand-receptor database. According to the webpage it is curated from other databases and the literature. Details are available in the publication.
> 0 votes

---
Tags: database, public-databases, protein-protein-interaction
---


# Question
Title: Metagenome simulation from a concatenated FASTA file
I am trying to simulate metagenome sequencing. So I will start with a file with a lot of concatenated genomes. From there, I would like to randomly extract 10,000 sequences of length 200bps.

I got help subtracting a random 200 sequence from a fasta file like this:

```
>JQ086372.1 Enterobacteria phage HK446, complete genome
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTCGCAGAAAGAGATTTTTCACCTTAAAGGCTTCGGATTCACCGGGCTTGTAGGCCTGTCAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

```

using this:

```
tail -n +2 file.fasta | tr -d '\n' > newfile
n=$(stat -c "%s" newfile)
r=$(shuf -i1-"$((n-200+1))" -n1)
newfile tail -c+"$r" | head -c200

```

And as an output I got: `GTTAAGGGCGCCGTTCTGCTTGTCTTTGCTGACATGTTTGAACACCGGACGGCACAAAGCGAAGTGCAGCTTTATGAGAATGCAGCCGCCGAACGCATGATGTTCATTCATCGCAACTGGCGCGGTAAATCTGAACCTGAGGAGGGCTCCTGATGGAACCTGGACGATTCAGGCACAGGGTAAAAATTCTCACCTTCACG`

But now I want to extract instead of 1 fragment of 200 characters; 10,0000 random fragments of 200 characters using the strategy above. For this, I am looking to exclude all those lines that begin with "\>" (p.e. \>JQ086372.1 Enterobacteria phage HK446, complete genome).The input file looks like this:

```
>JQ086372.1 Enterobacteria phage HK446, complete genome
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTCGCAGAAAGAGATTTTTCACCTTAAAGGCTTCGGATTCACCGGGCTTGTAGGCCTGTCAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

>NC_010391.1 Salmonella phage Fels-1, complete genome
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCTCCATATACTGTATAGACCATCTGAGCGTTTGAGTGCCCCATTTGTGAAGCAATAAAGTTCGGGTTCGCGCCAGCGGCAAGCGCCCAGCATGCACCGATTTTTTTAGTTTTAAGGTGTTGACATTAGGTATGTCGGGACTGGTATGCTTTCCTGTGTCGCAGCCCGGCGCGTCTCAATGCAGATTCCCATATCCTGTTCAGGGAAGGTGCAGCATAATGTGCTTTGGATTTACCTGAGCGATCGGTAATTTGCGGATTGAAGACAAAGGTAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

>NC_007577.1 Prochlorococcus marinus str. MIT 9312, complete genome
AATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTACAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC

```

So for my new output, I will end up with a file with 10,000 subsets of 200 characters that would look like this:

```
>1
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTTCAGTCTGCCAATATGCGTTCCGACAGCAGATAAAC

>2
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCTCCATATACTGTATAGACCATCTGCAGATA

>3
AATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAAT

```

Makes sense? Any help?

Thanks in advance

# Answer
> 1 votes
Given that you want to simulate short-read sequencing from reference genomes, there are read simulator programs that already exist. For example, see this answer.

---
Tags: fasta, ngs, file-formats, shell, subset
---


title: Segment cell organelles with pixellib

text:
# Question
Title: Segment cell organelles with pixellib
I've got some images of cell organelles and I really want to avoid labeling them by hand. The images all look something like the image below. Is there already an existing model specifially for cells? I tried it with existing models like `mask_rcnn_coco.h5` or `deeplabv3_xception65_ade20k.h5` but as they aren't built for cells the results aren't good.

My code is pretty basic right now:

```
from pixellib.instance import instance_segmentation

from PIL import Image

    def segmentImage():
        segment_image = instance_segmentation()

        segment_image.load_model("mask_rcnn_coco.h5")
        segment_image.segmentImage("cell.jpg", output_image_name="cell_seg.jpg", show_bboxes=True)

```

Cell example image:

Example segmentation result:

# Answer
I am the creator of PixelLib. PixelLib supports instance segmentation of 80 classes of objects using pre-trained coco model. Cell organelles is not part of the classes supported. The only solution is to train your own custom model. PixelLib supports custom training of a model using your own labelled dataset. Read this article I published on how to train a custom model with PixelLib.
> 2 votes

---
Tags: python, single-cell, tensorflow
---


# Question
Title: What does the 'suppressed' folder mean on the NCBI ftp genomes website?
I wish to understand the term `suppressed` in folder nomenclature on the NCBI ftp genomes website?

Example: here ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/Oscillatoriales\_cyanobacterium\_C42\_A2020\_001/all\_assembly\_versions/

Normally, I'd see folders with names like `representative` or `latest assembly`, which are pretty straightforward. But what does `suppressed` refer to in this case?

Thank you for your time!

# Answer
> 2 votes
Here you have the explanation:

> Suppressed genomes: Several hundred RefSeq bacterial genomes that did not pass assembly or annotation quality validation have been suppressed(RefSeq genome quality criteria). An FTP report file listing the most recent suppressions will be provided soon in the RefSeq FTP site.

---
Tags: phylogenetics, webservice
---


# Question
Title: Check if a string is ensembl gene id or not in python
I have many `scRNAseq` datasets saved as `hdf5` files (hierarchical data format\*), and when I open them keys sometimes are confused, so I need a way of determining from string values themselves whether the array contains `gene symbols` or `ensembl ids`. What would be a regular expression in python to check for it?

\*, A database format for quickly accessing large amounts of data, which is accessed using `keys`. The format is quicker than pickle/unpickle and can be viewed via a hdf5 viewer, a bit of info is here

# Answer
> 2 votes
Fixed it myself actually, not that hard. I found `regular expression` for `ensembl ids` here:

https://www.biostars.org/p/99142/

And used the `re` package in `python`:

```
import re

ensembl_regexp = 'ENS[A-Z]+[0-9]{11}'

is_ensembl = np.all([bool(re.match(ensembl_regexp, gene)) for gene in gene_symbols])
    if is_ensembl:
        gene_symbols = list(f['ensembl_ids'])

```

Where `gene_symbols` is a `numpy` array that I need to check whether it has `ensembl ids`, and if yes, switch to using `gene symbols` which are going to be under `ensembl_ids` key due to a mistake of someone who made the `hdf5` dataset.

# Answer
> 2 votes
The answer by terdon is probably good enough, but some commenters have noted it will allow something like `ENSABCDEFGHIJKL12345678901`.

It's helpful to reference the official documentation of Ensembl Stable Ids:

> Stable IDs are created in the form ENS\[species prefix\]\[feature type prefix\]\[a unique eleven digit number\]. For example a mouse gene might be ENSMUSG###########. This means that we can immediately tell from a stable ID what kind of feature they refer to and what species they are in. The prefixes are listed on our prefix page.

The two places I look first for regexes are identifiers.org and wikidata. If these are wrong, best to report that upstream. Here are their ensembl regexes:

The regex used by the identifiers.org registry (for any ensembl ID) is quite long:

```
^((ENS[FPTG]\\d{11}(\\.\\d+)?)|(FB\\w{2}\\d{7})|(Y[A-Z]{2}\\d{3}[a-zA-Z](\\-[A-Z])?)|([A-Z_a-z0-9]+(\\.)?(t)?(\\d+)?([a-z])?))$

```

Meanwhile Wikidata has this regex for an Ensembl gene ID:

```
(ENS(|MUS|RNO)G\d{11})|(Y[A-P][LR]\d{3}[CW](-[A-G])?)

```

---
Tags: python, ensembl
---


# Question
Title: ATAC seq peak annotation : Which peak to choose for downstream analysis
After doing peak call I'm annotating the peaks using chipseeker tool, which I want to take further downstream analysis. Now unlike RNA seq where i have single gene and their respective expression. But here the peaks for a single gene multiple peaks. As an example have filtered promoter peaks.

```
seqnames  start    end width strand     V4 annotation geneChr geneStart geneEnd geneLength geneStrand geneId
1     chr1 826797 828101  1305      *  Peak1   Promoter       1    826832  852225      25394          1 643837
2     chr1 869647 870324   678      *  Peak4   Promoter       1    868240  870201       1962          2 284593
3     chr1 876462 877140   679      *  Peak5   Promoter       1    868071  876903       8833          2 284593
4     chr1 877219 878105   887      *  Peak6   Promoter       1    874529  877234       2706          2 284593
5     chr1 923459 924318   860      * Peak14   Promoter       1    923928  939291      15364          1 148398
6     chr1 924379 926062  1684      * Peak15   Promoter       1    925150  935793      10644          1 148398
       transcriptId distanceToTSS         ENSEMBL    SYMBOL                                    GENENAME
1 ENST00000623808.3             0 ENSG00000228794 LINC01128 long intergenic non-protein coding RNA 1128
2 ENST00000432963.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
3 ENST00000446136.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
4 ENST00000635557.1             0 ENSG00000230368    FAM41C family with sequence similarity 41 member C
5 ENST00000420190.6             0 ENSG00000187634    SAMD11    sterile alpha motif domain containing 11
6 ENST00000437963.5             0 ENSG00000187634    SAMD11    sterile alpha motif domain containing 11

```

Now if we see first 6 rows we have two genes with three different peaks. When i load and see in the igv browser I do see the there are subtle differences in the accessibility signals.

How to address this issue ,in other words which peak to consider or take all the peak and consider their accessibility as one.

Any suggestion or help would be really appreciated.

# Answer
> 2 votes
I'll consider all peaks around promoter of particular gene as one peak. Moreover I'll interpret this as an open chromatin region = active gene. Why? Why there's no one wide peak which represents open chromatin in promoter of particular gene? Why is divided into several narrower peaks? There might be some transcription factors or other proteins binded to DNA, hence you see peak (open, free DNA) - break (valley, TF binded) - peak(again free DNA where Tn5 cut DNA and incorporate adapters) in IGV. In my study I compared ATAC-seq and H3k27Ac. In such case usually you may notice that peaks from chip and atac never literally overlap. Because ATAC represents open chromatin (DNA without any binded proteins) and chip histone modification (DNA is actually wound on histone proteins), hence we have something like this figure C, first we see high ATAC peak - open chromatin, then H3k27ac peak - lysine acetylated = gene activated. In other words ATAC-seq won't show you peaks around histones, binded polymerase, or transcription factors.

# Answer
> 2 votes
You can see from your table that the different peaks correspond to different promoters for different transcripts, for example, FAM41C has ENST00000432963.1, ENST00000446136.1, ENST00000635557.1. And the fact that they are all within the TSS (see distanceToTSS) suggest these peaks are actually falling onto actually TSS.

The underlying issue is are these TSS or do will we mis classify the peaks that are inside genes or other non-coding elements as TSS peaks. Nothing to do with TF, histones etc.

I cannot see the number of reads in those peaks, but if one of them is really strong, and the rest almost non-existent it would make sense.

How can we improve this, I think it depends on the aim of your analysis. Since you are looking at TSS, I guess it makes sense to be able to link it to something that is actually expressed. You can do:

1. To use an updated annotation (i see that some of the transcripts come from old hg37 annotation)
2. you have prior knowledge which transcripts might be expressed from RNA seq data, I would use that to reduce the number of transcripts

---
Tags: r, atac-seq
---


# Question
Title: Align the HIV-1 Gag protein with the Gag protein of Visna virus
**Align the HIV-1 Gag protein with the Gag protein of Visna virus**

**Background** Visna virus is a lentivirus causing encephalitis in sheep

The problem is align the follow two protein Gag of viruses HIV-1 (GenBank KU749412.1) and Gag Visna virus (GenBank L06906.1) that are:

```
"MGARASILRGEKLDTWERIRLRPGGKKHYMLKHLVWASRELDRFALNPGLLETLEGCKQIIKQLQPALQTGTEELRSLFNTVATLYCVHAGIPVKDTKEALDKIEEEQNKCQQKAKQAEAAAAGKISLNYPIVQNLQGQMVHQAISPRTLNAWVKVIEEKAFSPEVIPMFTALSEGATPQDLNTMLNTVGGHQAAMQMLKDTINEEAAEWDRLHPVHAGPVAPGQMREPRGSDIAGTTSTLAEQIAWITGNPPVPVGEIYKRWIILGLNKIVRMYSPXSILDIRQGPKEPFRDYVDRFFKTLRAEQATQDVKNWMTDTLLVQNANPDCKTILRALGPGATLEEMMTACQGVGGPGHKARVLAEAMSQTTNAIMMQRSNFKGPKRTIKCFNCGKEGHLARNCRAPRKKGCWKCGKEGHQMKDCNERQANFLGKIWPSHKGRPGNFLQSRPEPTAPPAESFRFEETTPAPKQEQKDKEPLTSLKSLFGSDPLSQ"
"MAKQGSKEKKGYPELKEVIKATCKIRVGPGKETLTEGNCLWALKTIDFIFEDLKTEPWTITKMYTVWDRLKGLTPEETSKREFASLQATLACIMCSQMGMKPETVQAAKGIISMKEGLHENKEAKGEKVEQLYPNLEKHREVYPIVNLQAGGRSWKAVESVVFQQLQTVAMQHGLVSEDFERQLAYYATTWTSKDILEVLAMMPGNRAQKELIQGKLNEEAERWVRQNPPGPNVLTVDQIMGVGQTNQQASQANMDQARQICLQWVITALRSVRHMSHRPGNPMLVKQKNTESYEDFIARLLEAIDAEPVTDPIKTYLKVTLSYTNASTDCQKQMDRTLGTRVQQATVEEKMQACRDVGSEGFKMQLLAQALRPQGKAGQKGVNQKCYNCGKPGHLARQCRQGIICHHCGKRGHMQKDCRQKKQQGNNRRGPRVVPSAPPML"

```

I propose that the better alignment method is the one that maximize the function

```
f(string1, string2)=(number of letter equal between strign1 and string2)/(common length of strings).

```

With that definition I found f(s1,s2)=0.25760 (and 0.2949 between POL protein) for the two strings above.

* What score might others find?
* Is that definition one of biological significative?
* Is there is a better one, what definition for align proteins?
* Are Visna and HIV viruses relatives?

# Answer
> 2 votes
Both are lentiviruses and `env` is a fast evolving gene. What you are looking at here is the location of indels and this will vary between alignment algorithms.

There will always be an element of error in attempting to align over this genetic distance. It depends what you end goal is, if it is to maximise homology than okay, however optical enhancement of alignments are common if your downstream analysis is phylogenetics. If it is protein structure homology style modelling, where you have the protein structure of HIV env experimentally determined mulitple times you want to use it to help estimate that of the cat lentivirus then optimising a scoring metric is okay.

What you have to understand is there is no 'right' method, in phylogenetics this doesn't matter as much because indels are excluded in the subsequent calculation, so it is less sensative to indel positioning errors and optical enhancement is used. You can also simply exclude regions that appear 'garbled'. For protein homology modeling this is a big problem, because the indels will make drastic changes to the structure prediction. So the only approach here in my opinion is to produce a separate structural prediction for every variation in indel alignment you have obtained. In the latter context getting the alignment right is everything and I suspect thats the output here, because its a classic surface antigen to understand, which will readily crystalise.

---
Tags: phylogenetics, sequence-alignment, phylogeny, proteins, virus
---


title: Biopython SeqIO check input file

text:
# Question
Title: Biopython SeqIO check input file
Hi I am trying to learnt python3 and Biopython,

I am trying to check imported fasta file before processing so far using:

```
sequences = SeqIO.parse("fake_fasta_file", "fasta")

if not any(sequa):
        print('\n  SOMETHING IS WRONG WITH YOUR INPUT FILE \n')
        exit()

```

works but as soon as I enter a '\>' in my fake\_fasta\_file , as an example

```
nihvivicòvo-l
,mnjobùJN
, Békibpç
, mnkn

nbigp
ù
mlkn

mlnù
>  ftfyfiuuogpigèi

```

any(sequences) returns True even if:

```
for seq_record in sequences:
       print(seq_record)

```

doesnt print anything

Any better way to check for real fasta file in Biopython ?

Or maybe it my fault using any() on the iterator instead of creating

a list of seq\_record and checking on that ?

# Answer
`any()` loops over the iterator until the first truthy item is found or the end is reached. The iterator is not reset to the start after this. Using a list will change that behaviour since a list is not exhausted like a iterator. So, remove your `any()` call and you will see a sequence record printed. Your `'fake_fasta_file.fa'` is actually recognised by Biopython as valid fasta:

```
>>> SeqIO.read('fake_fasta_file.fa', 'fasta')
SeqRecord(seq=Seq(''), id='ftfyfiuuogpigèi', name='ftfyfiuuogpigèi', description='  ftfyfiuuogpigèi', dbxrefs=[])

```

Biopython ignores lines before the first line beginning with `'>'` (as possible comment lines) and then starts collecting records, in this case an empty sequence record is still a record. This behaviour is deliberate from Biopython, although you could argue if it's best or not, it's unlikely to change for historical reasons.

What is *real* fasta? (philosophical pause) If you want to validate a particular definition of FASTA format then you can write your own parser depending on the exact specification and rules you want.

Edit, e.g.

```
from Bio import SeqIO

records = SeqIO.parse("myfile.fa", "fasta")
first = next(records, None)
if first is None:
    print("No fasta records")
elif first.seq == "":
    print("Empty sequence")
elif all(c not in "ATCGN" for c in first.seq.upper()):
    print("Non ATCGN characters")

```

etc.
> 1 votes

---
Tags: fasta, python, seqio
---