# Question
Title: How can I predict the folding of a synthetic protein?
I want to predict the folding of a chimeric antigen receptor sequence. These sequences are completely made out of small parts of different proteins. Do you know if AlphaFold can be used for this?

# Answer
I'd consult references about AlphaFold. For example, wikipedia has a "Limitations" page on the tool, though the page doesn't mention anything about synthetic proteins.

The hesitation that I would have is that the tool relies to some extent on coevolutionary models of amino acids to infer interactions at both training and prediction steps, so throwing a synthetic protein at it may confuse it. See e.g. slide 10 here.

In an interview, they acknowledge this rather directly:

> Now, if you have a protein which is completely new, which does not have an evolutionary history—a synthetic protein, which has just been designed— then our system would not be as accurate as it would be for existing proteins.

I am not totally clear on whether this is due to learning a latent representation of constraints on similar proteins, or whether the tool attempts to directly find homologs. The paper methods suggest to me that it is the former, but I'm not sure. They describe training on only chunks of a big MSA. If that's the case, I'd guess that they still can use information about the natural chunks of your synthetic protein, and aren't expecting to be able to do a full-length protein alignment.

So I think that you can expect the prediction to be *relatively* low quality. It is still probably much better than any alternative way of predicting structure!

As always, not a bad idea to post an issue on github asking this question.
> 1 votes

---
Tags: protein-structure
---


title: Remove low quality reads

text:
# Question
Title: Remove low quality reads
I want to remove reads from FASTQ file that contain homopolymers \> 10bp and remove reads with \<35 average quality score across the entire read.

To remove homopolymers \> 10bp, I tried this on a Linux machine, but it only removes the sequence line:

```
zcat file.fastq.gz |
  awk '!/A{10,}/&&!/C{10,}/&&!/G{10,}/&&!/T{10,}/ {print}' cleaned_file.fastq

```

# Answer
> 4 votes
Use an off the shelf tool for read preprocessing.

Here is one:

```
./fastq_qual_trimmer -i test.fq -m 35 -H 10

```

That one does exactly what you want but seems a little old/unmaintained, so here is another that does more or less the same thing:

```
fastq-mcf --qual-mean 35 --homopolymer-pct {X} adapters.fa reads.fq

```

where `{X}` is 10 / read length, adapters.fa is an adapter file (which I believe can be empty or filled with dummy sequences).

You could also use a library like biopython or dnaio to write a quick script to do this, but it hardly seems worth it.

---
Tags: fastq, awk
---


# Question
Title: expansions and contractions from OrthoFinder
The input file looks like this, and the complete file can be found here:

```
OG0000008 NbL07g11380.1 NbL19g07810.1 NbL19g09170.1 NbL19g19070.1 NbQ01g01670.1 NbQ01g03330.1 NbQ01g04070.1 NbQ01g04670.1 NbQ01g05120.1 NbQ01g05870.1 NbQ01g06940.1 NbQ01g07580.1 NbQ01g08860.1 NbQ01g10050.1 NbQ01g10360.1 NbQ01g14200.1 NbQ01g14790.1 NbQ01g16080.1 NbQ01g17760.1 NbQ01g19270.1 NbQ01g19310.1 NbQ01g19390.1 NbQ01g21260.1 NbQ01g21330.1 NbQ01g21740.1 NbQ01g21910.1 NbQ01g23100.1 NbQ01g24620.1 NbQ01g25340.1 NbQ01g26060.1 NbQ01g26320.1 NbQ02g00750.1 NbQ02g03100.1 NbQ02g03420.1 NbQ02g03610.1 NbQ02g03680.1 NbQ02g05120.1 NbQ02g07460.1 NbQ02g08170.1 NbQ02g08330.1 NbQ02g09220.1 NbQ02g09400.1 NbQ02g10620.1 NbQ02g11310.1 NbQ02g14330.1 NbQ02g14460.1 NbQ02g14520.1 NbQ02g15320.1 NbQ02g17090.1 NbQ02g17130.1 NbQ02g20290.1 NbQ02g23070.1 NbQ02g23420.1 NbQ02g24450.1 NbQ02g24480.1 NbQ02g26700.1 NbQ03g00830.1 NbQ03g01970.1 NbQ03g04460.1 NbQ03g06900.1 NbQ03g09530.1 NbQ03g10620.1 NbQ03g12760.1 NbQ03g13450.1 NbQ03g15540.1 NbQ03g15640.1 NbQ03g17180.1 NbQ03g20740.1 NbQ03g21510.1 NbQ03g24670.1 NbQ04g01350.1 NbQ04g01720.1 NbQ04g08420.1 NbQ04g09090.1 NbQ04g10450.1 NbQ04g11470.1 NbQ04g12120.1 NbQ04g14130.1 NbQ04g15440.1 NbQ04g15860.1 NbQ04g16450.1 NbQ04g16620.1 NbQ04g17760.1 NbQ04g19040.1 NbQ04g20020.1 NbQ05g03320.1 NbQ05g04660.1 NbQ05g05970.1 NbQ05g07500.1 NbQ05g08900.1 NbQ05g09760.1 NbQ05g10830.1 NbQ05g11150.1 NbQ05g11340.1 NbQ05g11510.1 NbQ05g11530.1 NbQ05g11780.1 NbQ05g16980.1 NbQ05g18190.1 NbQ05g21710.1 NbQ05g23400.1 NbQ06g01110.1 NbQ06g01430.1 NbQ06g04200.1 NbQ06g04440.1 NbQ06g05330.1 NbQ06g05770.1 NbQ06g05820.1 NbQ06g06700.1 NbQ06g08620.1 NbQ06g09190.1 NbQ06g10460.1 NbQ06g15220.1 NbQ06g15330.1 NbQ06g15700.1 NbQ06g16320.1 NbQ06g16590.1 NbQ06g17590.1 NbQ06g17670.1 NbQ06g20050.1 NbQ07g01030.1 NbQ07g02010.1 NbQ07g04350.1 NbQ07g04900.1 NbQ07g05610.1 NbQ07g06200.1 NbQ07g07110.1 NbQ07g07690.1 NbQ07g08640.1 NbQ07g10390.1 NbQ07g11920.1 NbQ07g14130.1 NbQ07g15590.1 NbQ07g15620.1 NbQ07g16910.1 NbQ07g17130.1 NbQ07g17950.1 NbQ08g00060.1 NbQ08g02240.1 NbQ08g02300.1 NbQ08g02310.1 NbQ08g03290.1 NbQ08g05330.1 NbQ08g09280.1 NbQ08g14890.1 NbQ08g15820.1 NbQ08g15950.1 NbQ08g19830.1 NbQ08g20150.1 NbQ08g22050.1 NbQ08g22620.1 NbQ09g02100.1 NbQ09g02620.1 NbQ09g03950.1 NbQ09g04200.1 NbQ09g06040.1 NbQ09g06640.1 NbQ09g08160.1 NbQ09g08330.1 NbQ09g09660.1 NbQ09g11220.1 NbQ09g13860.1 NbQ09g15180.1 NbQ09g15310.1 NbQ09g16530.1 NbQ09g17900.1 NbQ09g18100.1 NbQ09g18720.1 NbQ09g19280.1 NbQ09g21840.1 NbQ10g00480.1 NbQ10g01350.1 NbQ10g02870.1 NbQ10g03640.1 NbQ10g03730.1 NbQ10g08070.1 NbQ10g09510.1 NbQ10g11010.1 NbQ10g11760.1 NbQ10g12050.1 NbQ10g12060.1 NbQ10g12910.1 NbQ10g19200.1 NbQ10g19930.1 NbQ10g20390.1 NbQ10g20730.1 NbQ10g21080.1 NbQ10g21140.1 NbQ10g24010.1 NbQ11g00310.1 NbQ11g01210.1 NbQ11g01370.1 NbQ11g04610.1 NbQ11g04800.1 NbQ11g06060.1 NbQ11g07820.1 NbQ11g08390.1 NbQ11g09100.1 NbQ11g09350.1 NbQ11g13660.1 NbQ11g13930.1 NbQ11g16260.1 NbQ11g17360.1 NbQ11g18430.1 NbQ11g21080.1 NbQ11g23280.1 NbQ11g23990.1 NbQ11g25050.1 NbQ12g03770.1 NbQ12g04850.1 NbQ12g07340.1 NbQ12g09080.1 NbQ12g10820.1 NbQ12g12070.1 NbQ12g14750.1 NbQ12g15000.1 NbQ12g15230.1 NbQ12g20380.1 NbQ12g21080.1 NbQ12g21830.1 NbQ12g23960.1 NbQ13g01300.1 NbQ13g02350.1 NbQ13g03860.1 NbQ13g04410.1 NbQ13g08800.1 NbQ13g09850.1 NbQ13g10370.1 NbQ13g11700.1 NbQ13g12420.1 NbQ13g15780.1 NbQ13g16040.1 NbQ13g23160.1 NbQ13g24120.1 NbQ13g24540.1 NbQ13g25080.1 NbQ13g25490.1 NbQ13g28240.1 NbQ13g29770.1 NbQ14g01070.1 NbQ14g03950.1 NbQ14g05360.1 NbQ14g05410.1 NbQ14g06880.1 NbQ14g07270.1 NbQ14g07500.1 NbQ14g10290.1 NbQ14g10770.1 NbQ14g14320.1 NbQ14g17890.1 NbQ14g18710.1 NbQ14g20960.1 NbQ14g22890.1 NbQ15g00150.1 NbQ15g02300.1 NbQ15g02330.1 NbQ15g02350.1 NbQ15g03230.1 NbQ15g06190.1 NbQ15g07120.1 NbQ15g07750.1 NbQ15g09000.1 NbQ15g09050.1 NbQ15g11920.1 NbQ15g12650.1 NbQ15g12840.1 NbQ15g15670.1 NbQ15g15930.1 NbQ15g18670.1 NbQ15g19070.1 NbQ15g20620.1 NbQ15g22880.1 NbQ15g23000.1 NbQ15g26060.1 NbQ16g00880.1 NbQ16g04360.1 NbQ16g06490.1 NbQ16g09100.1 NbQ16g11020.1 NbQ16g11560.1 NbQ16g13810.1 NbQ16g13820.1 NbQ16g17040.1 NbQ16g17130.1 NbQ16g17340.1 NbQ16g18390.1 NbQ16g18430.1 NbQ16g23100.1 NbQ16g23570.1 NbQ16g24270.1 NbQ16g25200.1 NbQ16g25830.1 NbQ16g25880.1 NbQ16g25990.1 NbQ16g26610.1 NbQ16g26660.1 NbQ16g28010.1 NbQ16g28180.1 NbQ17g01150.1 NbQ17g01180.1 NbQ17g01570.1 NbQ17g01950.1 NbQ17g05460.1 NbQ17g05540.1 NbQ17g05980.1 NbQ17g07990.1 NbQ17g08300.1 NbQ17g09330.1 NbQ17g09400.1 NbQ17g10090.1 NbQ17g11220.1 NbQ17g13030.1 NbQ17g15460.1 NbQ17g16690.1 NbQ17g20980.1 NbQ17g22370.1 NbQ17g25040.1 NbQ17g28730.1 NbQ18g02140.1 NbQ18g02740.1 NbQ18g05440.1 NbQ18g06120.1 NbQ18g07470.1 NbQ18g12320.1 NbQ18g12530.1 NbQ18g12850.1 NbQ18g13840.1 NbQ18g14420.1 NbQ18g14930.1 NbQ18g15730.1 NbQ18g17750.1 NbQ18g17850.1 NbQ18g21060.1 NbQ19g01040.1 NbQ19g05480.1 NbQ19g06450.1 NbQ19g06510.1 NbQ19g08330.1 NbQ19g11840.1 NbQ19g11880.1 NbQ19g13750.1 NbQ19g14190.1 NbQ19g14210.1 NbQ19g14920.1 NbQ19g18540.1 NbQ19g19870.1 NbQ19g21020.1 NbQ19g21220.1 NbQ19g22080.1 NbQ19g22800.1 NbQ19g24690.1 NbQ19g24730.1 rna19561
OG0000001 Capann_59V1aChr01g048170.1 NbL01g00020.1 NbL01g00940.1 NbL01g02330.1 NbL01g03550.1 NbL01g03650.1 NbL01g04410.1 NbL01g04920.1 NbL01g06850.1 NbL01g16120.1 NbL01g19150.1 NbL01g20140.1 NbL01g20930.1 NbL01g22230.1 NbL01g24190.1 NbL01g24280.1 NbL01g24300.1 NbL02g00570.1 NbL02g00900.1 NbL02g01270.1 NbL02g02110.1 NbL02g02210.1 NbL02g02470.1 NbL02g03180.1 NbL02g04740.1 NbL02g04750.1 NbL02g06120.1 NbL02g06860.1 NbL02g07280.1 NbL02g07680.1 NbL02g07740.1 NbL02g09780.1 NbL02g11320.1 NbL02g12670.1 NbL02g13080.1 NbL02g14050.1 NbL02g14190.1 NbL02g15010.1 NbL02g15890.1 NbL02g16190.1 NbL02g16730.1 NbL02g17070.1 NbL02g17360.1 NbL02g18820.1 NbL02g19340.1 NbL02g20100.1 NbL02g23950.1 NbL02g24800.1 NbL03g01610.1 NbL03g01680.1 NbL03g01890.1 NbL03g02230.1 NbL03g02600.1 NbL03g03410.1 NbL03g04990.1 NbL03g05400.1 NbL03g08030.1 NbL03g08250.1 NbL03g08690.1 NbL03g10230.1 NbL03g11060.1 NbL03g13030.1 NbL03g14960.1 NbL03g15110.1 NbL03g16690.1 NbL03g16900.1 NbL03g18260.1 NbL03g18950.1 NbL03g21180.1 NbL03g21210.1 NbL03g21530.1 NbL03g22960.1 NbL03g24430.1 NbL04g01140.1 NbL04g01490.1 NbL04g02030.1 NbL04g02560.1 NbL04g03700.1 NbL04g04160.1 NbL04g05240.1 NbL04g05420.1 NbL04g05850.1 NbL04g12420.1 NbL04g12640.1 NbL04g13650.1 NbL04g13780.1 NbL04g14310.1 NbL04g16260.1 NbL04g17750.1 NbL04g18380.1 NbL04g18870.1 NbL04g19030.1 NbL04g19630.1 NbL05g00320.1 NbL05g03060.1 NbL05g03300.1 NbL05g04060.1 NbL05g07620.1 NbL05g08630.1 NbL05g09580.1 NbL05g10060.1 NbL05g11400.1 NbL05g12280.1 NbL05g13170.1 NbL05g16020.1 NbL05g17530.1 NbL05g17730.1 NbL05g18340.1 NbL05g18590.1 NbL05g18600.1 NbL05g19640.1 NbL05g20640.1 NbL05g21000.1 NbL05g21640.1 NbL05g22610.1 NbL06g00640.1 NbL06g00660.1 NbL06g02210.1 NbL06g03150.1 NbL06g03680.1 NbL06g04910.1 NbL06g07950.1 NbL06g09970.1 NbL06g11480.1 NbL06g12220.1 NbL06g12400.1 NbL06g12460.1 NbL06g12850.1 NbL06g13120.1 NbL06g14450.1 NbL06g14780.1 NbL06g16990.1 NbL06g17200.1 NbL06g17760.1 NbL06g20380.1 NbL07g02100.1 NbL07g02540.1 NbL07g02970.1 NbL07g03110.1 NbL07g04840.1 NbL07g05350.1 NbL07g06580.1 NbL07g07530.1 NbL07g08450.1 NbL07g09380.1 NbL07g09870.1 NbL07g10730.1 NbL07g10850.1 NbL07g11080.1 NbL07g12450.1 NbL07g12710.1 NbL07g13110.1 NbL07g13920.1 NbL07g14240.1 NbL07g15520.1 NbL07g16220.1 NbL07g17480.1 NbL08g01820.1 NbL08g02750.1 NbL08g02930.1 NbL08g03510.1 NbL08g03620.1 NbL08g03850.1 NbL08g03970.1 NbL08g04040.1 NbL08g06150.1 NbL08g06410.1 NbL08g06680.1 NbL08g06730.1 NbL08g07620.1 NbL08g08450.1 NbL08g08640.1 NbL08g09910.1 NbL08g10160.1 NbL08g11760.1 NbL08g12570.1 NbL08g13630.1 NbL08g13890.1 NbL08g15050.1 NbL08g15340.1 NbL08g18010.1 NbL08g18420.1 NbL08g19080.1 NbL08g19190.1 NbL09g00900.1 NbL09g02160.1 NbL09g02330.1 NbL09g02470.1 NbL09g04150.1 NbL09g05210.1 NbL09g07010.1 NbL09g09070.1 NbL09g10290.1 NbL09g10500.1 NbL09g11220.1 NbL09g13490.1 NbL09g15290.1 NbL09g15830.1 NbL09g17240.1 NbL09g19250.1 NbL09g19460.1 NbL09g20190.1 NbL09g21040.1 NbL09g21520.1 NbL09g23460.1 NbL10g00080.1 NbL10g03710.1 NbL10g04330.1 NbL10g04560.1 NbL10g05200.1 NbL10g06320.1 NbL10g07510.1 NbL10g07960.1 NbL10g08670.1 NbL10g08970.1 NbL10g11120.1 NbL10g11340.1 NbL10g11820.1 NbL10g13720.1 NbL10g14560.1 NbL10g14770.1 NbL10g16430.1 NbL10g18140.1 NbL10g18380.1 NbL10g19280.1 NbL10g19690.1 NbL10g21210.1 NbL10g22680.1 NbL10g23160.1 NbL10g23560.1 NbL10g24210.1 NbL11g00680.1 NbL11g00970.1 NbL11g01230.1 NbL11g01270.1 NbL11g01520.1 NbL11g01530.1 NbL11g02920.1 NbL11g03540.1 NbL11g03990.1 NbL11g05630.1 NbL11g08950.1 NbL11g08980.1 NbL11g09510.1 NbL11g10840.1 NbL11g11030.1 NbL11g11230.1 NbL11g12430.1 NbL11g13300.1 NbL11g15430.1 NbL11g16390.1 NbL11g16410.1 NbL11g17320.1 NbL11g18090.1 NbL11g21310.1 NbL11g21470.1 NbL11g21780.1 NbL11g21820.1 NbL11g22270.1 NbL11g22310.1 NbL11g23180.1 NbL11g24100.1 NbL12g00130.1 NbL12g01810.1 NbL12g02230.1 NbL12g02720.1 NbL12g02760.1 NbL12g04120.1 NbL12g04550.1 NbL12g06630.1 NbL12g07830.1 NbL12g09170.1 NbL12g10580.1 NbL12g12090.1 NbL12g12490.1 NbL12g12630.1 NbL12g12800.1 NbL12g13320.1 NbL12g13460.1 NbL12g14430.1 NbL12g14970.1 NbL12g15490.1 NbL12g17460.1 NbL12g18190.1 NbL12g18590.1 NbL12g19900.1 NbL12g20690.1 NbL12g22040.1 NbL12g22560.1 NbL13g00350.1 NbL13g01440.1 NbL13g02400.1 NbL13g03210.1 NbL13g03360.1 NbL13g04070.1 NbL13g05250.1 NbL13g08460.1 NbL13g09010.1 NbL13g09140.1 NbL13g10290.1 NbL13g11570.1 NbL13g13370.1 NbL13g14910.1 NbL13g18680.1 NbL13g19510.1 NbL13g23520.1 NbL13g24010.1 NbL13g24190.1 NbL13g24460.1 NbL13g26310.1 NbL13g26640.1 NbL13g26860.1 NbL13g27260.1 NbL13g27960.1 NbL14g02460.1 NbL14g02750.1 NbL14g08750.1 NbL14g08910.1 NbL14g09120.1 NbL14g09540.1 NbL14g09920.1 NbL14g11070.1 NbL14g11150.1 NbL14g12570.1 NbL14g14530.1 NbL14g14860.1 NbL14g15240.1 NbL14g15460.1 NbL14g16620.1 NbL14g16910.1 NbL14g17940.1 NbL14g21150.1 NbL14g21750.1 NbL14g21910.1 NbL15g00790.1 NbL15g01170.1 NbL15g02310.1 NbL15g04220.1 NbL15g05970.1 NbL15g06340.1 NbL15g06440.1 NbL15g07020.1 NbL15g07370.1 NbL15g07470.1 NbL15g09010.1 NbL15g13210.1 NbL15g14550.1 NbL15g14600.1 NbL15g17290.1 NbL15g18170.1 NbL15g19710.1 NbL15g21840.1 NbL15g21930.1 NbL15g23410.1 NbL15g23420.1 NbL15g23430.1 NbL15g25130.1 NbL15g25200.1 NbL16g01760.1 NbL16g02140.1 NbL16g04460.1 NbL16g05010.1 NbL16g05020.1 NbL16g06780.1 NbL16g07540.1 NbL16g07980.1 NbL16g09760.1 NbL16g10610.1 NbL16g12320.1 NbL16g13510.1 NbL16g14420.1 NbL16g15690.1 NbL16g17420.1 NbL16g17790.1 NbL16g17880.1 NbL16g18730.1 NbL16g18940.1 NbL16g19440.1 NbL16g20980.1 NbL16g23180.1 NbL16g23610.1 NbL16g23660.1 NbL16g23910.1 NbL16g24550.1 NbL16g24640.1 NbL16g25300.1 NbL16g25630.1 NbL16g26710.1 NbL17g01590.1 NbL17g02070.1 NbL17g02120.1 NbL17g02920.1 NbL17g03040.1 NbL17g03540.1 NbL17g03700.1 NbL17g03800.1 NbL17g05400.1 NbL17g07510.1 NbL17g08450.1 NbL17g08930.1 NbL17g10090.1 NbL17g14370.1 NbL17g14600.1 NbL17g15390.1 NbL17g15900.1 NbL17g16000.1 NbL17g16910.1 NbL17g17480.1 NbL17g18240.1 NbL17g20020.1 NbL17g20830.1 NbL17g21220.1 NbL17g21690.1 NbL17g25960.1 NbL18g00030.1 NbL18g00150.1 NbL18g00310.1 NbL18g00670.1 NbL18g00700.1 NbL18g01630.1 NbL18g02650.1 NbL18g04460.1 NbL18g05210.1 NbL18g05690.1 NbL18g07270.1 NbL18g07440.1 NbL18g07500.1 NbL18g09090.1 NbL18g09810.1 NbL18g09880.1 NbL18g10500.1 NbL18g10990.1 NbL18g12070.1 NbL18g13060.1 NbL18g17480.1 NbL19g00230.1 NbL19g04470.1 NbL19g04700.1 NbL19g07770.1 NbL19g07870.1 NbL19g09260.1 NbL19g09870.1 NbL19g13480.1 NbL19g14360.1 NbL19g14400.1 NbL19g14720.1 NbL19g17700.1 NbL19g18860.1 NbL19g22750.1 NbQ13g00370.1 NbQ14g10310.1 NbQ17g06050.1

```

This script appears to have some problems with the file:

```
import pandas as pd

orthofinder_output = "../OrthoFinder-res/Results_Jan23/Orthogroups/Orthogroups-fixed.txt"
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None)

# Extract gene family information
expansions = {}
contractions = {}
for i, row in orthogroups.iterrows():
    if len(row) > 2:
        expansions[row[0]] = row[2:]
    else:
        contractions[row[0]] = row[1]

# Print expansions and contractions
print("Expansions:", expansions)
print("Contractions:", contractions)

```

I got the following error:

```
geneFamilyExpansionsContractions.py:
 20: DtypeWarning:
   Columns (2,3,4,5,6,...,964,965,966,967,968) have mixed types.
   Specify dtype option on import or set low_memory=False.
  orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None)
...
Expansions:
...
Name: 63241, Length: 967, dtype: object, 'OG0063242':
2      NaN
3      NaN
4      NaN
5      NaN
6      NaN
      ...
964    NaN
965    NaN
966    NaN
967    NaN
968    NaN
Name: 63242, Length: 967, dtype: object}
Contractions: {}

```

How can it be fixed?

# Answer
> 2 votes
The answer is actually right there in the error message (warning message, actually):

> DtypeWarning: Columns (2,... 968) have mixed types. Specify dtype option on import or set low\_memory=False.

Setting `low_memory-False` makes the error go away:

```
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None, low_memory=False)

```

According to this random blog post I found when I searched for "pandas DtypeWarning: have mixed types", this is because:

> Pandas is really nice because instead of stopping altogether, it guesses which dtype a column has. However, by default, Pandas has the low\_memory=True argument. This means that the CSV file gets split up in multiple chunks and the guess is done for every chunk, resulting in a column with multiple dtypes. Pandas is so kind to let you know it was confused and something might have happened.

The same blog also suggests that this isn't really the right approach, and that makes sense, since this is forcing Pandas to read the entire file and that might be a problem for huge files or low memory machines. So, reading on, the blog suggests the other option provided by the error, setting the `dtype` manually.

Now, I am afraid I have never used Pandas and don't really know much Python either, but a quick glance at the Pandas documentation suggests that the `object` data type is the most versatile and should work for everything, including strings. Since your file is nothing but a collection of strings, as far as I can tell, this seems like the best approach.

So I changed your script to this:

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

import pandas as pd
import sys

orthofinder_output = sys.argv[1]
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None, low_memory=False)
# Extract gene family information
expansions = {}
contractions = {}
for i, row in orthogroups.iterrows():
    if len(row) > 2:
        expansions[row[0]] = row[2:]
    else:
        contractions[row[0]] = row[1]

# Print expansions and contractions
print("Expansions:", expansions)
print("Contractions:", contractions)

```

And it seems to work just fine with no warnings shown.

# Answer
> 1 votes
The error right at the end is

> Specify `dtype` option on import or set `low_memory=False`.

**Edit** Oh wait @terdon has already spotted this! :-D ... great minds and all that :)

In addition, the `pandas` print out it states ...

```
dtype object

```

This isn't good. `Name: 63242,` is a bit unusual BTW worth checking.

Looking at the values ...

```
NaN

```

... that is really not good.

What I think is happening is it's been assigned to a numerical datatype, but it's categorical data. Hence `NaN` because like the alphabet can be none-base 10 numbers, in data compression used a lot, but it's got to be consistent pattern (I genuinely don't know if pandas can spot/work in none-base 10 number systems). `pandas` is really `C` (probably `C++`) - its a static type backend - so those datatypes are the key to this.

It makes sense because your data is categorical data and `NaN` means its managed to confuse a large amount of it for e.g. `int64`.

Looking the documentation the alternative is

```
orthogroups = pd.read_csv(orthofinder_output, sep=' ', header=None, dtype='category')

```

The documentation is

> low\_memorybool, default True Internally process the file in chunks, resulting in lower memory use while parsing, but possibly mixed type inference. To ensure no mixed types either set False, or specify the type with the dtype parameter. Note that the entire file is read into a single DataFrame regardless, use the chunksize or iterator parameter to return the data in chunks. (Only valid with C parser).

What you want is to issue a `df.dtypes` command once it is loaded and the output will define the datatype per column,

```
'A'    int64
'B'    int32
'C'   category

```

If you have mixed dtypes it *may* (I stress *may*) be able to correct it via,

```
df.astype('int64') # this is not your data
df.astype('category') # this is your data

```

HOWEVER, I don't think this will work, because `NaN` indicates its trashed the data. Anyway the data is definitely categorical BTW.

Anyway, good question and interesting stuff.

---

As a bit of gossip, you can also switch to categorical ordered data I think this will work for dataframes, this example is `Series`. There's other ways to order the data BTW.

```
from pandas.api.types import CategoricalDtype
my_dtype = CategoricalDtype(\
    categories=[2, 1], ordered=True)
df.dtype(my_dtype)

```

The output is `1` and `2` now switch places because 2 \< 1. Categorial datatypes look to be the right datatype for your inputs.

---
Tags: python, proteins, pandas, orthofinder
---


# Question
Title: Exon number between GTF file and ENSEMBL site does not match
Based on the ENSEMBL GTF file, I counted the number of exons per gene as follows:

```
awk '$3 == "exon" {print}' Homo_sapiens.GRCh38.108.gtf | gffread -F --keep-exon-attrs --table "gene_id","transcript_id",@numexons | sort -k1,1 | datamash -sg1 max 3 > gene_lists/exons_per_gene.txt

```

An extract of the output is:

```
ENSG00000034677 11
ENSG00000034693 12
ENSG00000034713 4
ENSG00000034971 4

```

However, when I check the ENSG00000034713 ENSEMBL site, I count 5 exons. Why is this different?

# Answer
*Answer from @terdon, @Ben\_Ensembl, converted from comments:*

Genes don't have exons, *transcripts* have exons. There are five transcripts for ENSG00000034713, one with one exon, two with three exons, and two with four exons. No transcripts with five exons though.

The count you retrieved from the GTF will be the number of unique exons for all transcripts of each gene (some unique to a single transcript, some common between two or more transcripts).

*\[please improve this answer if you can by editing it\]*
> 2 votes

---
Tags: ensembl, gtf
---


# Question
Title: Can 5' end 10x be used instead of 3' end
I have PBMCs of some patients who have recovered from cancer, and I want to look at their immune landscapes.

I got confused when looking at available kits: can I use 5' end kits from 10x instead of 3' end based?

I tried to search for information on this, but I am not clear if 5' can substitute for 3' or if they are just too different.

# Answer
> 3 votes
I'm not entirely clear on your question, but if it is about the 10X 3' vs. 5' solutions then either can be used to sequence PBMCs. The principle difference is that 3' captures the 3' end of the transcript and 5' captures the 5' end. The differences are more important for multi-modal single cell data, e.g. only the 5' can be combined with VDJ sequencing and for CITE-seq, different antibody panels are compatible with each tech. There are also some more subtle differences and technical nuances.

# Answer
> 3 votes
It depends on the task. If you want to recover TCR and BCR repertoires along with gene expression - 5' is a must. Also, for an overview repertoire analysis you dont even have to perform target enrichment and use 5' RNA-seq daata to directly recover the repertoire.

Its easy to do with MiXCR software in one line command:

```
mixcr analyze 10x-5gex-full-length \
 --species hsa \
 sample_R1.fastq.gz \
 sample_R2.fastq.gz \
 sample_result

```

Though MiXCR currently only supports repertoire analysis it will also support gene expression soon.

More on that here: https://docs.milaboratories.com/mixcr/reference/overview-built-in-presets/#10xgenomics

Disclaimer: I work for the company that develops MiXCR.

---
Tags: scrnaseq, single-cell, 10x
---


title: Trim reads 1kb upstream of sequence

text:
# Question
Title: Trim reads 1kb upstream of sequence
I need a quick way to trim multiple reads in a FASTA file. I need to trim everything that is 1kbp upstream of this sequence `AAGAGATGTTCAATCGTTTAAACAAATTCCAAGCTGCTTTAGCTTTGGCCCTTTACTCTCA`.

I figure a quick python script might be the way to go but I'm not sure if there's a tool that already does that (i.e. Trimmomatic or something like it) that would be able to do that for me.

Thank you!

# Answer
> 2 votes
Here's one way using . Note that the `Seq` object has a number of methods that act just like those of a Python string. One of these is the `find()` method, which returns the index of the first occurrence of the sequence if found. Replace with `rfind()` if you want the last occurrence if that makes more sense. This solution uses a generator to avoid storing the entire list of trimmed sequences in memory:

```
import argparse
import contextlib
import gzip
import sys

from Bio import SeqIO

```

```
def trim_fasta(fileobj, sequence, num_bases):
    """
    Trim num_bases upstream of the sequence provided
    """
    for record in SeqIO.parse(fileobj, 'fasta'):
        index = record.seq.find(sequence)
        if index != -1 and index > num_bases:
            record = record[index-num_bases:]
        yield record

```

```
def open_file(filename, mode='r'):
    if filename.endswith('.gz'):
        return gzip.open(filename, mode)
    else:
        return open(filename, mode)

```

```
def get_argument_parser():

    parser = argparse.ArgumentParser(add_help=False)

    group = parser.add_argument_group('input options')
    group.add_argument(
        "fasta",
        type=str,
        metavar="FILE",
        help="The input FASTA file",
    )
    group.add_argument(
        "sequence",
        type=str,
        metavar="STR",
        help="The sequence to find and trim",
    )

    group = parser.add_argument_group('trimming options')
    group.add_argument(
        "-n",
        "--num_bases",
        type=int,
        metavar="INT",
        default=1000,
        help="Trim INT bases upstream of input sequence",
    )

    group = parser.add_argument_group('output options')
    group.add_argument(
        "-o",
        "--output",
        type=str,
        metavar="FILE",
        default='-',
        help="Write the filtered output to FILE (default: stdout)",
    )
    group.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )

    group = parser.add_argument_group('additional options')
    group.add_argument(
        "-h",
        "--help",
        action="help",
        help="Show this help message and exit",
    )

    return parser

```

```
def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with contextlib.ExitStack() as stack:
        if args.output != '-':
            out = stack.enter_context(open_file(args.output, 'wt' if args.force else 'xt'))
        else:
            out = sys.stdout

        input_fasta = stack.enter_context(open_file(args.fasta))
        trimmed_reads = trim_fasta(input_fasta, args.sequence, args.num_bases)

        SeqIO.write(trimmed_reads, out, "fasta-2line")

if __name__ == '__main__':
    main()

```

---

---

```
usage: trim_reads.py [-n INT] [-o FILE] [-f] [-h] FILE STR

input options:
  FILE                  The input FASTA file
  STR                   The sequence to find and trim

trimming options:
  -n INT, --num_bases INT
                        Trim INT bases upstream of input sequence

output options:
  -o FILE, --output FILE
                        Write the filtered output to FILE (default: stdout)
  -f, --force           Overwrite the output file if it exists

additional options:
  -h, --help            Show this help message and exit

```

---
Tags: python, trimming, multi-fasta
---


# Question
Title: Find open reading frames in a DNA sequence
I have a fast file containing DNA sequences. I've already cleaned the file and has only DNA sequences with no gaps or errors or white space characters. I want to find open reading frames in the fasta file that are longer than 30 bp.

Here's a little description: An open reading frame is defined by a sequence that

1. start with `ATG`,
2. ends with `TAG`, `TAA`, or `TGA`,
3. is a multiple of three.

For each open reading frame, I want to print the start position in the sequence, the length of the open reading frame, and the open reading frame sequence to the screen. Then I want to print the total number of open reading frames found. I would appreciate some help.

I can't use libraries.

Also, here's the DNA sequence I am working with:

```
GCATGCAATACAGTGACATATATATATACCCTAACACTACCCTAACCCTACCCTATTTCA
ACCCTTCCAACCTGTCTCTCAACTTACCCTCACATTACCCTACCTCTCCACTTGTTACCC
TGTCCCATTCAACCATACCACTCCCAACCACCATCCATCCCTCTACTTACTACCACCAAT
CAACCGTCCACCATAACCGTTACCCTCCAATTAGCCATATTCAACTTCACTACCACTTAC
CCTGCCATTACTCTACCATCCACCATCTGCTACTCACCATACTGTTGTTCTACCCTCCAT
ATTGAAACGTTAACAAATGATCGTAAATAATACACATATACTTACCCTACCACTTTATAC
CACCACACATCACATGCCATACTCACCTTCACTTGTATACTGATATGCCATACGCACACG
GATGCTACAGTATATACCACTCTCAAACTTACCCTACTCTCACATTCTACTCCACTCCAT
GACCCATCTCTCACTAAATCAGTACTAAATGCACCCACATCATTATGCACGGCACTTGCC
TCAGCGGTCTATACCCTGAGCCATTTACCCATAACTCCCACGATTATCCACATTTTAATA
TCTATATCTCATTCGGCGGGCCCAAATATTGTATAACTGCTCTTAATACATACGTTATAC
CACTTTTGCACCATATACTAACCACTCAAATTATATACACTTATGCCAATATAACCAAAA
AATTACCACTAAAATCACCTAAACATAAAAATTATTTATCTTTCAACTTTACGAAATAAA
CACACTCAATTGCGTATCTATACCACCATGACGTCATTAACGTAAAAGTTCCTTAATATT
ACCATTTGCTTGAACGGATACCATTTCAGAATATTTCTAACTTTCACAGACCATACATTA
GAATAATATGCCACCTCACTGTCGTAACACTCTATATTCACCGAGAAACAATACGGTAGT
GGCTCAAACTCATGCCGGTGCTATGATACAATTGTATCCTATTTCCATTCTCATATGCTA
TCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCAT
ACGAAGAGACTATATCTAAAGAAGACGATACAGTGATATGTACGTTGTTTTTGTAGAATT
ATAATGAAACGTCAAATAACCCTACTATATTATAACTTATCAGCGCCGTATACTAAAACG
GACGTTACGATATTGTCTCACTTCATCTTACCACCCTCTATCTTATTGTTGATAAAACAC
TAACCCCTCAGCTTTATTTCTAGTTACAGTTACAACAAACTATCCCAAACCATAAATCTT
AATATTTTAGGTGTCAAAAAATGAGGATCTCCAAATGAGAGTTTGGTACCATGACTTGTA
ACTCCACTACCCTGATCTGCAATCTTGTTCTTAGAAGTGACGCATACTCTATATGGCCCG
ACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGT
GCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTT
ATTACCGTCCACGCTTTTTCTACAATAGTGTAAAAGTTTCTTTCTTATGTTCATCGTATT
CATAAAATGCTTCACGAACACCGTCATTGATCAAATAGGTTTATAATATTAATATACATT
TATATAATCGGCGGTATTTATATCATCAAAAAAAGTAGTTTTTTATTTTATTTTTTCATT
ACTTTTCACTGTCTATGGATTTTCATTCGTAAAGGCATCACTCCCTAGTTTGCGATAGTG
TAGATACCGTCCTTGGATAGAGCACTGGAGATGGCTGGCTTTAATCTGCTGGAGTACCAT
GGAACACCGGTGATCATTCTGGTCACTTGGTCTGGGGCAATACCAGTCAACATGGTGGTG
AAGTCACCGTAGTTGAAAACGGCTTCAGCAACTTCGACTGGGTAGGTTTCAGTTGGGTGG
GCGGCTTGGAACATGTAGTATTGGGCCAAGTGAGCTCTGATATCAGAGACGTAGACACCC
AATTCCACCAAGTTGACTCTTTCGTCAGATTGAGCTAGAGTGGTGGTTGCAGAAAGCAGT
AGCAGCGATGGCAGCGACACCAGCAGCGATTGAAGTTAATTTGACCATTGTATTTGTTTT
GTTTGTTAGTGCTGATGTAATCTTAACAAGAAATAGTGAAATGAAAGCGCATACCTCAAA
GGCATATAGTTGAAGCAGCTCTATTTATACCCGTTCCTCCATCTGTCATCACTACTTAAA
CGATTCGTTAACAGACGCTCATTTAGCACCTCACATATTCTCCATATCTCATCTTTCACA
CAATCTCATTATCTCTATGGAGATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAG
GTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGG
GAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAA
GTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGC
TTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTT
GTTACTCATTCCTGAAGGTTGAATAAAAAATAGCAACCATGTGCAGGGAGTCGTATATTG
TTAGATTCTAGAGACTTGTACGCATCATCAAAGCTGTAAATAGAATAAACATACGCAAGG
CGTCAAAAGTGCATAGTTAAGAAAATTCCTGACATGTGAAAATATGTGTTTATGAAATGT
GTCAAGGCCCGTCTATAGCGTAGTTAACCCCTCTGCAGGAGTAAGTGACTTTTTTTACGC
TCAAAAGGCAACGAGGGCACATACTTAAAAGTCATTTTCAAACACATCTGCAGTTTGCAA
CGACAGATAACAATATTATGATAGGATGGTATGATGTTATTAGCTGCCACATATTTTTAA
TAATAGCATTAGTCACGTCTCTTCAATTGTTGGGATGAAACTCTAAAATATCATTCCTTT
AGTAGTATTCCAGTTACCAGTATATTATCACATGCCGAAAAAGAAGATGACATAAAGATC
GACAAACAGTCTTCAAATATAATGGAAGCTGGAATGCAAGGATTGATAATGTAACAGGAT
ACTGAATGACAAAGTATAAATGAAAAAAAAAAAAAAAAAAAAGTAGTAATACTATTATGT
GGAAATACCGATTCCATTTTGAGGATTCCAATTGTTGGAATAAAAATCAACTATCATCTA
CTAACTAGTATTTACGTTACTAGTATATTATCATATACGGTGTTAGAAGATGACGCAAAT
GATGAGAAATAGTCATCTAAATTAGTGGAAGCTGAAACGCAAGGATTGATAATGTAATAG
GATCAATGAATATTAACATATAAAACGATGATAATAATATTTATAGAATTGTGTAGAATT
GCAGATTCCCTTTTATGGATTCCTAAATCCTGGAGGAGAACTTCTAGTATATCTACATAC
CTAATATTATAGCCTTAATCACAATGGAATCCCAACAATTACATCAAAATCCACATTCTC
TACACCAATACCATCGACGAGAGCTTCTAGTAAATTGTATACATAACAGTATAACCCTTA
CCAACAATGGAATCTCAAAGATTATTAAATTATTCACAGACTCTGAGGATTCGGGTAAAA
TAGGGTATTTAACTGGTTACCGGAAAGGTTTAGAAAATTCGTGGAGGGTTGGCCGAGTGG
TCTAAGGCGGCAGACTTAAGATCTGTTGGACGGTTGTCCGCGCGAGTTCGAACCTCGCAT
TCATTTTCCGGATCCACTTTATAGTGGGTCACTTCATCTCTTGTAGGGTTCACTGAAAGA
TAATCAATCAAATTATCATATGACCAATGGGTCACTTTTTTGATGAGACCCATGTTTATT
CTTCTATACGTTCGTGATACTGCACTTACGACTACCAGTAACATCAGAATAACCATAGAG
CTTATGCTATCTTGAAATATATCATCGCCGTAAAGTATGATGTGTAGCGCATTAGTCTTC
TGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGG
CTCTACTCTATATAAATTCGACATCCATTGATTTTCTCAATAAATGAGTTACCCAAGTAA
GCTTTCCATTGATACAAGTGATCTACATTCTTGCGACGCCAAATTATAAAGCACTAAAAA
TCATATCATACCCAATTGCGGGCAGCACATGTATATATTATACACTACTTATAACTAACC
TTCTTCATATATAACAATATGCCTGATAATACTAATGCATTTAAATCATATCATGAAAGA
AAATTACATGGGTTTTATTGACATAATTGCATTTAGAATACATAAAATTCTAAAGAATTA
ATATATCCAAAGTATTAGACATAACCAAGAATAATAGTGAATAATTTTAGATTTTGTTAC
ATATAATTCTGCTTGCCTATCTCTTCCACTCTTTTCAAAACGTTGCATGTAAGCGTTACT
AATATTCCGCTTTATTTTGTTGCAATTCCTAATTTTTTCATTACATTATCTTGCGAGTAC
GGAAGCGATTAACGTTCTCCCAATAGAAGGAACAAACATAGATATTGAAGTTTTACTGCT
TTTGCTTACCTGACCTTTTTCAAATTTAATTTTTTCCCGCTAATAAGACCATAAACTACC
CCGAACCAAATTCTAAAAGATAGTCAGCTGGATTAGAGTTGTCATCTCCAAACATTAATT
TTGCATTATCTTCGGCTTCAATCAAATCGCCTGATAAGAACTCCTTTAATTCTTCATGAA
TGTTTGTATGTGGATGACTCTCCATAGTGCCAGCATGATTGTGGTTACCGACCGAATCAT
ATAACGGTGGCTCCCAATTGTGCAATTCAGCCTTACTATTTTGTTCATTCATCAAAGCAT
TTGGGACAGATCTAATATCTATAATTCTTTCCTCACTATTCTCGCTATTATTTTGCCCCG
AACTGGCATGGTGGTTATTGGTAAAAGGAGATAATGCTGCGACAGAACTTTTCTTCTCTT
CCAATTCTTTTATCTTCGCCAATAACGCCTTCTTTTTTCGTGCTTGTATTTCTAAAATTT
CGGCCAGGTATTGTAAATATTCGACCGTTCTATCCAAAATGATACCCTTATTTGGTTTGA
TTTGTTTACCTAGGTCATCGTAATTCAATAAAGATGGTGGAACCAACTGGCCGAGTTCTT
TTATCTTTTGCTTTATTAATTCTCTTCTTCTCCTTTCGACGGCATTATGAAACTCTCTTT
TGCGCCTCAGTTTCTCATCAGAAGTTAACCCGCCTAAAATCTTTGGAACACTTCCAGGTC
CTATATTTTCAGTCATATTGCTACTTATTGAAGTGTGTCTTGTTCTGGGTGTGTTTATGC
TACCATGCCTAAAAGAAGATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATC
TTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAG
AAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAATTGATGCTAGAAGATAAAG
AGGAGCTCAAATCATCAGTTAAGTTATACATTGTGTCATCGGTGCCATGCTGAAATAAAA
AGTCATGTGGTAACTCAGCTTTTAAGCCTTCTTGTTGCTGTAATGGAGTATTCATTGCTC

```

# Answer
ok I have a Python attempt, used `Biopython` only to load your DNA sequence; working only on your sequence ( not calculating the reverse complement) and as per your instruction an ORF is something in between an ATG and a stop codon, so I am using every ATG as a start even if preceded by another ATG, lenght of ORFs is \>30bp meaning \>90 bases. Code:

```
from Bio import SeqIO
# record = SeqIO.read("sample_test.fasta", "fasta")
# record = SeqIO.read("sample_test_2.fasta", "fasta")

record = SeqIO.read("sample.fasta", "fasta")

# print(record.seq)
# print(len(sequence))

sequence = record.seq

orfs = []

for i in range(len(sequence)-2):

    if sequence[i]+sequence[i+1]+sequence[i+2] == 'ATG':

        print('ATG at : ', i)

        orf = 'ATG'
        x = 3
        while sequence[i+x:i+3+x] not in ['TAG', 'TAA', 'TGA'] and i+x < len(sequence):
            orf += sequence[i+x:i+3+x]
            x += 3
        if len(orf) > 90:
            orfs.append((i, orf , len(orf)))

cnt = 1
for i in orfs:

    print('\nsequence number : ', cnt ,'\nstart : ', i[0] ,'\nframe : ',  ((len(sequence)-len(sequence[i[0]:]))%3)+1 ,
          '\nsequence : ', i[1],'\nlenght sequence :', i[2],'\n\n')
    cnt += 1

```

Result, start is with a sequence that goes from 0 to len(sequence)-1:

```
sequence number :  1
start :  2
frame :  3
sequence :  ATGCAATACAGTGACATATATATATACCCTAACACTACCCTAACCCTACCCTATTTCAACCCTTCCAACCTGTCTCTCAACTTACCCTCACATTACCCTACCTCTCCACTTGTTACCCTGTCCCATTCAACCATACCACTCCCAACCACCATCCATCCCTCTACTTACTACCACCAATCAACCGTCCACCATAACCGTTACCCTCCAATTAGCCATATTCAACTTCACTACCACTTACCCTGCCATTACTCTACCATCCACCATCTGCTACTCACCATACTGTTGTTCTACCCTCCATAT
lenght sequence : 300

sequence number :  2
start :  524
frame :  3
sequence :  ATGCACGGCACTTGCCTCAGCGGTCTATACCCTGAGCCATTTACCCATAACTCCCACGATTATCCACATTTTAATATCTATATCTCATTCGGCGGGCCCAAATATTGTATAACTGCTCTTAATACATACGTTATACCACTTTTGCACCATATACTAACCACTCAAATTATATACACTTATGCCAATATAACCAAAAAATTACCACTAAAATCACCTAAACATAAAAATTATTTATCTTTCAACTTTACGAAA
lenght sequence : 252

sequence number :  3
start :  907
frame :  2
sequence :  ATGCCACCTCACTGTCGTAACACTCTATATTCACCGAGAAACAATACGGTAGTGGCTCAAACTCATGCCGGTGCTATGATACAATTGTATCCTATTTCCATTCTCATATGCTATCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCATACGAAGAGACTATATC
lenght sequence : 189

sequence number :  4
start :  982
frame :  2
sequence :  ATGATACAATTGTATCCTATTTCCATTCTCATATGCTATCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCATACGAAGAGACTATATC
lenght sequence : 114

sequence number :  5
start :  1143
frame :  1
sequence :  ATGAAACGTCAAATAACCCTACTATATTATAACTTATCAGCGCCGTATACTAAAACGGACGTTACGATATTGTCTCACTTCATCTTACCACCCTCTATCTTATTGTTGATAAAACAC
lenght sequence : 117

sequence number :  6
start :  1370
frame :  3
sequence :  ATGACTTGTAACTCCACTACCCTGATCTGCAATCTTGTTCTTAGAAGTGACGCATACTCTATATGGCCCGACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGT
lenght sequence : 183

sequence number :  7
start :  1432
frame :  2
sequence :  ATGGCCCGACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTTATTACCGTCCACGCTTTTTCTACAA
lenght sequence : 153

sequence number :  8
start :  1457
frame :  3
sequence :  ATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGT
lenght sequence : 96

sequence number :  9
start :  1485
frame :  1
sequence :  ATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTTATTACCGTCCACGCTTTTTCTACAATAGTG
lenght sequence : 105

sequence number :  10
start :  1830
frame :  1
sequence :  ATGGCTGGCTTTAATCTGCTGGAGTACCATGGAACACCGGTGATCATTCTGGTCACTTGGTCTGGGGCAATACCAGTCAACATGGTGGTGAAGTCACCG
lenght sequence : 99

sequence number :  11
start :  2107
frame :  2
sequence :  ATGGCAGCGACACCAGCAGCGATTGAAGTTAATTTGACCATTGTATTTGTTTTGTTTGTTAGTGCTGATGTAATCTTAACAAGAAATAGTGAAATGAAAGCGCATACCTCAAAGGCATATAGT
lenght sequence : 123

sequence number :  12
start :  2356
frame :  2
sequence :  ATGGAGATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAGGTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght sequence : 303

sequence number :  13
start :  2362
frame :  2
sequence :  ATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAGGTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght sequence : 297

sequence number :  14
start :  2407
frame :  2
sequence :  ATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght sequence : 252

sequence number :  15
start :  2462
frame :  3
sequence :  ATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTT
lenght sequence : 93

sequence number :  16
start :  3447
frame :  1
sequence :  ATGATAATAATATTTATAGAATTGTGTAGAATTGCAGATTCCCTTTTATGGATTCCTAAATCCTGGAGGAGAACTTCTAGTATATCTACATACCTAATATTA
lenght sequence : 102

sequence number :  17
start :  3563
frame :  3
sequence :  ATGGAATCCCAACAATTACATCAAAATCCACATTCTCTACACCAATACCATCGACGAGAGCTTCTAGTAAATTGTATACATAACAGTATAACCCTTACCAACAATGGAATCTCAAAGATTATTAAATTATTCACAGACTCTGAGGATTCGGGTAAAATAGGGTATTTAACTGGTTACCGGAAAGGTTTAGAAAATTCGTGGAGGGTTGGCCGAGTGGTC
lenght sequence : 219

sequence number :  18
start :  4056
frame :  1
sequence :  ATGATGTGTAGCGCATTAGTCTTCTGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGGCTCTACTCTATA
lenght sequence : 96

sequence number :  19
start :  4059
frame :  1
sequence :  ATGTGTAGCGCATTAGTCTTCTGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGGCTCTACTCTATA
lenght sequence : 93

sequence number :  20
start :  5073
frame :  1
sequence :  ATGCTGCGACAGAACTTTTCTTCTCTTCCAATTCTTTTATCTTCGCCAATAACGCCTTCTTTTTTCGTGCTTGTATTTCTAAAATTTCGGCCAGGTATTGTAAATATTCGACCGTTCTATCCAAAA
lenght sequence : 126

sequence number :  21
start :  5253
frame :  1
sequence :  ATGGTGGAACCAACTGGCCGAGTTCTTTTATCTTTTGCTTTATTAATTCTCTTCTTCTCCTTTCGACGGCATTATGAAACTCTCTTTTGCGCCTCAGTTTCTCATCAGAAGTTAACCCGCCTAAAATCTTTGGAACACTTCCAGGTCCTATATTTTCAGTCATATTGCTACTTATTGAAGTGTGTCTTGTTCTGGGTGTGTTTATGCTACCATGCC
lenght sequence : 216

sequence number :  22
start :  5456
frame :  3
sequence :  ATGCTACCATGCCTAAAAGAAGATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTA
lenght sequence : 120

sequence number :  23
start :  5478
frame :  1
sequence :  ATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght sequence : 144

sequence number :  24
start :  5502
frame :  1
sequence :  ATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght sequence : 120

sequence number :  25
start :  5508
frame :  1
sequence :  ATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght sequence : 114

sequence number :  26
start :  5511
frame :  1
sequence :  ATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght sequence : 111

```

I am sure there is some awk, perl, regex way to get the same in much less lines.
> 1 votes

# Answer
using an unwanted (regex library) I can steal code from https://stackoverflow.com/questions/13114246/how-to-find-a-open-reading-frame-in-python , code :

```
import re

from Bio import SeqIO
# record = SeqIO.read("sample_test.fasta", "fasta")
# record = SeqIO.read("sample_test_2.fasta", "fasta")

record = SeqIO.read("sample.fasta", "fasta")

# print(record.seq)
# print(len(sequence))

sequence = str(record.seq)

pattern = re.compile(r'(?=(ATG(?:...)*?)(?=TAG|TGA|TAA))')

revcompseq = sequence[::-1].maketrans("ATGC", "TACG") #reverse complement

# print (pattern.findall(sequence)) #forward search
# print (pattern.findall(sequence[::-1].translate(revcompseq))) #backward search

a = pattern.findall(sequence)

# print(type(a))

cnt = 1
for i in a:
    if len(i) > 90:
        print('\nnumero = ',cnt,'\n',i ,'\nlenght = ', len(i),'\n---------------------------\n\n')
        cnt += 1

```

output :

```
numero =  1
 ATGCAATACAGTGACATATATATATACCCTAACACTACCCTAACCCTACCCTATTTCAACCCTTCCAACCTGTCTCTCAACTTACCCTCACATTACCCTACCTCTCCACTTGTTACCCTGTCCCATTCAACCATACCACTCCCAACCACCATCCATCCCTCTACTTACTACCACCAATCAACCGTCCACCATAACCGTTACCCTCCAATTAGCCATATTCAACTTCACTACCACTTACCCTGCCATTACTCTACCATCCACCATCTGCTACTCACCATACTGTTGTTCTACCCTCCATAT
lenght =  300
---------------------------

numero =  2
 ATGCACGGCACTTGCCTCAGCGGTCTATACCCTGAGCCATTTACCCATAACTCCCACGATTATCCACATTTTAATATCTATATCTCATTCGGCGGGCCCAAATATTGTATAACTGCTCTTAATACATACGTTATACCACTTTTGCACCATATACTAACCACTCAAATTATATACACTTATGCCAATATAACCAAAAAATTACCACTAAAATCACCTAAACATAAAAATTATTTATCTTTCAACTTTACGAAA
lenght =  252
---------------------------

numero =  3
 ATGCCACCTCACTGTCGTAACACTCTATATTCACCGAGAAACAATACGGTAGTGGCTCAAACTCATGCCGGTGCTATGATACAATTGTATCCTATTTCCATTCTCATATGCTATCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCATACGAAGAGACTATATC
lenght =  189
---------------------------

numero =  4
 ATGATACAATTGTATCCTATTTCCATTCTCATATGCTATCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCATACGAAGAGACTATATC
lenght =  114
---------------------------

numero =  5
 ATGAAACGTCAAATAACCCTACTATATTATAACTTATCAGCGCCGTATACTAAAACGGACGTTACGATATTGTCTCACTTCATCTTACCACCCTCTATCTTATTGTTGATAAAACAC
lenght =  117
---------------------------

numero =  6
 ATGACTTGTAACTCCACTACCCTGATCTGCAATCTTGTTCTTAGAAGTGACGCATACTCTATATGGCCCGACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGT
lenght =  183
---------------------------

numero =  7
 ATGGCCCGACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTTATTACCGTCCACGCTTTTTCTACAA
lenght =  153
---------------------------

numero =  8
 ATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGT
lenght =  96
---------------------------

numero =  9
 ATGGAAGGACAAAGTGCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTTATTACCGTCCACGCTTTTTCTACAATAGTG
lenght =  105
---------------------------

numero =  10
 ATGGCTGGCTTTAATCTGCTGGAGTACCATGGAACACCGGTGATCATTCTGGTCACTTGGTCTGGGGCAATACCAGTCAACATGGTGGTGAAGTCACCG
lenght =  99
---------------------------

numero =  11
 ATGGCAGCGACACCAGCAGCGATTGAAGTTAATTTGACCATTGTATTTGTTTTGTTTGTTAGTGCTGATGTAATCTTAACAAGAAATAGTGAAATGAAAGCGCATACCTCAAAGGCATATAGT
lenght =  123
---------------------------

numero =  12
 ATGGAGATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAGGTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght =  303
---------------------------

numero =  13
 ATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAGGTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght =  297
---------------------------

numero =  14
 ATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGGGAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGCTTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTTGTTACTCATTCCTGAAGGT
lenght =  252
---------------------------

numero =  15
 ATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAAGTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTT
lenght =  93
---------------------------

numero =  16
 ATGATAATAATATTTATAGAATTGTGTAGAATTGCAGATTCCCTTTTATGGATTCCTAAATCCTGGAGGAGAACTTCTAGTATATCTACATACCTAATATTA
lenght =  102
---------------------------

numero =  17
 ATGGAATCCCAACAATTACATCAAAATCCACATTCTCTACACCAATACCATCGACGAGAGCTTCTAGTAAATTGTATACATAACAGTATAACCCTTACCAACAATGGAATCTCAAAGATTATTAAATTATTCACAGACTCTGAGGATTCGGGTAAAATAGGGTATTTAACTGGTTACCGGAAAGGTTTAGAAAATTCGTGGAGGGTTGGCCGAGTGGTC
lenght =  219
---------------------------

numero =  18
 ATGATGTGTAGCGCATTAGTCTTCTGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGGCTCTACTCTATA
lenght =  96
---------------------------

numero =  19
 ATGTGTAGCGCATTAGTCTTCTGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGGCTCTACTCTATA
lenght =  93
---------------------------

numero =  20
 ATGCTGCGACAGAACTTTTCTTCTCTTCCAATTCTTTTATCTTCGCCAATAACGCCTTCTTTTTTCGTGCTTGTATTTCTAAAATTTCGGCCAGGTATTGTAAATATTCGACCGTTCTATCCAAAA
lenght =  126
---------------------------

numero =  21
 ATGGTGGAACCAACTGGCCGAGTTCTTTTATCTTTTGCTTTATTAATTCTCTTCTTCTCCTTTCGACGGCATTATGAAACTCTCTTTTGCGCCTCAGTTTCTCATCAGAAGTTAACCCGCCTAAAATCTTTGGAACACTTCCAGGTCCTATATTTTCAGTCATATTGCTACTTATTGAAGTGTGTCTTGTTCTGGGTGTGTTTATGCTACCATGCC
lenght =  216
---------------------------

numero =  22
 ATGCTACCATGCCTAAAAGAAGATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTA
lenght =  120
---------------------------

numero =  23
 ATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght =  144
---------------------------

numero =  24
 ATGACGATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght =  120
---------------------------

numero =  25
 ATGATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght =  114
---------------------------

numero =  26
 ATGGTGATCTTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAGAAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAAT
lenght =  111
---------------------------

```

I can get the start position in my sequence of all the found ORFs with :

```
b = [[m.start(1), m.end(1)] for m in pattern.finditer(sequence)]

print(b, len(b))

```

result a list of tuples :

```
[(2, 302), (316, 400), (373, 400), (404, 479), (421, 505),
(478, 505), (508, 556), (524, 776), (702, 711), (807, 822),
(907, 1096), (971, 983), (982, 1096), (1014, 1032), (1046, 1055),
(1064, 1067), (1117, 1141), (1143, 1260), (1340, 1355), (1354, 1378),
(1370, 1553), (1432, 1585), (1457, 1553), (1485, 1590), (1606, 1663),
 (1626, 1647), (1754, 1769), (1830, 1929),

.........................................

```

but cannot filter for `m.end() - m.start() > 90` in the one liner `b`

any help ??

**EDITED :**

```
b = [[m.start(1), m.end(1)] for m in pattern.finditer(sequence) if (m.end(1) - m.start(1)) > 90]

```

should work only cant get found string with `m.group()`

as an example for first match I get:

`2, 302, <re.Match object; span=(2, 2), match=''>` with `m.group()` = `''`

**EDITED 2:**

I can get the found substring using  `m.groups()[0]`
> 1 votes

# Answer
The code is,

```
from Bio.Seq import Seq
frame = 2
myseq = myseq[frame:]
bioseq = Seq(myseq)
fullaa = bioseq.translate()
for aa in fullaa.split("*"):
    print (aa)

```

The sequence is below the `-----`

The output looks like open readingframe2. The first protein is:

> MQYSDIYIYPNTTLTLPYFNPSNLSLNLPSHYPTSPLVTLSHSTIPLPTTIHPSTYYHQSTVHHNRYPPISHIQLHYHLPCHYSTIHHLLLTILLFYPPY

This meets with the requirements ... the remaining output is below ... it looks like a eukaryotic with introns. It feels like a home-work question.

---

**Edit** Blasted it, its a yeast protein

> RecName: Full=Putative uncharacterized protein YBL109W \[Saccharomyces cerevisiae S288C\] Sequence ID: P89495.3

The hit is 100%

---

Remaining output

```
NVNK
S
IIHIYLPYHFIPPHITCHTHLHLYTDMPYAHGCYSIYHSQTYPTLTFYSTP
PISH
ISTKCTHIIMHGTCLSGLYPEPFTHNSHDYPHFNIYISFGGPKYCITALNTYVIPLLHHILTTQIIYTYANITKKLPLKSPKHKNYLSFNFTK
THSIAYLYHHDVINVKVP
YYHLLERIPFQNISNFHRPYIRIICHLTVVTLYIHRETIR
WLKLMPVL
YNCILFPFSYAIRNILKA
LMHL
SCM
HYSYEETISKEDDTVICTLFL
NYNETSNNPTIL
LISAVY
NGRYDIVSLHLTTLYLIVDKTLTPQLYF
LQLQQTIPNHKS
YFRCQKMRISK
EFGTMTCNSTTLICNLVLRSDAYSIWPDATRQKMKKEAATHFYGRTKCCEVIRFQFHYCLLDILC
LYYRPRFFYNSVKVSFLCSSYS
NASRTPSLIK
VYNINIHLYNRRYLYHQKK
FFILFFHYFSLSMDFHS
RHHSLVCDSVDTVLG
STGDGWL
SAGVPWNTGDHSGHLVWGNTSQHGGEVTVVENGFSNFDWVGFSWVGGLEHVVLGQVSSDIRDVDTQFHQVDSFVRLS
SGGCRKQ
QRWQRHQQRLKLI
PLYLFCLLVLM
S
QEIVK
KRIPQRHIVEAALFIPVPPSVITT
TIR
QTLI
HLTYSPYLIFHTISLSLWRCSCF
TNHTSFIGFVCGVLFYGTHVYSYAQNVGMPIIGVPGGALKNPFARLLSFRFQSKRISEF
IWTLVQKLTV
AIIWCVLNGFHAKQHLCYSFLKVE
KIATMCRESYIVRF
RLVRIIKAVNRINIRKASKVHS
ENS
HVKICVYEMCQGPSIA
LTPLQE
VTFFTLKRQRGHILKSHFQTHLQFATTDNNIMIGWYDVISCHIFLIIALVTSLQLLG
NSKISFL

YSSYQYIITCRKRR
HKDRQTVFKYNGSWNARIDNVTGY
MTKYK
KKKKKKSSNTIMWKYRFHFEDSNCWNKNQLSSTN
YLRY
YIIIYGVRR
RK

EIVI
ISGS
NARIDNVIGSMNINI
NDDNNIYRIV
NCRFPFMDS
ILEENF
YIYIPNIIALITMESQQLHQNPHSLHQYHRRELLVNCIHNSITLTNNGISKIIKLFTDSEDSGKIGYLTGYRKGLENSWRVGRVV
GGRLKICWTVVRASSNLAFIFRIHFIVGHFISCRVH
KIINQIII
PMGHFFDETHVYSSIRS
YCTYDYQ
HQNNHRAYAILKYIIAVKYDV
RISLLHSYYQWYTVWYMLSNHEVLALLYINSTSIDFLNK
VTQVSFPLIQVIYILATPNYKALKIISYPIAGSTCIYYTLLITNLLHI
QYA

Y
CI
IIS
KKITWVLLT
LHLEYIKF
RINISKVLDITKNNSE
F
ILLHIILLAYLFHSFQNVACKRY
YSALFCCNS
FFHYIILRVRKRLTFSQ
KEQT
ILKFYCFCLPDLFQI
FFPANKTINYPEPNSKR
SAGLELSSPNINFALSSASIKSPDKNSFNSS
MFVCG
LSIVPA
LWLPTESYNGGSQLCNSALLFCSFIKAFGTDLISIILSSLFSLLFCPELAWWLLVKGDNAATELFFSSNSFIFANNAFFFRACISKISARYCKYSTVLSKMIPLFGLICLPRSS
FNKDGGTNWPSSFIFCFINSLLLLSTAL
NSLLRLSFSSEVNPPKIFGTLPGPIFSVILLLIEVCLVLGVFMLPCLKEDERKLPARNDDDGDLTFGEYVEDTDAGPKLCGL
ENDDEYVFGVIISELMLEDKEELKSSVKLYIVSSVPC
NKKSCGNSAFKPSCCCNGVFIA

```

---

```
myseq = "GCATGCAATACAGTGACATATATATATACCCTAACACTACCCTAACCCTACCCTATTTCA\
ACCCTTCCAACCTGTCTCTCAACTTACCCTCACATTACCCTACCTCTCCACTTGTTACCC\
TGTCCCATTCAACCATACCACTCCCAACCACCATCCATCCCTCTACTTACTACCACCAAT\
CAACCGTCCACCATAACCGTTACCCTCCAATTAGCCATATTCAACTTCACTACCACTTAC\
CCTGCCATTACTCTACCATCCACCATCTGCTACTCACCATACTGTTGTTCTACCCTCCAT\
ATTGAAACGTTAACAAATGATCGTAAATAATACACATATACTTACCCTACCACTTTATAC\
CACCACACATCACATGCCATACTCACCTTCACTTGTATACTGATATGCCATACGCACACG\
GATGCTACAGTATATACCACTCTCAAACTTACCCTACTCTCACATTCTACTCCACTCCAT\
GACCCATCTCTCACTAAATCAGTACTAAATGCACCCACATCATTATGCACGGCACTTGCC\
TCAGCGGTCTATACCCTGAGCCATTTACCCATAACTCCCACGATTATCCACATTTTAATA\
TCTATATCTCATTCGGCGGGCCCAAATATTGTATAACTGCTCTTAATACATACGTTATAC\
CACTTTTGCACCATATACTAACCACTCAAATTATATACACTTATGCCAATATAACCAAAA\
AATTACCACTAAAATCACCTAAACATAAAAATTATTTATCTTTCAACTTTACGAAATAAA\
CACACTCAATTGCGTATCTATACCACCATGACGTCATTAACGTAAAAGTTCCTTAATATT\
ACCATTTGCTTGAACGGATACCATTTCAGAATATTTCTAACTTTCACAGACCATACATTA\
GAATAATATGCCACCTCACTGTCGTAACACTCTATATTCACCGAGAAACAATACGGTAGT\
GGCTCAAACTCATGCCGGTGCTATGATACAATTGTATCCTATTTCCATTCTCATATGCTA\
TCCGCAATATCCTAAAAGCATAACTGATGCATCTTTAATCTTGTATGTGACACTACTCAT\
ACGAAGAGACTATATCTAAAGAAGACGATACAGTGATATGTACGTTGTTTTTGTAGAATT\
ATAATGAAACGTCAAATAACCCTACTATATTATAACTTATCAGCGCCGTATACTAAAACG\
GACGTTACGATATTGTCTCACTTCATCTTACCACCCTCTATCTTATTGTTGATAAAACAC\
TAACCCCTCAGCTTTATTTCTAGTTACAGTTACAACAAACTATCCCAAACCATAAATCTT\
AATATTTTAGGTGTCAAAAAATGAGGATCTCCAAATGAGAGTTTGGTACCATGACTTGTA\
ACTCCACTACCCTGATCTGCAATCTTGTTCTTAGAAGTGACGCATACTCTATATGGCCCG\
ACGCGACGCGCCAAAAAATGAAAAAAGAAGCAGCGACTCATTTTTATGGAAGGACAAAGT\
GCTGCGAAGTCATACGCTTCCAATTTCATTATTGTTTATTGGACATACTCTGTTAGCTTT\
ATTACCGTCCACGCTTTTTCTACAATAGTGTAAAAGTTTCTTTCTTATGTTCATCGTATT\
CATAAAATGCTTCACGAACACCGTCATTGATCAAATAGGTTTATAATATTAATATACATT\
TATATAATCGGCGGTATTTATATCATCAAAAAAAGTAGTTTTTTATTTTATTTTTTCATT\
ACTTTTCACTGTCTATGGATTTTCATTCGTAAAGGCATCACTCCCTAGTTTGCGATAGTG\
TAGATACCGTCCTTGGATAGAGCACTGGAGATGGCTGGCTTTAATCTGCTGGAGTACCAT\
GGAACACCGGTGATCATTCTGGTCACTTGGTCTGGGGCAATACCAGTCAACATGGTGGTG\
AAGTCACCGTAGTTGAAAACGGCTTCAGCAACTTCGACTGGGTAGGTTTCAGTTGGGTGG\
GCGGCTTGGAACATGTAGTATTGGGCCAAGTGAGCTCTGATATCAGAGACGTAGACACCC\
AATTCCACCAAGTTGACTCTTTCGTCAGATTGAGCTAGAGTGGTGGTTGCAGAAAGCAGT\
AGCAGCGATGGCAGCGACACCAGCAGCGATTGAAGTTAATTTGACCATTGTATTTGTTTT\
GTTTGTTAGTGCTGATGTAATCTTAACAAGAAATAGTGAAATGAAAGCGCATACCTCAAA\
GGCATATAGTTGAAGCAGCTCTATTTATACCCGTTCCTCCATCTGTCATCACTACTTAAA\
CGATTCGTTAACAGACGCTCATTTAGCACCTCACATATTCTCCATATCTCATCTTTCACA\
CAATCTCATTATCTCTATGGAGATGCTCTTGTTTCTGAACGAATCATACATCTTTCATAG\
GTTTCGTATGTGGAGTATTGTTTTATGGCACTCATGTGTATTCGTATGCGCAGAATGTGG\
GAATGCCTATTATAGGGGTGCCGGGGGGTGCCTTGAAAAACCCTTTTGCGCGCCTGTTAA\
GTTTCCGTTTTCAGTCAAAAAGAATATCCGAATTTTAGATTTGGACCCTCGTTCAGAAGC\
TTACTGTCTAAGCCATCATTTGGTGTGTCCTAAACGGTTTCCATGCAAAGCAACATCTTT\
GTTACTCATTCCTGAAGGTTGAATAAAAAATAGCAACCATGTGCAGGGAGTCGTATATTG\
TTAGATTCTAGAGACTTGTACGCATCATCAAAGCTGTAAATAGAATAAACATACGCAAGG\
CGTCAAAAGTGCATAGTTAAGAAAATTCCTGACATGTGAAAATATGTGTTTATGAAATGT\
GTCAAGGCCCGTCTATAGCGTAGTTAACCCCTCTGCAGGAGTAAGTGACTTTTTTTACGC\
TCAAAAGGCAACGAGGGCACATACTTAAAAGTCATTTTCAAACACATCTGCAGTTTGCAA\
CGACAGATAACAATATTATGATAGGATGGTATGATGTTATTAGCTGCCACATATTTTTAA\
TAATAGCATTAGTCACGTCTCTTCAATTGTTGGGATGAAACTCTAAAATATCATTCCTTT\
AGTAGTATTCCAGTTACCAGTATATTATCACATGCCGAAAAAGAAGATGACATAAAGATC\
GACAAACAGTCTTCAAATATAATGGAAGCTGGAATGCAAGGATTGATAATGTAACAGGAT\
ACTGAATGACAAAGTATAAATGAAAAAAAAAAAAAAAAAAAAGTAGTAATACTATTATGT\
GGAAATACCGATTCCATTTTGAGGATTCCAATTGTTGGAATAAAAATCAACTATCATCTA\
CTAACTAGTATTTACGTTACTAGTATATTATCATATACGGTGTTAGAAGATGACGCAAAT\
GATGAGAAATAGTCATCTAAATTAGTGGAAGCTGAAACGCAAGGATTGATAATGTAATAG\
GATCAATGAATATTAACATATAAAACGATGATAATAATATTTATAGAATTGTGTAGAATT\
GCAGATTCCCTTTTATGGATTCCTAAATCCTGGAGGAGAACTTCTAGTATATCTACATAC\
CTAATATTATAGCCTTAATCACAATGGAATCCCAACAATTACATCAAAATCCACATTCTC\
TACACCAATACCATCGACGAGAGCTTCTAGTAAATTGTATACATAACAGTATAACCCTTA\
CCAACAATGGAATCTCAAAGATTATTAAATTATTCACAGACTCTGAGGATTCGGGTAAAA\
TAGGGTATTTAACTGGTTACCGGAAAGGTTTAGAAAATTCGTGGAGGGTTGGCCGAGTGG\
TCTAAGGCGGCAGACTTAAGATCTGTTGGACGGTTGTCCGCGCGAGTTCGAACCTCGCAT\
TCATTTTCCGGATCCACTTTATAGTGGGTCACTTCATCTCTTGTAGGGTTCACTGAAAGA\
TAATCAATCAAATTATCATATGACCAATGGGTCACTTTTTTGATGAGACCCATGTTTATT\
CTTCTATACGTTCGTGATACTGCACTTACGACTACCAGTAACATCAGAATAACCATAGAG\
CTTATGCTATCTTGAAATATATCATCGCCGTAAAGTATGATGTGTAGCGCATTAGTCTTC\
TGCATTCATATTATCAGTGGTATACCGTTTGGTATATGCTTTCCAATCATGAAGTGTTGG\
CTCTACTCTATATAAATTCGACATCCATTGATTTTCTCAATAAATGAGTTACCCAAGTAA\
GCTTTCCATTGATACAAGTGATCTACATTCTTGCGACGCCAAATTATAAAGCACTAAAAA\
TCATATCATACCCAATTGCGGGCAGCACATGTATATATTATACACTACTTATAACTAACC\
TTCTTCATATATAACAATATGCCTGATAATACTAATGCATTTAAATCATATCATGAAAGA\
AAATTACATGGGTTTTATTGACATAATTGCATTTAGAATACATAAAATTCTAAAGAATTA\
ATATATCCAAAGTATTAGACATAACCAAGAATAATAGTGAATAATTTTAGATTTTGTTAC\
ATATAATTCTGCTTGCCTATCTCTTCCACTCTTTTCAAAACGTTGCATGTAAGCGTTACT\
AATATTCCGCTTTATTTTGTTGCAATTCCTAATTTTTTCATTACATTATCTTGCGAGTAC\
GGAAGCGATTAACGTTCTCCCAATAGAAGGAACAAACATAGATATTGAAGTTTTACTGCT\
TTTGCTTACCTGACCTTTTTCAAATTTAATTTTTTCCCGCTAATAAGACCATAAACTACC\
CCGAACCAAATTCTAAAAGATAGTCAGCTGGATTAGAGTTGTCATCTCCAAACATTAATT\
TTGCATTATCTTCGGCTTCAATCAAATCGCCTGATAAGAACTCCTTTAATTCTTCATGAA\
TGTTTGTATGTGGATGACTCTCCATAGTGCCAGCATGATTGTGGTTACCGACCGAATCAT\
ATAACGGTGGCTCCCAATTGTGCAATTCAGCCTTACTATTTTGTTCATTCATCAAAGCAT\
TTGGGACAGATCTAATATCTATAATTCTTTCCTCACTATTCTCGCTATTATTTTGCCCCG\
AACTGGCATGGTGGTTATTGGTAAAAGGAGATAATGCTGCGACAGAACTTTTCTTCTCTT\
CCAATTCTTTTATCTTCGCCAATAACGCCTTCTTTTTTCGTGCTTGTATTTCTAAAATTT\
CGGCCAGGTATTGTAAATATTCGACCGTTCTATCCAAAATGATACCCTTATTTGGTTTGA\
TTTGTTTACCTAGGTCATCGTAATTCAATAAAGATGGTGGAACCAACTGGCCGAGTTCTT\
TTATCTTTTGCTTTATTAATTCTCTTCTTCTCCTTTCGACGGCATTATGAAACTCTCTTT\
TGCGCCTCAGTTTCTCATCAGAAGTTAACCCGCCTAAAATCTTTGGAACACTTCCAGGTC\
CTATATTTTCAGTCATATTGCTACTTATTGAAGTGTGTCTTGTTCTGGGTGTGTTTATGC\
TACCATGCCTAAAAGAAGATGAAAGGAAACTTCCTGCACGAAATGACGATGATGGTGATC\
TTACCTTTGGGGAATATGTGGAAGATACGGATGCTGGGCCCAAGCTTTGTGGATTATAAG\
AAAATGATGATGAATATGTGTTTGGTGTCATCATATCAGAATTGATGCTAGAAGATAAAG\
AGGAGCTCAAATCATCAGTTAAGTTATACATTGTGTCATCGGTGCCATGCTGAAATAAAA\
AGTCATGTGGTAACTCAGCTTTTAAGCCTTCTTGTTGCTGTAATGGAGTATTCATTGCTC"

```
> 0 votes

---
Tags: python, homework
---


# Question
Title: Error while using SummarizedExperiment() in R
I'm tryig to perform RKM normalization on raw counts for RNA-Seq Data:

```
counts_data <- exprs181517[,-1]
#Setting title col to rownames for pheno data
rownames(pheno.data181517) <- NULL
rownames(counts_data) <- exprs181517[,1]
phenoColData <- tibble::column_to_rownames(pheno.data181517, "title")
colnames(counts_data) <- rownames(phenoColData)
#Defining the design column
designCol <- as.data.frame(phenoColData$characteristics_ch1)
rownames(designCol) <- rownames(phenoColData)

#Creating a summarized experiment
se <- SummarizedExperiment(list(counts=counts_data), colData=phenoColData)

DESeqDataSetFromMatrix(se, design ~ designCol)

```

but I keep getting this error:

```
Error in h(simpleError(msg, call)) :
  error in evaluating the argument 'x' in selecting a method for function 'as.matrix':
    no method for coercing this S4 class to a vector

```

# Answer
> 2 votes
DESeqDataSetFrom**Matrix** but a SummarizedExperiment is not a matrix. Use the function `DESeqDataSet()` instead and make sure the `design` argument is followed by a `=`.

---
Tags: rna-seq, bioconductor, normalization
---


title: BLAST(n): No hits found

text:
# Question
Title: BLAST(n): No hits found
I am currently exploring the BLAST program, just for testing purpose i generated two FASTA files, containing two genes A and B, such that B is just a motif of repeated 'G's that occurs in A.

file A.fna:

```
>A
ACCAACAGACGGGGGGGGGGGGGAGAAAACCGGCCCCCGCAGGGACGCCACCCAGAGGGAAGGG

```

file B.fna:

```
>B
GGGGGGGGGGGG

```

So by my understanding BLAST should find the local alignment of B in A easily? But the program returns 'No hits found'? Maybe I am doing something wrong here?

I executed the following:

```
makeblastdb -in A.fna -dbtype 'nucl'
makeblastdb -in B.fna -dbtype 'nucl'
blastn -db A.fna -query B.fna

```

output of blastn:

```
BLASTN 2.7.1+

Reference: Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb Miller (2000), "A greedy algorithm for aligning DNA sequences", J Comput Biol 2000; 7(1-2):203-14.

Database: A.fna
           1 sequences; 64 total letters

Query= B

Length=12

***** No hits found *****

Lambda      K        H
    1.33    0.621     1.12

Gapped Lambda      K        H
    1.28    0.460    0.850

Effective search space used: 480

  Database: A.fna
    Posted date:  May 3, 2018  5:17 PM   Number of letters in database: 64   Number of sequences in database:  1

Matrix: blastn matrix 1 -2 Gap Penalties: Existence: 0, Extension: 2.5

```

Thanks for any suggestions (:

# Answer
There are three possible problems that come to mind.

1. Masking

Blast will mask low complexity regions by default. Since your sequence is nothing but Gs, it is a safe bet that it is being masked, so no hits will be found for it.
2. Score/e-value thresholds

Another source of complication is that even if a match is found, that match will have very bad scores. Both the actual score of the alignment and the e-value will be very bad. Since this is such a simple sequence, it will always score badly.
3. Word-size

The way blast works is (simplifying a bit) by finding a match for N residues (the word size) and then extending that match if extending increases the score. If your query sequence is shorter than the word size, no match will ever be found.

With all that said, the specific problem in your case seems to be a combination of the word-size and masking. If you run your query setting the word-size to 11 and disabling masking, as shown below, you should get results.

```
blastn -db A.fna -query B.fna -word_size 11 -dust no

```

That should return two hits (two because the word size is one less than the query size, so it can align perfectly at two locations):

```
> A
Length=64

 Score = 23.3 bits (12),  Expect = 5e-05
 Identities = 12/12 (100%), Gaps = 0/12 (0%)
 Strand=Plus/Plus

Query  1   GGGGGGGGGGGG  12
           ||||||||||||
Sbjct  11  GGGGGGGGGGGG  22

 Score = 23.3 bits (12),  Expect = 5e-05
 Identities = 12/12 (100%), Gaps = 0/12 (0%)
 Strand=Plus/Plus

Query  1   GGGGGGGGGGGG  12
           ||||||||||||
Sbjct  12  GGGGGGGGGGGG  23

```

Note that the default word size for BLASTn searches is usually 11. I don't know why the blast+ version (2.2.31+) I have on my Ubuntu machine didn't use that as a default. The manual says it should. But I got no hits when I didn't explicitly set the word size to 11.
> 13 votes

# Answer
I came across this thread while searching for a solution to a similar problem I was experiencing. I am inputting markers that are 20 basepairs in length as the query sequence and haven't been getting hits. As terdon mentioned in their excellent reply, setting word size to something smaller like 11 helps for these shorter query sequences.

I didn't get any results on my 20 bp query sequences when using blastn against the sunflower genome, but I got plenty when setting word size to 11. So, why does setting word size for blastn to 11 improve the results despite documentation stating that default word size is 11? Well it turns out that this is very misleading. This webpage explains that the default **algorithm** used in blastn is actually megablast, not blastn. The default word length for this default megablast algorithm is actually 28!

The documentation is quite unclear about this and many websites providing tutorials for blastn also seem to miss this key detail. Hope this helps anyone else who struggles with this!
> 3 votes

---
Tags: sequence-alignment, blast
---


# Question
Title: Is there a smaller representative version of refseq with fewer than 30K files in it?
Does any one have a condensed NCBI reference seq database with all the gbff files, which covers mostly all the bacteria and archaea. The current NCBI ReqSeq contains more than 30,000 gbff files. I am looking for a much more curated refseq with less than 30K files in it.

# Answer
proGenomes includes a collection of genomes that are representative of different types - https://progenomes.embl.de/download.cgi

Or Genome Updater could be used to filter on unique taxonomy IDs and genome completeness etc.

https://github.com/pirovc/genome\_updater
> 1 votes

---
Tags: bacteria, refseq
---


# Question
Title: How to set Entrez Direct - ESearch NCBI entrez to pipe or ignore errors?
I am using the following one-liner with `NCBI entrez` to query their databases from the terminal (see Entrez Direct: E-utilities on the Unix Command Line:

```
esearch -db assembly -query ${species name} |
    xtract -pattern ENTREZ_DIRECT -element Count

```

This works a charm but occasionally it sends an xml error that messes up the output.

Here is an example of the error:

```
xml
FAILURE ( Fri Jan 27 01:12:44 PM CST 2023 )
nquire -url https://eutils.ncbi.nlm.nih.gov/entrez/eutils/ esearch.fcgi -retmax 0 -usehistory y -db assembly -term "Crepis commutata" -tool edirect -edirect 18.7 -edirect_os Linux -email {}
<ErrorList>
  <PhraseNotFound>commutata</PhraseNotFound>
  <PhraseNotFound>Crepis commutata[All Fields]</PhraseNotFound>
<ErrorList>
SECOND ATTEMPT

```

Entrez, unfortunately, does not have a man page or `--help` entry on my installation for some reason and I am not quite sure where I should change the one-liner to change this behavior.

Does anyone have any clues?

# Answer
This is because there are not assemblies for *Crepis commutate* within NCBI's assembly database.

The cross check is here:

https://www.ncbi.nlm.nih.gov/assembly/?term=Crepis+commutata

The output is a bit clearer than Entrez

> The following terms were not found in Assembly: commutata, Crepis commutata\[All Fields\]. No items found.

All you need do to correct it is issue a `next` statement in the a loop of species above the code you've presented.

```
if [[ "$species" == 'Crepis commutata' ]]; then
  continue
fi
esearch -db assembly -query "$species" | xtract -pattern ENTREZ_DIRECT -element Count

```

---

To answer the comments if this is done via a Python subprocess, you could simply issue the following code (which has not been checked)

```
import subprocess

for species in List:
    try:
        subprocess.run(["esearch", "-db" "assembly", "-query", f"${species}", "|", "xtract", "-pattern", "ENTREZ_DIRECT", "-", "element"], shell=True")
    except ValueError as ve:
        print(f'Species {species}, is not present.')
        continue

```
> 1 votes

---
Tags: database, entrez, efetch, ncbi
---


# Question
Title: How is the "canonical" version (`_entity_poly.pdbx_seq_one_letter_code`) obtained in the PDB?
I encountered the dichotomy in the context of PDBx/mmcif files, say, `6OSQ`: each chain has a `_entity_poly.pdbx_seq_one_letter_code` but also `_entity_poly.pdbx_seq_one_letter_code_can` which i assume stands for "canonical". The *raw* and the *canonical* sequences.

The difference seems to be that all the ligands, hetatoms and other junk (`(PSU)`, `(MP4)`, `(MG4)`) have been removed from the the raw sequence to make it canonical i.e. it can now be packaged as `fasta`.

My question is where does one find documentation for this process? Is making it amenable to fasta tools the primary motivation? What kind of hetero atoms *can* occur there? Where is the logic written down for replacing them or deleting them when deriving the canonical sequence? Is the job of the structure's depositor or the PDB? Is there a tool for this?

Below is chain `1` of PDB`6OSQ`. First is the raw seq, then is the canonical one.

```
1 polyribonucleotide no yes
;GGUUAAGCGACUAAGCGUACACGGUGGAUGCCCUGGCAGUCAGAGGCGAUGAAGGACGUGCUAAUCUGCGAUAAGCGUCG
GUAAGGUGAUAUGAACCGUUAUAACCGGCGAUUUCCGAAUGGGGAAACCCAGUGUGUUUCGACACACUAUCAUUAACUGA
AUCCAUAGGUUAAUGAGGCGAACCGGGGGAACUGAAACAUCUAAGUACCCCGAGGAAAAGAAAUCAACCGAGAUUCCCCC
AGUAGCGGCGAGCGAACGGGGAGCAGCCCAGAGCCUGAAUCAGUGUGUGUGUUAGUGGAAGCGUCUGGAAAGGCGCGCGA
UACAGGGUGACAGCCCCGUACACAAAAAUGCACAUGCUGUGAGCUCGAUGAGUAGGGCGGGACACGUGGUAUCCUGUCUG
AAUAUGGGGGGACCAUCCUCCAAGGCUAAAUACUCCUGACUGACCGAUAGUGAACCAGUACCGUGAGGGAAAGGCGAAAA
GAACCCCGGCGAGGGGAGUGAAAAAGAACCUGAAACCGUGUACGUACAAGCAGUGGGAGCACGCUUAGGCGUGUGACUGC
GUACCUUUUGUAUAAUGGGUCAGCGACUUAUAUUCUGUAGCAAGGUUAACCGAAUAGGGGAGCCGAAGGGAAACCGAGUC
UUAACUGGGCGUUAAGUUGCAGGGUAUAGACCCGAAACCCGGUGAUCUAGCCAUGGGCAGGUUGAAGGUUGGGUAACACU
AACUGGAGGACCGAACCGACUAAU(1MG)(PSU)(5MU)GAAAAAUUAGCGGAUGACUUGUGGCUGGGGGUGAAAGGCCA
AUCAAACCGGGAGAUAGCUGGUUCUCCCCGAAAGCUAUUUAGGUAGCGCCUCGUGAAUUCAUCUCCGGGGGUAGAGCACU
GUUUCGGCAAGGGGGUCAACCCGACUUACCAACCCGAUGCAAACUGCGAAUACCGGAGAAUGUUAUCACGGGAGACACAC
GGCGGG(PSU)GCUAACGUCCGUCGUGAAGAGGGAAACAACCCAGACCGCCAGCUAAGGUCCCAAAGUCAUGGUUAAGUG
GGAAACGAUGUGGGAAGGCCCAGACAGCCAGGAUGUUGGCUUAGAAGCAGCCAUCAUUUAAAGAAAGCGUAAUAGCUCAC
UGGUCGAGUCGGCCUGCGCGGAAGAUGUAACGGGGCUAAACCAUGCACCGAAGCUGCGGCAGCGACGCUUAUGCGUUGUU
GGGUAGGGGAGCGUUCUGUAAGCCUGCGAAGGUGUGCUGUGAGGCAUGCUGGAGGUAUCAGAAGUGCGAAUGCUGACAUA
AGUAACGAUAAAGCGGGUGAAAAGCCCGCUCGCCGGAAGACCAAGGGUUCCUGUCCAACGUUAAUCGGGGCAGGGUGAGU
CGACCCCUAAGGCGAGGCCGAAAGGCGUAGUCGAUGGGAAACAGGUUAAUAUUCCUGUACUUGGUGUUACUGCGAAGGGG
GGACGGAGAAGGCUAUGUUGGCCGGGCGACGGUUGUCCCGGUUUAAGCGUGUAGGCUGGUUUUCCAGGCAAAUCCGGAAA
AUCAAGGCUGAGGCGUGAUGACGAGGCACUACGGUGCUGAAGCAACAAAUGCCCUGCUUCCAGGAAAAGCCUCUAAGCAU
CAGGUAACAUCAAAUCGUACCCCAAACCGACAC(6MZ)GGUGGUCAGGUAGAGAAUACCAAGGCGCUUGAGAGAACUCGG
GUGAAGGAACUAGGCAAAAUGGUGCCGUAACUUCGGGAGAAGGCACGCUGAUAUGUAGGUGAGGUCCCUCGCGGAUGGAG
CUGAAAUCAGUCGAAGAUACCAGCUGGCUGCAACUGUUUAUUAAAAACACAGCACUGUGCAAACACGAAAGUGGACGUAU
ACGGUGUGACGCCU(2MG)CCCGGUGCCGGAAGGUUAAUUGAUGGGGUUAGCGCAAGCGAAGCUCUUGAUCGAAGCCCCG
GUAAACGGCGGCCG(PSU)AAC(3TD)A(PSU)AACGGUCCUAAGGUAGCGAAA(5MU)UCCUUGUCGGGUAAGUUCCGA
C(5MC)UGCACGAAUGGCGUAAUGAUGGCCAGGCUGUCUCCACCCGAGACUCAGUGAAAUUGAACUCGCUGUG(6MZ)AG
AUGCAGUGUACCCGCGGCAAGACGGAAAGACCCCGU(G7M)AACCUUUACUAUAGCUUGACACUGAACAUUGAGCCUUGA
UGUGUAGGAUAGGUGGGAGGCUUUGAAGUGUGGACGCCAGUCUGCAUGGAGCCGACCUUGAAAUACCACCCUUUAAUGUU
UGAUGUUCUAACGUUGACCCGUAAUCCGGGUUGCGGACAGUGUCUGGUGGGUAGUUUGACUG(OMG)GGCGGUCUCCUCC
UAAAGAGUAACGGAGGAGCACGAAGGUUGGCUAAUCCUGGUCGGACAUCAGGAGGUUAGUGCAAUGGCAUAAGCCAGCUU
GACUGCGAGCGUGACGGCGCGAGCAGGUGCGAAAGCAGGUCAUAGUGAUCCGGUGGUUCUGAAUGGAAGGGCCAUCGCUC
AACGGAUAAAAGGUACUCCG(2MG)GGAUAACAGGC(PSU)GAUACCGCCCAAGAGUUCAUAUCGACGGCGGUGUUUGGC
A(OMC)CUCG(2MA)(PSU)GUCGGCUCAUCACAUCCUGGGGCUGAAGUAGGUCCCAAGGGUAUGGC(OMU)GUUCGCCA
UUUAAAGUGGUACGCGAGC(PSU)GGGUUUAGAACGUCGUGAGACAGU(PSU)CGGUCCCUAUCUGCCGUGGGCGCUGGA
GAACUGAGGGGGGCUGCUCCUAGUACGAGAGGACCGGAGUGGACGCAUCACUGGUGUUCGGGUUGUCAUGCCAAUGGCAC
UGCCCGGUAGCUAAAUGCGGAAGAGAUAAGUGCUGAAAGCAUCUAAGCACGAAACUUGCCCCGAGAUGAGUUCUCCCUGA
CCCUUUAAGGGUCCUGAAGGAACGUUGAAGACGACGACGUUGAUAGGCCGGGUGUGUAAGCGCAGCGAUGCGUUGAGCUA
ACCGGUACUAAUGAACCGUGAGGCUUAACCU
;
;GGUUAAGCGACUAAGCGUACACGGUGGAUGCCCUGGCAGUCAGAGGCGAUGAAGGACGUGCUAAUCUGCGAUAAGCGUCG
GUAAGGUGAUAUGAACCGUUAUAACCGGCGAUUUCCGAAUGGGGAAACCCAGUGUGUUUCGACACACUAUCAUUAACUGA
AUCCAUAGGUUAAUGAGGCGAACCGGGGGAACUGAAACAUCUAAGUACCCCGAGGAAAAGAAAUCAACCGAGAUUCCCCC
AGUAGCGGCGAGCGAACGGGGAGCAGCCCAGAGCCUGAAUCAGUGUGUGUGUUAGUGGAAGCGUCUGGAAAGGCGCGCGA
UACAGGGUGACAGCCCCGUACACAAAAAUGCACAUGCUGUGAGCUCGAUGAGUAGGGCGGGACACGUGGUAUCCUGUCUG
AAUAUGGGGGGACCAUCCUCCAAGGCUAAAUACUCCUGACUGACCGAUAGUGAACCAGUACCGUGAGGGAAAGGCGAAAA
GAACCCCGGCGAGGGGAGUGAAAAAGAACCUGAAACCGUGUACGUACAAGCAGUGGGAGCACGCUUAGGCGUGUGACUGC
GUACCUUUUGUAUAAUGGGUCAGCGACUUAUAUUCUGUAGCAAGGUUAACCGAAUAGGGGAGCCGAAGGGAAACCGAGUC
UUAACUGGGCGUUAAGUUGCAGGGUAUAGACCCGAAACCCGGUGAUCUAGCCAUGGGCAGGUUGAAGGUUGGGUAACACU
AACUGGAGGACCGAACCGACUAAUGUUGAAAAAUUAGCGGAUGACUUGUGGCUGGGGGUGAAAGGCCAAUCAAACCGGGA
GAUAGCUGGUUCUCCCCGAAAGCUAUUUAGGUAGCGCCUCGUGAAUUCAUCUCCGGGGGUAGAGCACUGUUUCGGCAAGG
GGGUCAACCCGACUUACCAACCCGAUGCAAACUGCGAAUACCGGAGAAUGUUAUCACGGGAGACACACGGCGGGUGCUAA
CGUCCGUCGUGAAGAGGGAAACAACCCAGACCGCCAGCUAAGGUCCCAAAGUCAUGGUUAAGUGGGAAACGAUGUGGGAA
GGCCCAGACAGCCAGGAUGUUGGCUUAGAAGCAGCCAUCAUUUAAAGAAAGCGUAAUAGCUCACUGGUCGAGUCGGCCUG
CGCGGAAGAUGUAACGGGGCUAAACCAUGCACCGAAGCUGCGGCAGCGACGCUUAUGCGUUGUUGGGUAGGGGAGCGUUC
UGUAAGCCUGCGAAGGUGUGCUGUGAGGCAUGCUGGAGGUAUCAGAAGUGCGAAUGCUGACAUAAGUAACGAUAAAGCGG
GUGAAAAGCCCGCUCGCCGGAAGACCAAGGGUUCCUGUCCAACGUUAAUCGGGGCAGGGUGAGUCGACCCCUAAGGCGAG
GCCGAAAGGCGUAGUCGAUGGGAAACAGGUUAAUAUUCCUGUACUUGGUGUUACUGCGAAGGGGGGACGGAGAAGGCUAU
GUUGGCCGGGCGACGGUUGUCCCGGUUUAAGCGUGUAGGCUGGUUUUCCAGGCAAAUCCGGAAAAUCAAGGCUGAGGCGU
GAUGACGAGGCACUACGGUGCUGAAGCAACAAAUGCCCUGCUUCCAGGAAAAGCCUCUAAGCAUCAGGUAACAUCAAAUC
GUACCCCAAACCGACACAGGUGGUCAGGUAGAGAAUACCAAGGCGCUUGAGAGAACUCGGGUGAAGGAACUAGGCAAAAU
GGUGCCGUAACUUCGGGAGAAGGCACGCUGAUAUGUAGGUGAGGUCCCUCGCGGAUGGAGCUGAAAUCAGUCGAAGAUAC
CAGCUGGCUGCAACUGUUUAUUAAAAACACAGCACUGUGCAAACACGAAAGUGGACGUAUACGGUGUGACGCCUGCCCGG
UGCCGGAAGGUUAAUUGAUGGGGUUAGCGCAAGCGAAGCUCUUGAUCGAAGCCCCGGUAAACGGCGGCCGUAACXAUAAC
GGUCCUAAGGUAGCGAAAUUCCUUGUCGGGUAAGUUCCGACCUGCACGAAUGGCGUAAUGAUGGCCAGGCUGUCUCCACC
CGAGACUCAGUGAAAUUGAACUCGCUGUGAAGAUGCAGUGUACCCGCGGCAAGACGGAAAGACCCCGUGAACCUUUACUA
UAGCUUGACACUGAACAUUGAGCCUUGAUGUGUAGGAUAGGUGGGAGGCUUUGAAGUGUGGACGCCAGUCUGCAUGGAGC
CGACCUUGAAAUACCACCCUUUAAUGUUUGAUGUUCUAACGUUGACCCGUAAUCCGGGUUGCGGACAGUGUCUGGUGGGU
AGUUUGACUGGGGCGGUCUCCUCCUAAAGAGUAACGGAGGAGCACGAAGGUUGGCUAAUCCUGGUCGGACAUCAGGAGGU
UAGUGCAAUGGCAUAAGCCAGCUUGACUGCGAGCGUGACGGCGCGAGCAGGUGCGAAAGCAGGUCAUAGUGAUCCGGUGG
UUCUGAAUGGAAGGGCCAUCGCUCAACGGAUAAAAGGUACUCCGGGGAUAACAGGCUGAUACCGCCCAAGAGUUCAUAUC
GACGGCGGUGUUUGGCACCUCGAUGUCGGCUCAUCACAUCCUGGGGCUGAAGUAGGUCCCAAGGGUAUGGCUGUUCGCCA
UUUAAAGUGGUACGCGAGCUGGGUUUAGAACGUCGUGAGACAGUUCGGUCCCUAUCUGCCGUGGGCGCUGGAGAACUGAG
GGGGGCUGCUCCUAGUACGAGAGGACCGGAGUGGACGCAUCACUGGUGUUCGGGUUGUCAUGCCAAUGGCACUGCCCGGU
AGCUAAAUGCGGAAGAGAUAAGUGCUGAAAGCAUCUAAGCACGAAACUUGCCCCGAGAUGAGUUCUCCCUGACCCUUUAA
GGGUCCUGAAGGAACGUUGAAGACGACGACGUUGAUAGGCCGGGUGUGUAAGCGCAGCGAUGCGUUGAGCUAACCGGUAC
UAAUGAACCGUGAGGCUUAACCU

```

# Answer
You can check specification of pdbx\_seq\_one\_letter\_code and pdbx\_seq\_one\_letter\_code\_can. These items are generated from the depositor-provided sequence on the PDB side. I suppose that the software used is MAXIT.

Note: these are not PDB files. They are distributed by the PDB organization, but the file format is mmCIF (or, more formally, PDBx/mmCIF).
> 2 votes

# Answer
So, i was little myopic when reading the PDB description the first few times, but roughly what i was looking for is this:

1. the non-canonical sequence(`pdbx_seq_one_letter_code`) makes explicit the ligands present in the given chain.
2. the full dictionary of these is called the CCD and contains ~20000 entries, is curated and updated weekly by the PDB. That's basically their database (glorified text file over ftp) against which they check whether a given ligand has already been recorded.
3. the process of chaffing them out of the EM dough as well as schema and more context canbe found here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4243272/pdf/bau116.pdf
> 2 votes

---
Tags: sequence-alignment, pdb, file-formats, structural-biology
---


title: somatic SNV tool for ONT samples

text:
# Question
Title: somatic SNV tool for ONT samples
I'm looking for a tool to call somatic single nucleotide variations on reads from an Oxford Nanopore Technologies (ONT) run. I have both tumor and normal samples and have already aligned them to a reference genome.

Popular tools like mutect2 have been developed for short reads, and according to this paper it is not recommended to use tools specifically designed for short reads with long reads.

I know of tools that work with ONT reads like clair3 which is used in the epi2melabs workflow for SNPs wf-human-snp, but clair3 takes only one sample at a time, so seems for germline calling only.

Another approach would be to simply call the SNVs in both tumor and normal samples independently with clair3, and then subtract normal SNVs from tumor SNVs, but according to what I have read in a few papers and online resources like this from gatk, this is a "naive approach", and a somatic caller should use both normal and tumor samples together.

* Do you know of any tool that takes both tumor and normal long read samples to detect somatic variants (similar to mutect2 but for long reads)
* Have you used tools for short reads like mutect2 with long reads?
* Have you attempted the "naive approach" of simply subtracting SNVs?

Thanks, Walter.

# Answer
ClairS is now released for long-read somatic variant calling. https://github.com/HKU-BAL/Clairs
> 1 votes

---
Tags: variant-calling, nanopore, long-reads, somatic, snv
---


# Question
Title: Is the potential skewing of FDR and ranking of gene set analysis with gene sets with hierarchical structure a concern?
When performing GSEA and many other gene set analysis, when gene sets with hierarchy such as KEGG collection or Reactome collection, the FDR could be skewed sometimes, and the ranking of gene sets by any metrics could be confusing because multiple subsets belonging to the same set would have similar ranking, pushing other gene sets down the list. For example, in reactome axon guidance (R-HSA-422475.6) contains semaphorin interaction (R-HSA-373755.1), which contains Sema4d in semaphorin signaling (R-HSA-400685.2). The enrichment of Sema4d gene set usually (though not necessarily) entails the enrichment of the other gene sets above it. That might lead to an over-adjustment of p-value if I am not mistaken.

I am just wondering if there are ways to mitigate it, or it is not a major concern that I should have?

# Answer
This is actually a known problem in the correction of nominal p-values in over-representation analyses that use the Gene Ontology categories (commonly but misleadingly called "GO enrichment" in many publications).

Most currently used methods for p-value adjustment (like BH) assume that p-values are independent, which is not true for a case like the one you describe - a set contained by another set is basically by definition not independent. Even if you didn't have a hierarchical structure, you would still have many overlaps across all your genesets, which would still mean that you cannot assume independence. So in both cases, but I think especially in your case, a BH/FWER correction would be statistically wrong (and possibly overly conservative).

There are a few ways to avoid this: one way would be to prune the directed acyclic graph of your genesets/categories, retaining only the categories that are enriched and have no enriched children. In your case, for instance, if "axon guidance" and "semaphorin interaction" are significant but "Sema4d in semaphorin signaling" is not, you would stop at "semaphorin signaling" removing the result for "axon guidance".

This is kind of a strong correction, which is guaranteed to lower the multiple testing correction burden (unless you have no significant genesets).

There are however more refined ways to address this:

One of them is to remove from set A (parent) any genes that belong to a significantly over-represented/enriched set A' (child). This method is termed the **elim method** and was described in this publication.

The same publication discusses another method called **weight** which assigns a weight to a significant category based on the significance of its children.

Another method is implemented in the gProfiler tool and introduced in this 2007 publication, which is the **g:SCS method**. Namely it calculates p-value thresholds based on empirical, randomized p-values for sets of different sizes. According to the authors it is more conservative than BH but less than Bonferroni.

The only issue with these approaches is that there are, to my knowledge, no ways to apply them to your data independently, i.e. they are coded within the respective tools and do not exist as stand-alone packages for multiple test correction. I hope I'm wrong and someone developed these versions, but for the time being the only thing you can do is to try to implement your own versions of these methods.

Another relevant, more recent publication is the SetRank one, which describes a method to do GSEA with a p-value calculation that accounts for overlaps, and has its own package on CRAN. If you prefer to apply a correction to your own results (like the ones you get from `fgsea`), however, I do not know whether this package has methods you can use or recycle.
> 2 votes

# Answer
Statisticians might kill me for it but I've always filtered long and redundant geneset lists to only contain what remotely makes biological sense in my setup. If I am investigating leukocyte transcriptomes I do not care about terms related to filtration processes in the kidney or signal transduction in the hypothalamus. I think the redundancy of common pathway databases in terms of genes contributing to many terms justifiers to eliminate certain terms and choose those that make sense and are not overly large (for example not more than 500 genes). Top-of-hierarchy terms like 'metabolism' eith more than 1000 genes are anyway unspecific and hard to interpret. That then reduces the multiple testing burden quite a lot.
> 2 votes

---
Tags: pathway, gsea
---


# Question
Title: How to Recalculate Allele Frequency from a VCF File?
I have a VCF file with 200 samples (mitochondrial genome of Plasmodium falciparum). Here is a pic to take a look at:

And a few relevant lines from the actual file:

```
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  FP0008-C   FP0009-C
Pf_M76611       19      .       A       T       33196.8 MissingVQSLOD;Mitochondrion     AC=9;AF=0.0003019;AN=29816;ANN=T|intragenic_variant|MODIFIER|PF3D7_MIT00100|PF3D7_MIT00100|gene_variant|PF3D7_MIT00100|||n.19A>T||||||;AS_InbreedingCoeff=0.0437;AS_QD=20.96;BaseQRankSum=0.141;CDS;DP=559138;ExcessHet=0.0026;FS=2.892;InbreedingCoeff=0.52;MLEAC=8;MLEAF=0.0002683;MQ=60;MQRankSum=0;QD=20.96;ReadPosRankSum=0.856;RegionType=Mitochondrion;SOR=0.789;set=genotypegvcfs  GT:AD:DP:GQ:PL  0/0:52,0:52:99:0,120,1800  0/0:54,0:54:99:0,120,1800

```

A small subset of samples I am working on can be found here: https://www.mediafire.com/file/6dkct6guk5zx83m/sample.vcf/file

I am trying to calculate allele frequency (AF) for each variant in the dataset. The population is assumed to be all the samples in the file.

In other words, even though the VCF file contains AF statistics from population data, I would like to know how to recalculate AF only from the VCF file. I believe I need to know how AC and AN are calculated first.

AF is calculated this way: allele count (AC) / allele number (AN)

I tried to calculate them from the information from GT (counting 0/0, 1/1 and 0/1 in all the 200 samples), but then I learned that this calculation is not suitable for mithocondrial data.

All lines in the VCF file contain an AF statistic:

```
$ grep -v 'AF=' sample.vcf | grep -v '^#'
<No output>

```

The information is provided under INFO column, but how were they calculated? If I do not have AC and AN to begin with, what should I do? Is there a way to calculate AC and AN from the other information provided in the dataset? such as the data mentioned in the FORMAT column? GT: genotype AD: unfiltered allele depth DP : read depth at this position for this sample (Integer) GQ : conditional genotype quality, encoded as a phred quality PL : the phred-scaled genotype likelihoods rounded to the closest integer

Thanks

# Answer
> 2 votes
It sounds like your main confusion is how to calculate AF. I'll answer that right away, but I'll also suggest a Python library that you can use to do this and much more.

You're correct, the AF is defined as the number of non-missing alternate alleles (AC) divided the number of non-missing alleles (AN). Calculating these are straightforward. Map each GT to its alternate alleles:

```
0/0 => 0
0/1 => 1
1/1 => 2
./0 => 0
./1 => 1
./. => 0

```

Now, compute the sum. That is the AC. For the AN, map each GT to its number of non-missing alleles:

```
0/0 => 2
0/1 => 2
1/1 => 2
./0 => 1
./1 => 1
./. => 0

```

Now, compute the sum. That is the AN. Now calculate AF as

```
AC / AN

```

All of this is kind of annoying to do though! I suggest you use Hail (disclaimer: I am a Hail developer).

```
import hail as hl

# Hail is not yet aware of Plasmodium falciparum reference genome
mt = hl.import_vcf('sample.vcf', reference_genome=None)

# let's peek at our data
mt.show(n_rows=10, n_cols=5)

# Let's get some information about the GT specifically:
mt.GT.summarize()

# We can use Hail's built in variant_qc method to calculate AF for us:
mt = hl.variant_qc(mt)

# Let's take a look at that
mt.variant_qc.show()

# We can overwrite the old AF field and export back to a VCF file if necessary
mt = mt.annotate_rows(
    # AF[1] is the frequency of the first alternate allele
    # Your dataset appears to have some multi-allelics. The other
    # alternate allele frequencies are available in AF[2], ...
    info = mt.info.annotate(
        AF = mt.variant_qc.AF[1]
    )
)
hl.export_vcf(mt, 'sample-updated.vcf.bgz')

```

---

This comment doesn't make sense to me:

> I learned that this calculation is not suitable for mithocondrial data.

The definition of AF, AC, and AN doesn't change depending on the type of config (autosome, sex chromosome, mitochondrial). Perhaps the issue is ploidy encoding? A GT in a VCF can be:

* 1-ploidy, e.g. `0` (reference), `1` (first alternate), `2` (second alternate)
* 2-ploidy, e.g. `0/0` (homozygous reference), `0/1` (heterozygous), `1/1` (homozygous first alternate)
* etc.

Sometimes people incorrectly encode a 1-ploidy config (for example the Y chromosome) as: `0/0` and `1/1` instead of `0` and `1`. That's not your fault though and it doesn't mean the definition of AF, AC, and AN change. It means someone gave you malformed data.

---
Tags: vcf, allele-frequency
---


# Question
Title: Is there a software that extract molecular descriptors of single stranded RNAs secondary structure from dot bracket notation in nucleotide resolution?
I want to generate descriptors from dot bracket notations of single stranded RNA/DNA secondary structures. So far i have come across this paper:https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0097696#s5 In that, there is given a list of descriptors derived from the secondary structure of single stranded RNA but it is not explained how they derived these descriptors. Is there a Software out there, that would accomplish that task? I wrote a script myself, to extract loop counts etc. from dot bracket notation but it would rather tedious to write a script all by myself for all descriptors possible.

I would appreciate every constructive thought on that.

Thanks!

**EDIT** Summary of paper for relevance

In this paper the secondary structures of single stranded RNA is parsed into many descriptors, whereby the secondary structure is given as dot bracket notation and where obtained by mfold. THis paper is nothing more than an example of the variety of descriptors that are describing the secondary structure of single stranded RNA/DNA present in dot bracket notation. My question now is: Is there a general accepted software, that takes as input the dot bracket notation of an single stranded RND/DNA Oligonucleotide and returns many different descriptors? Such alike rdkit would do it for small organic molecules when given the SMILES string as input. I hope this clarified my need.

# Answer
> 0 votes
The authors use

> Kapetanovic IM (2008) Computer-aided drug discovery and development (CADDD): in silico-chemico-biological approach. Chem Biol Interact 171: 165–176.

This generates the descriptors, which are normally used in small molecule interactions and they obtained a moderate (in reality low for ML) correlation. That is the innovative part of the study. "Decision trees" used alone are not recommended because it will over tighten onto the data. The reason they used them was to maximise the fit - the risk is over tightening. When its says "the algorithm successfully ...", well thats using a Decision Tree, any other ML algorithm will be lower.

The very nice part of the paper is the combination of wet-lab and analytical work, that balance is often not logistically trivial to achieve.

What is not clear is they talk about 6 descriptors out of 60 and also 27 descriptors out of 60. Which is the one they use? I would suggest the paper was not rigorously reviewed/edited to permit this sort of confusion for the reader. It could be (suspect) they performed two different analyses, course grain (6) and fine grain (27), but I don't think that is described very well.

I will not critique further.

**Summary** The paper is commendable for using a small molecular physiochemical/ thermodynamic description algorithm onto RNA data. Thats cool. Wrapping it around a Decision Tree alone, I wouldn't agree with - nothing wrong with DT's but they can't be used in isolation.

> QSAR Model

> In order to find the mathematical relationship between aptamers’ structure and activity in the training set, two functions of MASS packages of R program were used (R 2.14.1 2011 The R Foundation). Firstly, a stepwise-selection with the function “stepAIC” (direction = both) focused on the most relevant 27 descriptors out of the 60 that can contribute to the prediction of the relative activity. Then, a second step with multiple-linear-regression function “lm” was applied on the selected descriptors to generate appropriate coefficients for the linear correlation equation. The quality of the output was evaluated according to the p-value of the null hypothesis, and multiple R2 indicated the correlation between actual and predicted values.

> QSAR

> Sixty descriptors quantifying structural and thermodynamic properties of the aptamers were produced and gradually the most relevant descriptors were selected to form a six descriptor- based algorithm, which can predict the binding of aptamers to influenza virus and thus anticipate its anti-viral activity. The algorithm successfully predicted the binding affinity in the training set of aptamers (R2 = 0.702) and the was validated with a test set (R2 = 0.66); moreover, the model’s sensitivity and specificity in selecting aptamers with enhanced relative binding were high, 89% and 87%, respectively.

---

Personally, I don't think that is the question the OP is asking - I think they are asking about a different technique in RNA secondary structure analysis to the technique authors performed here. That however is simply a guess. It's an interesting paper, there is no question about that - it's not how I personally would do it - but it's interesting. The issue is the thermodynamics has been rigorously worked out historically for RNA folding - it's a large body of work - and I'm sure the one or two reviewers would have pointed this out. It's complicated and if the authors are using their predictions to recycle back to the wet-lab .... thats what it's all about.

---

Normally, folding is calculated using graph theory ... FYI.

---

From there comments .... whats the program ... its Kapetanovic (2008) as cited above, which will be a classic small molecular binding algorithm the authors called `QSAR`.

In cases there's any misunderstanding - we wouldn't use small molecule interaction based calculations (MD) to understand RNA folding/binding because there's an alternative body of theory for this specific calculation, partly thermodynamics, partly a separate factor. ML is heavily used but not via an MD style calculation - which is what the authors are doing. Thus if you use the method there is a risk of criticism.

It's complicated - I understand the authors perspective, they are saying classical style hybridisation its just not that simple. It might be complicated but hybridisation shouldn't be ignored as the empirical assumption.

---
Tags: rna, secondary-structure
---


title: blood contamination rna seq

text:
# Question
Title: blood contamination rna seq
I have tumor rna seq samples. I want to find out percentage of blood contamination within these samples. What is the best possible way to find % of blood contamination within the rna seq samples?

# Answer
> 5 votes
I assume by "blood contamination" you are interested in tumor purity and immune cell infiltration. For that you can use a deconvolution tool for bulk RNA-seq like CIBERSORT or quanTIseq

---
Tags: rna-seq
---


# Question
Title: Differential gene expression: how do I correctly compare three groups to each other?
I'm trying to figure out the right way to do differential gene expression analysis with a discrete variable with three groups. The context is, I want to assess differential expression as a function of whether samples are WT, have heterozygous loss, or homozygous loss for a particular gene. A key part of it is comparing the differences between the WT and het-loss samples with the differences between het-loss and hom-loss samples.

I tried doing just doing differential expression with DESeq2, treating it as a categorical variable, but I got a strange result: differential expression between het and hom samples seemed to usually go in the opposite direction as in wt vs. het (so if a gene was upregulated in het relative to WT, it was downregulated in hom samples relative to het). This was very counterintuitive; why would losing the second copy of a gene consistently have the opposite effect on overall gene expression that losing the first copy has? So to test if something was going awry, I did the same analysis, but randomly assigning samples into three groups, and found the same thing: the log fold change comparing group 3 to 2 was strongly negatively correlated with the log FC comparing group 2 to 1, even though the groups are entirely random. So this must be some artifact of the model being used (note that I observe the same thing when used limma voom instead of DESeq 2). I include an example here from data loaded from TCGA via TCGAbiolinks:

```
library(TCGAbiolinks)
library(DESeq2)

#Loading KIRC data from TCGA:
query_TCGA = GDCquery(
  project = "TCGA-KICH",
  data.type = "Gene Expression Quantification",
  data.category = "Transcriptome Profiling", # parameter enforced by GDCquery
  experimental.strategy = "RNA-Seq",
  workflow.type = "STAR - Counts")
download <- GDCdownload(query = query_TCGA)
dat <- GDCprepare(query = query_TCGA, save = FALSE)
rna <- as.data.frame(SummarizedExperiment::assay(dat))

#Randomly assigning samples to groups 1, 2, and 3 and running DESeq2:
set.seed(1234)
xvar <- factor(sample(1:3, ncol(rna), replace = TRUE), levels = c(1,2,3))
columnData <- data.frame('Patient'=colnames(rna),'xvar'=xvar)
columnData$xvar <- factor(columnData$xvar, levels = c('1','2','3'))
deseq2Data <- DESeqDataSetFromMatrix(countData=rna, colData=columnData, design= ~     xvar)
deseq2Data <- deseq2Data[rowSums(counts(deseq2Data)) > 5, ]
deseq2Data <- DESeq(deseq2Data)
res.1_v_3 <- as.data.frame(results(deseq2Data, contrast = c('xvar','1','3')))
res.1_v_2 <- as.data.frame(results(deseq2Data, contrast = c('xvar','1','2')))
res.2_v_3 <- as.data.frame(results(deseq2Data, contrast = c('xvar','2','3')))

#Merging comparisons of 1 v 2, 1 v 3, and 2 v 3 into one data frame:
merged <- merge(data.frame('Gene'=rownames(res.1_v_2),     'log2FC_1_v_2'=res.1_v_2$log2FoldChange, 'padj_1_v_2'=res.1_v_2$padj),
                data.frame('Gene'=rownames(res.2_v_3),     'log2FC_2_v_3'=res.2_v_3$log2FoldChange, 'padj_2_v_3'=res.2_v_3$padj), by='Gene')
merged <- merge(merged, data.frame('Gene'=rownames(res.1_v_3),     'log2FC_1_v_3'=res.1_v_3$log2FoldChange,
                                   'padj_1_v_3'=res.1_v_3$padj), by='Gene')
#Color coding for illustration:
merged$color <- 'gray'
merged$sig_1_v_2 <- merged$padj_1_v_2 < .05
merged$sig_1_v_3 <- merged$padj_1_v_3 < .05
merged$sig_2_v_3 <- merged$padj_2_v_3 < .05
merged$color[which(merged$padj_1_v_2 < .05 & merged$padj_2_v_3 > .05)] <- 'blue'
merged$color[which(merged$padj_1_v_2 > .05 & merged$padj_2_v_3 < .05)] <- 'red'
merged$color[which(merged$padj_1_v_2 < .05 & merged$padj_2_v_3 < .05)] <- 'purple'
maxim <- max(abs(c(merged$log2FC_1_v_2, merged$log2FC_2_v_3)))
coef <- cor(merged$log2FC_1_v_2, merged$log2FC_2_v_3)
p <- try(round(chisq.test(table(merged[,c('sig_1_v_2','sig_2_v_3')]))$p.value,     digits = 4), silent = TRUE)
if(class(p) == 'try-error'){p <- NA}

plot(merged$log2FC_1_v_2, merged$log2FC_2_v_3, col=merged$color, pch=16,     cex=.5,main='Log2 FC',
     xlim=c(-maxim1,maxim), ylim=c(-maxim,maxim),xlab=c('LogFC:     Het/WT'),ylab=c('LogFC: Biallelic/Het'))
abline(a=0,b=coef1,col='red')
legend('topright',legend = c('Coefficient','Significant Group 2 v. Group 1',
                             'Significant Group 3 vs. group 2','Significant for     both comparisons'),
       lty=c(1, NA, NA,NA), pch=c(NA, 16, 16,16),     col=c('red','blue','red','purple'))
text(x=-maxim*.8, y=-maxim, labels=paste('p=',p,' (chi-sq. test',sep=''))

```

Is there something particular I'm doing wrong? And what would be the correct way to do differential expression analysis with an independent variable that is essentially an ordinal variable? Thanks.

# Answer
I will try to answer this open-ended question with general statements towards the analysis of larger cohorts:

**Larger cohorts are usually not homogeneous.**

Based on the provided annotations (`colData(dat)`) there are cancer samples **and** normal controls in this dataset, and they're not equal in size. It is 66 `Primary solid Tumor` and 25 `Solid Tissue Normal`. Running a PCA, for example via `plotPCA()` after running `vst()` will demonstrate that this is the main source of variation in this dataset, separating along PC1. That means a random split of the dataset will still have more or fewer cancers or normals in one or the other groups, and at that sample size and based on the expected strong expression differences between cancer and normal this will get DEGs almost certainly.

Here the PCA via `plotPCA(vst(dds, blind=TRUE), "definition")`

That having said, you need to inspect the data you deal with in depth before throwing any test on it. If not you might test some sort of confounders or factors that bias the analysis without noticing it. Based on the PCA there is also **notable heterogeneity** in the cancer samples, and the annotations list additional factors that can contribute to this, be it **gender**, **smoking status**, **tumor stage**, **age** and others. You would need to ensure that your random assignment does not nest any of these factors or the tumor/normal type with the random groups.

Data exploration is key. There was a recent paper that got quite some attention on Twitter claiming that common testing frameworks (incl DESeq2) lead to exaggerated false positive rates. However, the DESeq2 author (analysis here) showed that this paper missed aforementioned proper data exploration, and the dataset infact had unaddressed technical variation and samples forming groups within that dataset. Addressing that (he did it via `RUVSeq`) then deemphasized the strong claims in that paper.

Summarizing: Large cohorts, especially human ones, often have internal nesting and confounding with other factors that are unwanted. That can be the aforementioned obvious ones like age or gender (sex) but also more subtle ones often not recorded, such as dietary status, unreported drug consumption, amount of physical exercise the person does (if any), genetic predisposition, previous disease or accidents/trauma and so on. That having said, I am quite sure that the "random" DEGs you see are nested unwanted variation rather than anything "wrong" with the testing frameworks.
> 3 votes

---
Tags: r, deseq2, differential-expression
---