# Question
Title: Error while loading a big file with Cytoscape
I installed Cytoscape 3.7.1 with java 1.8.0\_191 in windows server 2012. I have a 7.4 GB `csv` file (about 1,500,000 reords) and when I tried to load it into the Cytoscape it throws an error `java.lang.OutOfMemoryErroe:Java heap space.`

My system has 500GB ram and then I changed `Cytoscape.vmoptions` file according to `-Xms102400M` but when I load my file, in task manager, the memory increased until 50GB and then Cytoscape crashed and again it threw the same error: `java.lang.OutOfMemoryErroe:Java heap space.`

My data is like this sample:

```
45677888,test1,3453453453,test2,3534523453,235412352,3453452345,test3,235423452,test4,1980/09/02,14:13,23523525,test5
45234288,test11,34234553453,test12,353434534553,23453452352,342422345,test13,23456543452,test14,1980/04/01,14:12,2323423425,test15
4243424888,test12,3235253453,test22,3533456343453,233534532352,2343452345,test33,23345342,test44,1980/03/4,14:11,23674575,test55

```

My file is in CSV format. col2 is source node, col1 and col3 are attributes for source node,col 4 is destination node, col5,col6,col7 are attributes for destination node, col8 until end are attributes for edge. I have about 200,000 nodes and 230000 edges.

How can I read this file in Cytoscape?

# Answer
> 1 votes
The first thing you can try is to separate the network from the meta information. Assuming you have access to a \*nix machine, you can run:

```
awk -F"," '{print $2,$4}' file.csv > network.csv

```

That will produce a file like this (based on your example data):

```
test1 test2
test11 test12
test12 test22

```

That should radically decrease the file size and will give you the best chance of loading it in Cytoscape. You can then similarly extract the attributes:

```
awk -F"," '{print $2,$1}' file.csv > source.column1.attributes
awk -F"," '{print $2,$3}' file.csv > source.column3.attributes
awk -F"," '{print $4,$5}' file.csv > target.column5.attributes
awk -F"," '{print $4,$6}' file.csv > target.column6.attributes
awk -F"," '{print $4,$7}' file.csv > target.column7.attributes
awk -F',' '{print $2,$4,$8}' file.csv > edge.column8.attributes
awk -F',' '{print $2,$4,$9}' file.csv > edge.column9.attributes
awk -F',' '{print $2,$4,$10}' file.csv > edge.column10.attributes
awk -F',' '{print $2,$4,$11}' file.csv > edge.column11.attributes
awk -F',' '{print $2,$4,$12}' file.csv > edge.column12.attributes
awk -F',' '{print $2,$4,$13}' file.csv > edge.column13.attributes
awk -F',' '{print $2,$4,$14}' file.csv > edge.column14.attributes

```

Now, try loading `network.csv` into Cytoscape and then load each of the node and edge attributes separately.

---

However, a network with 200,000 nodes and 230,000 edges is just too large to be able to manipulate in any useful way in a GUI. It will just look like a horrible hairball. I really urge you to do some preprocessing and pare the network down to something more manageable. I used to work with a human PPI network of ~13000 nodes and ~15000 edges. While I *could* load that into Cytoscape, that was already far too big to be able to do anything useful with. I could use the data analysis aspects of Cytoscape, but not really the visualization.

# Answer
> 1 votes
I work a lot with large networks and I agree with what most other posters have said, in that it is hardly useful to explore a network like this visually and Cytoscape - while it is one of best GUIs in general - is pretty bad with large graphs. Here are some alternatives you could try *if you want to go down that road*:

* Gephi also probably won't handle that number of nodes/edges well, but it may not crash. Development on it has stagnated though.
* Graphia I haven't tried this, but it is supposed to be better at handling large graphs. This project is in its infancy, so don't expect too much

If you care about the analysis of the graph itself, I'd go with NetworKit or SNAP which are much better suited to large graphs.

If you want to visualize the networks, something that tends to perform well here is:

1. Graph embeddings with VERSE
2. Dimensionality reduction of the embeddings with UMAP
3. Plot with Datashader

This can handle absolutely massive graphs if you have the compute - which you seem to have. Here's one I am just working on (albeit small compared to yours with 30K nodes and 800K edges) - from that exact pipeline. Took ~10 minutes read to render on my desktop:

Finally, as others have mentioned, graph partitioning might be the way to go. I really like InfoMap. It partitions the graph into a hierarchical structure that you can explore top down.

# Answer
> 0 votes
That's certainly a lot of nodes and edges, however, there is a preference called viewThreshold that should prevent the creation of a view if the number of nodes+edges is larger than the set value. I've just tested it on a large network I have (9000 nodes and 4,000,000 edges) and it doesn't appear to be working properly. I'll mark this as a bug. It does, however, render, but it's not really a useful network on Cytoscape 3.7. It will actually be usable on Cytoscape 3.8 based on early testing, however. If you want to make your network available somewhere, I'm happy to try to see why it's not loading.

-- scooter

---
Tags: networks
---


title: Discard sequences contain Ns

text:
# Question
Title: Discard sequences contain Ns
I have Fasta file contain multiple sequences, some of those sequences contain N's or n's, I'm trying to discard all the sequences and related header that contain Ns or ns in the file.

I have tried multiple python code but return an error.

```
import sys
from Bio import SeqIO
handle = open(sys.argv[1], "rU")
filtered = [record for record in SeqIO.parse(handle, "fasta") if record.seq.count('N') == 0]
output_handle = open("N_removed.fasta", "w")
SeqIO.write(filtered, output_handle, "fasta")
output_handle.close()
handle.close()

```

Error:

```
NsScript.py:3: DeprecationWarning: 'U' mode is deprecated
  handle = open(sys.argv[1], "rU")

```

and

```
import sys
from Bio import SeqIO
for record in SeqIO.parse(sys.argv[1], "fasta"):
    if record.seq.count('N') == 0:
        print(record.format("fasta")

```

Error:

```
/RemoveNs.py", line 6

                                        ^
SyntaxError: unexpected EOF while parsing

```

# Answer
The 1st code snippet is not an error, it's a deprecation warning. You can remove that warning by dropping the `"U"` and just specifying `"r"` as the mode. (The `"U"` is for universal newlines, which is now the default behaviour since python 3, and should be specified by a separate `newline` arg anyway).

The 2nd code snippet is a trivial error because you're missing a closing parantheses `)` on the final line.
> 2 votes

# Answer
If you're not married to Python, here's a one liner in bioawk:

```
bioawk -c fastx '{ if($seq !~ /[Nn]/) {print ">"$name; print $seq}}' old_file.fa > new_file.fa

```
> 2 votes

# Answer
With Visual Studio Code you can delete them using a regular expression. Replace `^>.*\n[^>]*N[^>]*` with an empty string You can also do that in a similar way using the sed Linux command. You can use regular expressions in Python too.
> 1 votes

---
Tags: filtering
---


# Question
Title: Download complete Phylogenetic Tree as JSON
I would like to integrate the phylogenetic tree into a website. Does anyone know where I can download the complete tree? Best would be a JSON, which looks something like this:

```
[{
  "level": "domain",
  "name": "eukaryota",
  "commonNames": {
    "en": ["eukaryotes"],
    "de": ["Eukaryoten", "Eukaryonten"]
  },
  "children": [
    {
      "level": "kingdom",
      "name": "plantae",
      "commonNames": {
        "en": ["plants"],
        "de": ["Pflanzen"]
      },
      "children": [...]
    },
    {
      "level": "kingdom",
      "name": "animalia",
      "commonNames": {
        "en": ["animals"],
        "de": ["Tiere"]
      },
      "children": [...]
    },
    ...
  ]
},
{
  "level": "domain",
  "name": "bacteria",
  "commonNames": {
    "en": ["bacteria"],
    "de": ["Bakterien"]
  },
  "children": [...]
}]

```

Other formats and structures are ok, too. The important thing is that the level of the nodes can be retrieved. In the end I only want to have the main levels Domain, Kingdom, Phylum, Class, Order, Family, Genus and Species in the tree. (https://en.wikipedia.org/wiki/Taxonomy\_(biology)#/media/File:Biological\_classification\_L\_Pengo\_vflip.svg) So if there are more nodes from other levels like Clade, Tribe or Subfamily, it must be possible to identify and remove them. Thanks a lot!

# Answer
NCBI maintains an SQL dump of their taxonomic database

You can see what that contains here

For smaller requests (depending on what you need), you might also consider the ENSEMBL REST API - although I don't know of a way to retrieve the complete tree (shoot them an email): https://rest.ensembl.org/taxonomy/classification/53399?content-type=application/json
> 4 votes

# Answer
It would be worth looking at the Tree of Life web project, it isn't downloadable but its a reference source that can be improved on and updated in your website.

ToL contains trees linking everything from bacteria, archaea through to the complete diversity of eukaryotes, but its a bit old with 2008 from my taxonomic knowledge being the last update for some of the entries, e.g. Plasmodium spp..
> 1 votes

---
Tags: phylogenetics, phylogeny, taxonomy
---


title: Apply trained PRS on another dataset

text:
# Question
Title: Apply trained PRS on another dataset
I am using PRSice to compute the PRS over a train set and want to use the coefficient used on the train set to apply it on another set which I will call the test set.

Once I compute the PRS I get a file called analysis.prsice

This is the few first lines

---

```
Pheno Set Threshold R2 P Coefficient Standard.Error Num_SNP

- Base 5e-08 0.00550157 0.0659631 -21.0513 11.4492 43

- Base 5.005e-05 0.00396637 0.113759 -61.8545 39.1108 184

- Base 0.00010005 0.00236122 0.220463 -68.9527 56.2742 286

- Base 0.00015005 0.00113035 0.394503 -59.4058 69.7677 382

- Base 0.00020005 0.00072253 0.495129 -55.2605 81.0064 471

- Base 0.00025005 0.000334179 0.642121 -41.5674 89.4432 538

- Base 0.00030005 7.5306e-05 0.825111 -21.391 96.8021 601

- Base 0.00035005 3.95773e-06 0.959527 5.25804 103.612 668

- Base 0.00040005 0.000175506 0.735018 37.723 111.455 743

```

---

I thought that I could apply the PRS on another dataset by multiplying the coefficient in the result by the allelic dosage of the corresponding SNPs. However, I am not sure how to know which are the corresponding SNPs of the coefficients.

Also, if there was a command that does that and it would be easier, could you please tell me how to do it?

Thanks

# Answer
> 0 votes
The coefficient is the third column from the right, you can tell because some are positive and other negative. When the value is negative the association is a negative correlation. When it is positive (last two) that is more like the result you are seeking. The further away from zero the stronger the correlation (positive or negative)

These are very large values for the coefficients. The SNP is the last column on the right. I would disgard the top two results because the probabilities are too borderline, e.g. if you used Bonferoni's correction.

There isn't enough information on this calculation nor the data set to make any further insights other than you have 7 SNPs showing PRS

---
Tags: snp, gwas, phylogenetics
---


# Question
Title: Can I eject the flow cell MinION from USB port when doing basecalling?
Can I eject the flow cell MinION from USB port when doing basecalling? Because its taking long and its not possible to wait for its finished, will it affects basecalling?

# Answer
You can, you will lose the transfer from the MinION. It would probably be better to just stop the run. The computer will continue to basecall untill it is done. All the raw reads are stored locally on the machine so it should be fine.
> 1 votes

# Answer
Yes, it should be fine, especially if the run is finished. MinKNOW should still continue basecalling even after the MinION is disconnected. Unfortunately I can't find good posts in the ONT community about that because search is currently broken. This post might be useful (found by searching my email records), but it's two years old:

> Hi, thank you for your post. You can certainly remove the flow cell and wash it without affecting the catch-up mode. Please see the downloads page here for the most up-to-date installer version.
> \- Monolina Binny, ONT

If you're likely to be doing more long runs in the future (\>6h), I'd recommend getting a cheap video card with CUDA capability (at least v10) to save time in the future.
> 0 votes

---
Tags: fastq
---


# Question
Title: Comparing AUCs: Discrimination of same Control from different Test Group - paired or unpaired? miRNA
I want to compare two AUCs using a bootstraping method from pROC package (roc.test).

I have done a classification using Logistic Regression and Serum miRNAs as prediction factors. The Controls are in both cases the same, the test group is different.

```
1. Cancer (high risk) vs. Control = AUC 1
2. Cancer (low risk) vs. Control = AUC 2

```

If I want to compare this two AUCs - should I use a paired or unpaired test?

I would tend to using a paired test, because of exact the same Control group (Control groups are correlated).

# Answer
Generally its paired. The bootstrapping personally I find concerning. I presume you are using it to augment your data set.

Bootstrapping at best tends to be a controversial statistical operation and could result in the non-significance you are observing. Anyway if its not significant between high and low - thats your answer.

You could of course perform a single ML classifying low, and high cancer risk against control. That would provide a method to understand the associations of miRNA.

---

Bootstrapping issue, I had a think about this. It depends on how the bootstrapping is done. In a paired sample this (I assume) samples with replacement at each paired position for both ROCs. If this is not the case - the test is junk. However, if this is the case, its actually a good test, because the paired structure of the data is always preserved no matter how many bootstrap samples are performed. So like if a 100 were done and the pairwise statistic (which we don't know what it is - but something like KS test) and its \>0.05 for all 100 samples - I don't think there's any doubt thats a robust result.

Again if the two ROCs are bootstrapped independently, in my opinion the test is junk.

**Reason** The concern of bootstrapping is it causes weird distributions which were not present in the original data. This is only okay if there's lots of experience about interpreting the result. For example, in phylogenetics there's decades of experience interpreting bootstrapping and a consensus about what it means. This often isn't the case for most other areas of statistics.

I note this particular issue is really ramping up the rep score ;0)
> 0 votes

---
Tags: statistics, cancer, design, logistic-regression
---


# Question
Title: Genome assembly of SRR12196449 with SPAdes
I am trying to assemble the run SRR12196449 with `SPAdes`. The description of their project is:

> This project expected to standardize a method for amplification and sequencing of FIV genome in a simple way, allowing a broader analysis to increase the knowledge on the biology and evolution of the virus and virus-host interaction.

This is from School of Veterinary Medicine; University of Sao Paulo.

I have tried in three ways:

1. Using a fasta file downloaded from NCBI (using Filtered Download method)

> SPAdes-3.14.0/spades.py --only-assembler -s sra\_data.fa -o raw-fiv1

2. Using a fastq file also downloaded from NCBI using Filtered Download method.

> SPAdes-3.14.0/spades.py --12 sra\_data.fastq.gz -o raw-fiv1-b

3. Downloading the original file and dumping to a fastq using sratools.

> SPAdes-3.14.0/spades.py --12 SRR12196449.fastq -o raw-fiv1-c

In the third case I get a very bad assembly, likely because it needs trimming. In the two former cases, I get an assembly with a max contig length of around 4k-5k bps. But when I use `quast` to evaluate the assembly I obtain a largest alignment of ~500bps. I use this as reference.
Furthermore, if make a blast search it finds that the 4k-5k contig matches other FIV sequences at ~90% . In between these other sequences, there is at least one (MF370550.1) submitted by the University of Sao Paulo (likely another run from another sample in the same project).

On the other hand, if I assemble the run using another software I obtain a contig that is an almost perfect match to my reference sequence (~99%).

I have a few questions.

* Am I using `SPAdes` correctly? Is there a better way to use it?
* Is it possible that the sequences that blast finds matching with the `SPAdes` contigs are indeed sequences that someone got using `SPAdes` and that's why these match (not because these are correct)?
* What other option could explain these circumstances that I could be missing?

EDIT: Some clarifications regarding some comments from @MaximilianPress. I can confirm that I used the same reference genome in both cases. I used that command to get the quast results:

```
quast-5.0.2/quast.py -r raw-fiv1/sequence.fasta SPAdes-3.14.0/raw-fiv1/contigs.fasta

```

The other assembler is an overlap-layout-consensus algorithm that I implemented. The methodology is similar: I use exactly the same input fasta file.

I have also tried using the `SPAdes` --plasmid flag and that doesn't improve much the result. Max contig length is now 3642 bps and largest alignment 737 bps, covering only 26% of the reference. BTW, I made also tests with other virus sequences. I am just interested to know if I am using `SPAdes` correctly so that I can compare the performance of my algorithm Here there is a link to additional cases with other genomes. There you can find all the data and results used (raw files, assembly results and quast results)

Another clarification: I have also got perfect assemblies using `SPAdes` for other runs corresponding to other viruses... I am following exactly the same procedure in this case.

# Answer
> 4 votes
**Update 2:**

I looked into this a little more, with the various data sources.

This is related in part to the answer submitted by OP juanjo75es, in addition to discussion on chat. I don't entirely understand the logic, but the general thrust seems to be that SPAdes makes weird assemblies somehow.

Some notes that I made:

**REFERENCE ASSEMBLIES**

* FIV sequence U11820.1 was deposited in 1996, before SPAdes existed. Unclear assembly method.
* FIV sequence MN630242 was deposited in 2020, but was assembled with CLC workbench.
* These two sequences align together ok with minimap2. Low exact identity, it's true. This is the PAF output:

```
MN630242.1      8977    194     8964    +       U11820.1        9462    575     9326    1130    8778    60      tp:A:P  cm:i:116        s1:i:1125
       s2:i:55 dv:f:0.1761

```
* These two assemblies are syntenic, and both have similar gene calls using prokka. GFFs:

```
MN630242.1      Prodigal:2.6    CDS     253     1605    .       +       0       ID=AMGANPBD_00001;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00001;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     1656    4868    .       +       0       ID=AMGANPBD_00002;eC_number=3.6.1.23;Name=dut;gene=dut;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:Q2YRG4;locus_tag=AMGANPBD_00002;product=Deoxyuridine 5'-triphosphate nucleotidohydrolase
MN630242.1      Prodigal:2.6    CDS     4861    5616    .       +       0       ID=AMGANPBD_00003;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00003;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     5891    8461    .       +       0       ID=AMGANPBD_00004;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00004;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     8626    8790    .       +       0       ID=AMGANPBD_00005;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00005;product=hypothetical protein

U11820.1        Prodigal:2.6    CDS     634     1983    .       +       0       ID=EANIPDKN_00001;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00001;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     1995    5246    .       +       0       ID=EANIPDKN_00002;eC_number=3.6.1.23;Name=dut;gene=dut;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:Q2YRG4;locus_tag=EANIPDKN_00002;product=Deoxyuridine 5'-triphosphate nucleotidohydrolase
U11820.1        Prodigal:2.6    CDS     5239    5994    .       +       0       ID=EANIPDKN_00003;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00003;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     6269    8830    .       +       0       ID=EANIPDKN_00004;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00004;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     8904    9152    .       +       0       ID=EANIPDKN_00005;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00005;product=hypothetical protein

```

I can share FAA files of the proteins if needed.

**READ ASSEMBLY**

I also assembled the indicated reads using SPAdes. For reference it is a **~9Kbp** virus, but this is a **4.3Mbp** assembly. There is a lot of non-virus sequence in there. The second-largest contig is a shuttle vector. the third-largest contig is cat (host). Many more are cat, so I think it's pretty cat-oriented. The original authors used CLC workbench, so I guess that just worked a lot better in this instance, even in the presence of all the contamination. Unclear why, it appears that CLC works similarly to OP's assembly tool.

I aligned these to MN630242.1 with minimap2. If I understand OP, they are unhappy about the mapping of these contigs to this genome reference. The identities are mostly high (though there is indeed some oddness with the largest contig, which only finds very low coverage- at the same time that there are overlapping contigs with very high ID?), and they cover the entire reference genome:

```
MN630242.1  8977    2155    3199    -   NODE_14_length_1054_cov_3786.620280 1054    2   1046    1021    1044    60  tp:A:P  cm:i:189    s1:i:1021   s2:i:378    dv:f:0.0031
MN630242.1  8977    7893    8626    -   NODE_25_length_743_cov_3814.258117  743 8   739 713 733 60  tp:A:P  cm:i:129    s1:i:713    s2:i:318    dv:f:0.0045
MN630242.1  8977    72  4904    -   NODE_1_length_4942_cov_25.814123    4942    34  4863    628 4832    41  tp:A:P  cm:i:7s1:i:628  s2:i:512    dv:f:0.1703
MN630242.1  8977    3376    3897    +   NODE_39_length_526_cov_5254.155388  526 3   524 512 521 0   tp:A:S  cm:i:8s1:i:512  dv:f:0.0023
MN630242.1  8977    6270    6673    +   NODE_82_length_409_cov_1.744681 409 3   406 357 403 0   tp:A:P  cm:i:53 s1:i:357    s2:i:351    dv:f:0.0185
MN630242.1  8977    6852    7223    -   NODE_87_length_381_cov_23.334646    381 6   376 355 371 2   tp:A:P  cm:i:5s1:i:355  s2:i:351    dv:f:0.0079
MN630242.1  8977    6187    6558    -   NODE_88_length_380_cov_2338.675889  380 7   378 351 371 0   tp:A:S  cm:i:5s1:i:351  dv:f:0.0085
MN630242.1  8977    6852    7308    +   NODE_56_length_471_cov_0.933140 471 15  471 351 456 0   tp:A:S  cm:i:43 s1:i:351    dv:f:0.0397
MN630242.1  8977    7332    7707    +   NODE_86_length_389_cov_3506.202290  389 5   380 334 375 17  tp:A:P  cm:i:5s1:i:334  s2:i:309    dv:f:0.0111
MN630242.1  8977    5181    5568    -   NODE_60_length_464_cov_1.005935 464 50  437 314 387 17  tp:A:P  cm:i:41 s1:i:314    s2:i:282    dv:f:0.0357
MN630242.1  8977    7277    7707    +   NODE_76_length_438_cov_0.919614 438 5   435 309 430 0   tp:A:S  cm:i:40 s1:i:309    dv:f:0.0437
MN630242.1  8977    6953    7285    +   NODE_94_length_345_cov_2.903670 345 8   340 301 332 0   tp:A:S  cm:i:50 s1:i:301    dv:f:0.0110
MN630242.1  8977    5744    6056    -   NODE_113_length_316_cov_2.169312    316 4   316 289 312 35  tp:A:P  cm:i:4s1:i:289  s2:i:244    dv:f:0.0096
MN630242.1  8977    5615    5927    +   NODE_108_length_322_cov_1103.400000 322 6   318 283 312 58  tp:A:P  cm:i:4s1:i:283  s2:i:209    dv:f:0.0143
MN630242.1  8977    6543    6813    -   NODE_147_length_279_cov_2592.519737 279 8   278 261 270 21  tp:A:P  cm:i:4s1:i:261  s2:i:236    dv:f:0.0032
MN630242.1  8977    8544    8872    +   NODE_100_length_333_cov_2259.189320 333 1   329 252 328 0   tp:A:P  cm:i:3s1:i:252  s2:i:245    dv:f:0.0318
MN630242.1  8977    7601    7853    -   NODE_200_length_255_cov_25.125000   255 0   252 250 252 12  tp:A:P  cm:i:3s1:i:250  s2:i:236    dv:f:0.0017
MN630242.1  8977    6648    6968    -   NODE_106_length_324_cov_2.538071    324 4   324 248 320 0   tp:A:P  cm:i:3s1:i:248  s2:i:248    dv:f:0.0303
MN630242.1  8977    5940    6193    -   NODE_162_length_266_cov_1809.287770 266 6   259 246 253 5   tp:A:P  cm:i:4s1:i:246  s2:i:240    dv:f:0.0033
MN630242.1  8977    5416    5664    -   NODE_264_length_255_cov_9.531250    255 5   253 243 248 0   tp:A:P  cm:i:4s1:i:243  s2:i:242    dv:f:0.0033
MN630242.1  8977    7221    7466    -   NODE_247_length_255_cov_12.539062   255 1   246 241 245 0   tp:A:P  cm:i:4s1:i:241  s2:i:241    dv:f:0.0032
MN630242.1  8977    4830    5071    -   NODE_1121_length_245_cov_11.711864  245 4   245 239 241 0   tp:A:P  cm:i:4s1:i:239  s2:i:235    dv:f:0.0042
MN630242.1  8977    7769    8013    +   NODE_677_length_251_cov_3.951613    251 6   249 238 244 7   tp:A:P  cm:i:4s1:i:238  s2:i:230    dv:f:0.0075
MN630242.1  8977    4901    5148    +   NODE_270_length_255_cov_8.726562    255 8   255 238 247 0   tp:A:P  cm:i:4s1:i:238  s2:i:234    dv:f:0.0073
MN630242.1  8977    8705    8948    -   NODE_382_length_254_cov_5.763780    254 7   250 234 243 15  tp:A:P  cm:i:4s1:i:234  s2:i:217    dv:f:0.0062
MN630242.1  8977    6065    6307    +   NODE_410_length_254_cov_3.488189    254 1   243 227 242 10  tp:A:P  cm:i:3s1:i:227  s2:i:210    dv:f:0.0122
MN630242.1  8977    5062    5239    -   NODE_7003_length_180_cov_9.207547   180 3   180 172 177 17  tp:A:P  cm:i:2s1:i:172  s2:i:157    dv:f:0.0066
MN630242.1  8977    11  136 +   NODE_21321_length_128_cov_1195.000000   128 2   127 125 125 3   tp:A:P  cm:i:2s1:i:125  s2:i:123    dv:f:0

```

Only a (relatively) few contigs align well (figure), as might be expected with heavy contamination, and of those some are overlapping:

To look into the odd behavior of the biggest contig, I ran prokka on it too and found the expected genes for its position, especially the largest gene in the virus, `dut` (Deoxyuridine 5'-triphosphate nucleotidohydrolase). I then took the protein sequences of the 3 assemblies and aligned them with clustalo. They're all quite similar:

```
CLUSTAL O(1.2.4) multiple sequence alignment

NODE1_KJHFFCBH_00001         ------MEKRPEIQIFVNGHPIKFLLDTGADITILNRKEFIIGNSIENGKQNMIGVGGGK
MN630242_AMGANPBD_00002      -------------MIFVNGYPIKFLLDTGADITILNRRDFQVKNSIENGRQNMIGVGGGK
U11820.1_EANIPDKN_00002      MGTTTTLERRLEIQIFVNGHPIKFLLDTGADITILNRKDFQIGNSIENGKQNMIGVGGGK
                                           *****:*****************::* : ******:**********

NODE1_KJHFFCBH_00001         RGTNYINVHLEIRDENYKTQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISEKI
MN630242_AMGANPBD_00002      RGTNYTNVHLEIRDENYKTQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISDKI
U11820.1_EANIPDKN_00002      RGTNYINVHLEIRDENYRMQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISEKI
                             ***** ***********: **************************************:**

NODE1_KJHFFCBH_00001         PIVKVRMKDPTQGPQVKQWPLSNEKIEALTEIVERLEQEGKVKRADPNNPWNTPVFAIKK
MN630242_AMGANPBD_00002      PIVKVKMKDPNKGPQIKQWPLSNEKIEALTEIVERLEKEGKVKRADPNNPWNTPVFAIKK
U11820.1_EANIPDKN_00002      PIVKVRMRDPIQGPQVKQWPLSNEKIEALTDIVERLESEGKVKRADPNNPWNTPVFAIKK
                             *****:*:** :***:**************:******.**********************

NODE1_KJHFFCBH_00001         KSGKWRMLIDFRVLNKLTDKGAEVQLGLPHPAGLQWKKQVTVLDIGDAYFTIPLDPDYAP
MN630242_AMGANPBD_00002      KSGKWRMLIDFRELNKLTEKGAEVQLGLPHPAGLQMKKQVTVLDIGDAYFTIPLDPDYAP
U11820.1_EANIPDKN_00002      KSGKWRMLIDFRVLNKLTDKGAEVQLGLPHPAGLQMKKQVTVLDIGDAYFTIPLDPDYAP
                             ************ *****:**************** ************************

NODE1_KJHFFCBH_00001         YTAFTLPRKNNAGPGKRYVWCSLPQGWVLSPLIYQSTLDNILQPYIKQNPELDIYQYMDD
MN630242_AMGANPBD_00002      YTAFTLPRKNNAGPGRRYIWCSLPQGWILSPLIYQSTLDNIIQPFIRQNPQLDIYQYMDD
U11820.1_EANIPDKN_00002      YTAFTLPRKNNAGPGRRYVWCSLPQGWVLSPLIYQSTLDNILQPFIKQNSELDIYQYMDD
                             ***************:**:********:*************:**:*:** :*********

NODE1_KJHFFCBH_00001         IYIGSNLSRQEHKQKVEELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWSIQQKQL
MN630242_AMGANPBD_00002      IYIGSNLNKKEHKEKVGELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWTIQQKQL
U11820.1_EANIPDKN_00002      IYIGSNLSKKEHKQKVEELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWSIQQKQL
                             *******.::***:** ************************************:******

NODE1_KJHFFCBH_00001         EIPERPTLNDLQKLAGKINWASQTIPKLSIKALTHMMRGDQKLDSIREWTEEAKKEVQKA
MN630242_AMGANPBD_00002      DIPEQPTLNELQKLAGKINWASQAIPDLSIKALTNMMRGNQNLNSIREWTKEARLEVQKA
U11820.1_EANIPDKN_00002      EIPERPTLNELQKLAGKINWASQTIPDLSIKELTNMMRGDQKLDSIREWTVEAKREVQKA
                             :***:****:*************:**.**** **:****:*:*:****** **: *****

NODE1_KJHFFCBH_00001         KEAIEKQAQLYYYDPNRELYAKISLVGPHQLCYQVYHKNPEQILWYGKMNRQKKRAENTC
MN630242_AMGANPBD_00002      KKAIEGQVQLGYYDPSKELYAKLSLVGPHQISYQVYQRNPERILWYGKMSRQKKKAENTC
U11820.1_EANIPDKN_00002      KEAIEKQAQLNYYDPNRGLYAKLSLVGPHQICYQVYQKNPEHILWYGKMNRQKKKAENTC
                             *:*** *.** ****.: ****:*******:.****::***:*******.****:*****

NODE1_KJHFFCBH_00001         DIALRACYKIREESIVRIGKEPVYEIPASREAWESNIIRSPYLKASPPEVEFIHAALNIK
MN630242_AMGANPBD_00002      DIALRACYKIREESIIRIGKEPKYEIPTSREAWESNLINSPYLKAPPPEVEYIHAALNIK
U11820.1_EANIPDKN_00002      DIALRACYKIREESIIRIGKEPMYEIPASREAWESNLIRSPYLKAPPPEVEFIHAALSIK
                             ***************:****** ****:********:*.****** *****:*****.**

NODE1_KJHFFCBH_00001         RALSMVQEVPILGAETWYIDGGRRQGKAARAAYWTNTGRWQVMEIEGSNQKAEVQALLMA
MN630242_AMGANPBD_00002      RALSMIKDAPILGAETWYIDGGRKLGKAAKAAYWTDTGKWQVMELEGSNQKAEIQALLLA
U11820.1_EANIPDKN_00002      RALSMIQDAPITGAETWYIDGSRKQGKAARAAYWTDTGKWQIMEIEGSNQKAEVQALLLA
                             *****:::.** *********.*: ****:*****:**:**:**:********:****:*

NODE1_KJHFFCBH_00001         LKEGPEEMNIITDSQYILNIMNQQPDLMEGIWQEVLEEMEKKIAIFIDWVPGHKGIPGNE
MN630242_AMGANPBD_00002      LKAGPEEMNIITDSQYVINIILQQPDMMEGIWQEVLEELEKKTAIFIDWVPGHKGIPGNE
U11820.1_EANIPDKN_00002      LKAGSEEMNIITDSQYILNIINQQPDLMEGLWQEVLEEMEKKIAIFIDWVPGHKGIPGNE
                             ** * ***********::**: ****:***:*******:*** *****************

NODE1_KJHFFCBH_00001         EVDKLCQTMMIIEGDGILNKRPEDAGYDLLAAQEIHILPGEVRIVPTRTRIMLPKGYWGL
MN630242_AMGANPBD_00002      EVDKLCQTMMIIEGDGILDKRSEDAGYDLLAAKEMHLLPGEVKVIPTGVKIMLPKGYWGL
U11820.1_EANIPDKN_00002      EVDKLCQTMMIIEGEGILDKRSEDAGYDLLAAQETHFLPGEVRIVPTKTRIMLPKGHWGL
                             **************:***:** **********:* *:*****:::** .:******:***

NODE1_KJHFFCBH_00001         IMGKSSIGNKGMDVLGGVIDEGYRGEIGVIMINLSKKSTTILEKQKVAQLIILPCKHESL
MN630242_AMGANPBD_00002      IIGKSSIGSKGLDVLGGVIDEGYRGEIGVIMINLSRKSITLLEQQKIAQLIILPCKHEVL
U11820.1_EANIPDKN_00002      IMGKSSIGSKGMDVLGGVIDEGYRGELGVIMINLTKKSITILEKQKVAQLIILPCRHESL
                             *:******.**:**************:*******::** *:**:**:********:** *

NODE1_KJHFFCBH_00001         EQGEIIMNSERGEKGFGSTGVFSSWVDRIEEAELNHEKFHSDPQYLRTEFNIPRIVAEEI
MN630242_AMGANPBD_00002      EQGKVIMDSERGDKGYGSTGVFSSWVDRIEEAEINHEKFHSDPQYLRTEFNLPKMVAEEI
U11820.1_EANIPDKN_00002      QQGEIQMDSERGEKGFGSTGVFSSWVDRIEEAELNHEKFHSDPQYLRTEFNLPRIVAEEI
                             :**:: *:****:**:*****************:*****************:*::*****

NODE1_KJHFFCBH_00001         KRKCPLCRIRGEQVEGKLKIGPGIWQMDCTHFNGKIIIVAIHVESGLLWAQIIPQETADC
MN630242_AMGANPBD_00002      RRKCPVCRIRGEQVGGQLKIGPGIWQMDCTHFDGKIILVAIHVESGYIWAQIISQETADC
U11820.1_EANIPDKN_00002      KRKCPLCRIRGEQVGGQLKIGPGIWQMDCTHFNGKIIIVAVHVESGFLWAQIIPQETAEC
                             :****:******** *:***************:****:**:***** :***** ****:*

NODE1_KJHFFCBH_00001         TLKAIMQLVSTHNVTEIQTDNGPNFKNQKIEGLLSYMGIKHKLGIPGNPQSQALVENANS
MN630242_AMGANPBD_00002      TVKAVLQLLSAHNVTELQTDNGPNFKNQKMEGALNYMGVKHKFGIPGNPQSQALVENVNQ
U11820.1_EANIPDKN_00002      TVKALLQLICAHNVTELQTDNGPNFKNQKMEGLLNYMGIKHKLGIPGNPQSQALVENANN
                             *:**::**:.:*****:************:** *.***:***:**************.*.

NODE1_KJHFFCBH_00001         TLKVWIQKFLPETTSLDNALALALHCLNFKQRGRLGKMAPYELYTQQESLRIQDYFSQLP
MN630242_AMGANPBD_00002      TLKAWIQKFLPETTSLENALALAVHCLNFKQRGRIGGMAPYELLAQQESLRIQEFFSKIP
U11820.1_EANIPDKN_00002      TLKAWIQKFLPETTSLDNALALALHCLNFKQRGRLGKMAPYELYIQQESLRIQDYFSQIP
                             ***.************:******:**********:* ******  ********::**::*

NODE1_KJHFFCBH_00001         QKLMMQWVYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIRRVPEPCTLPE
MN630242_AMGANPBD_00002      QKLQAQWIYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIKRVPEPCALPE
U11820.1_EANIPDKN_00002      QKLMMQWVYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIRRVPEPCTLPE
                             ***  **:*****************************************:******:***

NODE1_KJHFFCBH_00001         GDE
MN630242_AMGANPBD_00002      GDE
U11820.1_EANIPDKN_00002      GDE
                             ***

```

Overall, the level of homology is quite high at the protein level. There's maybe a little evidence that this contig's protein is closer to U11820.1 than to MN630242 (though there are also several positions where U11820.1 is the outgroup). So I'm not sure why pieces of this contig are so hard to align at the DNA level without further investigation, but honestly overall of these are looking like very similar viral sequences, as might be expected.

I also aligned the big contig to both U11820.1 and MN630242 at the DNA level in 3-way `clustalo` alignment it is ok. Not sure what the difference is there. There is no obvious reason to prefer one reference over the other, according to my eye. When I have `clustalo` output a clustering solution the big contig is essentially equidistant between the two references, so I don't buy that U11820.1 is "better" as a reference:

```
(
MN630242.1:0.116151
,
(
NODE1_rc:0.0504856
,
U11820.1:0.0504856
):0.0656651
)
;

```

As I noted in the comments, it's unfortunate that QUAST and SPAdes aren't working very well. I don't know what the deal is there. But it seems like overall the assembly is quite close to both references. I don't have any intuition for why SPAdes is doing what it's doing, or why alignment at the DNA level is acting weird (lots of artificial gaps?). But I don't think it's a pathology of SPAdes particularly, it seems like it's just something weird that happened with this dataset interacting with DNA aligners. Possibly something went strangely with read deposition?

**Update:**

See OP's self-answer as well (and also discussion in the comments).

---

Based on the (extremely informative!) follow-up updates that you have provided, I think that we can tentatively answer your questions:

The tl;dr is that you can check your reads to make sure they aren't weird and you can check your output sequence to see if it's weird.

* *Am I using SPAdes correctly? Is there a better way to use it?*

1. I think that you are using SPAdes correctly, based on your experience getting ~finished viral genomes in the past with the same workflow. (THis is also my experience using SPAdes on viral genomes with Illumina PE data, is that it "just works".)
    2. It is possible that you are using *some other tools* in your workflow in a suboptimal way. For example, maybe your trimming is not quite correct (it does seem based on your answer that you are already trimming). The question in my mind is then whether there is some non-intuitive thing in the reads that should still be trimmed.
    3. I would suggest running FASTQC on your reads for this run and also your reads for the other runs that have worked, to see if something weird is going on. FASTQC directly checks for adapter sequence, weird GC composition, quality scores etc. That will give you a lot of information.
* *Is it possible that the sequences that blast finds matching with the SPAdes contigs are indeed sequences that someone got using SPAdes and that's why these match (not because these are correct)?*

1. This is formally possible, of course. Our main way of knowing that you're getting the right genome is... does it look like other genomes that have been generated in the same way?
    2. If you are really worried about this, I would suggest directly inspecting the sequence to see what the matches are. If you like, you can then blast the matches to see if they
    3. Another approach would be to annotate the genome and make sure that it "looks like" a phage genome. E.g. it has the expected genes from FIV. Prokka is very easy to run and works ok on phage.
    4. A slightly self-serving suggestion is to use PhageTerm to check your reads against your final genomes to ensure that they behave like phage- e.g. they have termini in the correct places, etc. (it does this with test coverage.) I am a contributor to PhageTerm, so I happen to know it a bit. Likely there are other options.
* *What other option could explain these circumstances that I could be missing?*

1. As I've suggested, I think that that leaves the reads. How sure are you that they are comparable to your other runs? They look like they should be sufficient from my glance at SRA, but maybe someone bumped the MiSeq while it was running. I've given suggestions above for checking with FASTQC.

# Answer
> 2 votes
After many considerations, I am going to accept the response from Maximilian Press. I see now that some viruses have high variability (HIV even 50% of the sequence). Therefore MN630242.1. and U11820.1 are apparently two strains. There are things I still don't understand but these are beyond the initial goal of my question. Concretely:

* Why SPAdes returns one strain and rnaSPAdes the other one.
* Why one strain matches 99,9% with MN630242.1 (and at least another assembly) while the other one has the (now) expected variability of this kind of virus.

I also want to point that apparently, Quast is not effective for viruses with such high variation.

I'm not going to delete this answer given that it responds to part of my question that is not responded in Maximilian's answer.

Therefore, directly answering my questions:

* **Am I using SPAdes correctly? Is there a better way to use it?** I am using SPAdes mostly in the correct way, similarly to how an average user would. But rnaSPAdes seems to be more appropriate for RNA viruses and it indeed works far better in this case.
* **Is it possible that the sequences that Blast finds matching with the SPAdes contigs are indeed sequences that someone got using SPAdes and that's why these match?** That could be possible in some cases but that doesn't mean these sequences are wrong as I initially considered. For whatever reason SPAdes seems to miss one strain (if that's what really is happening)
* **What other option could explain these circumstances that I could be missing?** See Maximilian Press answer.

---
Tags: assembly, covid-19, genome-sequencing, software-usage, virus
---


# Question
Title: Trying to create a .bam file without the need for a .sam file
I'm trying to use the code specified in this link to create a .bam file without the need for a .sam file.

Here is the code I'm using:

```
bwa mem -R "@RG\tID:$SAMPLE\tPL:illumina\tSM:$SAMPLE" \
    $REF_PATH/hg38.fasta \
    mapping/R1.fastq.gz \
    mapping/R2.fastq.gz | \
    samtools sort - | \
    tee mapping/$SAMPLE\.bam | \
    samtools index - mapping/$SAMPLE\.bam.bai \
    2>> log

```

and I get the following error:

```
[M::bwa_idx_load_from_disk] read 0 ALT contigs
`[E::main_mem] fail to open file ` '.

```

The ref path and file locations work when creating and converting a .sam to a .bam file. Any ideas?

# Answer
This is an error from `bwa` saying that it can't find your input files:

Here's a checklist:

1. Check that you have no whitespace after your `\` newline breaks. I.e. `\` is the *very last character*
2. Check that you have no whitespace in any path names
3. *Quote your variables*: `$REF_PATH` -\> `"$REF_PATH"`
4. Use `set -x` at the top of your bash script to debug.
5. In general, read this
> 5 votes

# Answer
Why make this so complicated with the tee command?

I don't think you even can index an unsorted bam, so you have to wait until sorting is done before you index. Indexing is really short compared to mapping and sorting, you aren't saving much by trying to do both at once.
> 0 votes

---
Tags: samtools, bash, bwa
---


# Question
Title: How to visualize called narrowPeak files in UCSC Genome browser?
I have this file:

and I get this error

I googled and found this solution but it still doesn't work and I get this error  Error File 'GSM2797523\_FOXA2\_IDR0.02\_narrowPeak\_try.bed' - line 1 of custom input: variable needs to start with letter '"Peak name=""Somite narrowPeak""" '

Any suggestions?

Many thanks

Found a way!!! opened it with text editor and not exel as before and added the line: track type=narrowPeak name="Somite narrowPeak"

# Answer
> 2 votes
Headers should not contain tabs:

```
track type=narrowPeak name="Somite narrowPeak" description="Somite narrowPeak"

```

Ensure that you have NO tabs on that line.

---
Tags: chip-seq, bed, ucsc
---


title: how could i find uniprot ac in blast?

text:
# Question
Title: how could i find uniprot ac in blast?
I have made a search in p-blast. I would like to search for the protein results in uniprot. how can I do it directly? I can go to GENBANK but i would like to go to UNIPROT. Thank you!!!

# Answer
> 1 votes
UniProt also provides a BLAST service (https://www.uniprot.org/blast), or you could use the Expasy BLAST too, at https://web.expasy.org/blast.

If you have performed your similarity search at NCBI or elsewhere, you can map the identifiers from that result to UniProt by using the UniProt IDmapping service at https://www.uniprot.org/uploadlists. This can also be used programmatically, cf https://www.uniprot.org/help/api\_idmapping

---
Tags: blast, uniprot
---


# Question
Title: how can i search in p-blast against UNIPROT(not only swissprot)
i need to search for one protein against UNIPROT in BLAST for a selected organism. How can I do it? (SWSISSPROT+TREMBL) THANKS!!

# Answer
> 0 votes
I don't know of any online tool enabling a search against TrEMBL. You probably will have to set up your own local BLAST environment, roughly following these steps:

1. download blast+ \- but there might be an easier way to install the software, depending on your operating system (e.g. on Ubuntu you can just `sudo apt install ncbi-blast+`)
2. download target databases
3. build local database

I hope this helps.

# Answer
> 0 votes
UniProt provides a BLAST service (https://www.uniprot.org/blast) which you can run against several UniProt databases, or you could use the Expasy BLAST too, at https://web.expasy.org/blast.

---
Tags: blast, uniprot
---


# Question
Title: >My counter is counting genotypic combination occurences more than once, how do I ensure it counts one combination and doesnt count it again?
```
mutations = ['A222V', 'D614G', 'E484Q', 'E780Q', 'G476S', 'L18F', 'N439K',
             'S477', 'S477N', 'T478I', 'V483A']

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

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

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

combination_labels = []
combination_counts = []

for lineage in lineages:

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

    for combination in combinations:

            combination_df = df[list(df)]

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

        #print(combination_df.shape[0])

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

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

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

```

The input CSV

```
lineage,A222V,D614G,E484Q,E780Q,G476S,L18F,N439K,S477,S477N,T478I,V483A
417941,0,1,0,0,0,0,0,0,0,0,0

```

Output CSV

```
combination,count,percentage
D614G,87355,90.7084929856806

```

My above code is used to count all occurences of combinations of spike protein mutations

* to determine the genotypes. However, looking at the values processed into the .csv file, i can see ' D614G' is being counted even when it is in combination with other combinations.

My question is how do i ensure once something has been counted (a row in the imported .csv file) it is not counted for a second time?

OR possibly how can I edit this code to prevent counting of singular mutations even when presented within a combination?

Thanks in advance

# Answer
Alright, so there are a number of problematic patterns in your code - as far as I understand what you are trying to do. Next time, try to post a reproducible example that people can use and more people will be willing to help.

### L18-21:

```
combination_labels = []
combination_counts = []

for lineage in lineages:

```

Declaring these two lists *before* the loop, then appending all your combination labels inside the loop means that these two lists will contain all labels and counts for all runs (all your lineages). If I am understanding your intentions correctly this is probably not something you want. The simple fix here is to just move them into the for-loop.

### L29-30:

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

```

Here you are iteratively checking if the value for each mutation is set to '1', but you are also *overwriting* your `combination_df` variable at each iteration of the loop. After this loop, `combination_df` will be whatever mutation last was set to '1' in your current combination. I will get back to solving this in my solution below.

### L37-44:

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

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

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

```

Since this part is again scoped inside the for-loop, each iteration (each combination) will overwrite the same `'Results_2/' + lineage.replace(".csv", "") + '_3.csv'` file. This you also want to move out of the `for combination in ...` loop.

---

A big concern is also performance. You have 2000-something combinations and 100,000 lines in your file. As you loop over combinations, you also process the DataFrame each time, when this is a proble you can solve in one iteration. It's not much data so it's still fine, but these habits are still good to develop early.

Now here's an idea for a rewrite that should work better. I don't have a csv from you to test with, so there might (proabably are) still things to fix, but hopefully it'll give you a starting point:

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

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

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

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

for lineage in lineages:

    out_df = pd.DataFrame({"Count": 0}, index=combinations)

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

    for index, row in df.iterrows():
        combination = '_'.join(sorted(df.columns[row==1]))
        out_df.loc[combination, "Count"] += 1

    out_df['Percentage'] = (out_df['Count'] / df.shape[0]) * 100
    out_df = out_df.sort_values('Percentage', ascending = False)

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

```
> 3 votes

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


# Question
Title: Are there computational tools to extract features of DNA sequences?
I am looking for tools to extract features from short DNA sequences. For example, entropy, complexity, GC-content, etc.

I have found the generateFeatures.py script from the PyFeat repo, but is there love a more widely used source code or a standard way to extract features from sequences biopython or similar?

Also, I think that since there are $2N$ encoded bits in a sequence of $N$ nucleotides, we have **at most** $2N$ independent features that could be extracted.

---

In addition, I am curios if there are any transformers models for DNA sequences.

---

***Edit***: In addition, you can use DeepHF's feature utils, which can be accessed and consumed as seen here.

# Answer
If you want to do this with biopython, the SeqUtils package could be a solution.
> 2 votes

---
Tags: python, phylogenetics, features
---


# Question
Title: Get list of urls of GSM data set of a GSE set
I have this GSE dataset ( GSE104279 ) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104279).

I want to make a table with set IDs and ftp urls to use it as a table in galaxy.org

I know that we can use ENA to get specific arrangement: https://www.ebi.ac.uk/ena/browser/

I tried to get:

```
SampleID    Group   URL

```

so I used : https://www.ebi.ac.uk/ena/browser/view/PRJNA412223

But nothing is showing.

Is there away to get these urls in the arrangement above?

# Answer
> 1 votes
You can use Entrez Direct for this as follows:

```
esearch -db gds -query 'GSE104279' \
  | esummary \
  | xtract -pattern DocumentSummary \
    -if 'entryType' -equals 'GSM' \
    -def 'NA' -element Accession title summary FTPLink

```

This will return a table with data similar to this:

```
GSM2580330      18_Z13_2_d5_Zika        cortical organoids_Zika_5d      ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2580nnn/GSM2580330/
GSM2580329      17_Z13_2_d5_Control     cortical organoids_mock_5d      ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2580nnn/GSM2580329/

```

To download sequence reads, you should follow links to SRA. Using Entrez Direct you can do this as follows:

```
esearch -db gds -query 'GSE104279' \
  | elink -target sra \
  | efetch -format runinfo

```

This will return a comma-delimited table containing SRA identifiers and an FTP path to the SRA data. These FTP paths won't lead you to FASTQ files though. You can pass the SRA run identifiers of the format SRR### (or ERR### or DRR###) to `fastq-dump` or `fasterq-dump` tools from the SRAToolkit to download data in FASTQ format.

---
Tags: sequence-alignment, sam, bioconductor, phylogenetics
---


# Question
Title: How to identify to each scaffold a read belongs to, inside a .sam file?
I have a fasta file assembly and combining it with the raw reads we produced a .bam file which I converted to .sam .

The .sam information lines look like this:

```
A00321:42:HLLVYDSXX:2:2302:6153:3505    99      NODE_1_length_3415511_cov_137.721502    16      60      128M    =       607     742     CGATTAGTCCGGCCAAATCGCCGTCGAGCGCAATGAACATAACGGTCTTGCCCTCAGCGCGCAGCGCATCGGCCTTGGCGTCGATTGTGGAGTGCTCGACGCCCATGATGTCCATCATAGCACCATTG        FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF        RX:Z:TTGAGGGTATAGTAGT   QX:Z:FFFFFFFFFFFFFFFF   TR:Z:GACACCG    TQ:Z:FFFFFFF    BC:Z:AGTTGCAG   QT:Z:FFFFFFFF   XS:i:-10        AS:i:0  XM:Z:0  AM:Z:0  XT:i:1  RG:Z:over_1kb:LibraryNotSpecified:1:unknown_fc:0        OM:i:60

```

Separated by mandatory fields it would be something like this:

```
QNAME: A00321:42:HLLVYDSXX:1:1644:2248:3881
FLAG: 99
RNAME: NODE_1_length_3415511_cov_137.721502
POS: 1
MAPQ: 60
CIGAR: 1S127M
RNEXT: =
PNEXT: 536
TLEN: 386
SEQ:  ATCGGGTCTGACACCGCGATTAGTCCGGCCAAATCGCCGTCGAGCGCAATGAACATAACGGTCTTGCCCTCAGCGCGCAGCGCATCGGCCTTGGCGTCGATTGTGGAGTGCTCGACGCCCATGATGTC
QUAL: FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

```

I'm actually interested in the meta data. I want to know how the RX: and BC: fields are distributed across the scaffolds in the original assembly.

I imagined the .sam file already contains the information about the assembly used to produce it. If I'm wrong, I'm sorry and please correct me, I'm just assuming.

What I want to do is, for each read in the .sam file, I find out its position in the assembled scaffold, and I record, Read\_ID,Scaffold\_ID,Read\_Position\_Inside\_Scaffold,RX,BC

Then I want to use that database to analyse the distribution of RX and BC inside each scaffold.

That's what I want.

Ultimately what I'm trying to do is evaluate the quality of my assemblies based on the Barcode distribution.

I'm good at programming and parsing, I'm just having trouble figuring out, where, inside the .sam file, can I find the scaffold and scaffold position of each read.

# Answer
> What I want to do is, for each read in the .sam file, I find out its position in the assembled scaffold, and I record, Read\_ID,Scaffold\_ID,Read\_Position\_Inside\_Scaffold,RX,BC

> I'm just having trouble figuring out, where, inside the .sam file, can I find the scaffold and scaffold position of each read.

You already listed all those:

* ReadID -\> `QNAME`
* ScaffoldID -\> `RNAME`
* Position -\> `Pos`

Assembly evaluation is a wide topic, but other software you can use includes BUSCO, QUAST or something like LTRi. Here's a more in-depth guide by the SciLife lab in Sweden. Usually, when aligning reads back to an assembly, you care about the mapped insert size distribution compared to the theoretical insert size distribution, which tells you about collapsing or exapnding regions (*i.e.* structural misassemblies). This is something that tools like FRC\_Align do.
> 1 votes

---
Tags: sam, fasta, samtools
---


title: How to find adapter sequence

text:
# Question
Title: How to find adapter sequence
I have this GSE dataset ( GSE104279 ) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104279)

I want to run cutadapt, but how would I find the adapter sequence so I can run cutadapt?

# Answer
The definitive answer is the sequencing institute/lab. They know what protocol/chemistry etc. was used.

If you don't have access to that a number of tools check for known adapter sequences. Run e.g. FASTQC, which will tell you the proportions of adapter sequence found. Tools like TrimGalore can also autodetect the most common adapters.

UPDATE: I checked 2M reads from one run of your dataset and it seems to me that it was uploaded already preprocessed:

Check all your samples, but I don't see a need to do any preprocessing if all of them look like this.
> 4 votes

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


# Question
Title: How can I add for several bed files the header : track type=narrowPeak name=“narrowPeak” preferably in python ,can handle with R
I want to create custom tracks from these files

I can add the line : track type=narrowPeak name=“narrowPeak” manually by opening it with text editor: track type=narrowPeak name=“narrowPeak”  but I wonder if there is no simple way to do it with code- to open all files that end with bed and add them this line. Didnt find this option in bed tools, If it is possible to add the name extracting it from the file name it would be amazing.

Many thanks

# Answer
> 1 votes
I don't think you need to go to R or Python to do that. In bash

```
set -euxo pipefail
for f in *.bed
do
  bn="$(basename ${f})"
  on="${bn%%.bed}_with_header.bed"
  echo 'track type=narrowPeak name="narrowPeak"' > "${on}"
  cat "${f}" >> "${on}"
done

```

EDIT: I just read you want to set the name based on the filename.

EDIT2: And also to add an arbitrary directory

```
set -euxo pipefail

INDIR="/path/to/bedfiles"

for f in $(find "${INDIR}" -name "*.bed")
do
  bn=$(basename "${f}")
  on="${bn%%.bed}_with_header.bed"
  dn=$(dirname "${f}")
  echo "track type=narrowPeak name=\"${bn}\"" > "${dn}"/"${on}"
  cat "${f}" >> "${dn}"/"${on}"
done

```

---
Tags: python, chip-seq, bed, ucsc
---


# Question
Title: Calculating average mito.percentage for each cluster (seurat)
I have a tricky data set with cells that will have a higher percent of mitochondria genes than "typical" data sets. I would like to look at the mito percentage in each cluster without any filtering (*e.g.*, `subset(seurat_object, subset = nFeature_RNA > 200 & nFeature_RNA < 3000 & percent.mito < 10)`)

```
FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)

```

Done!

I would also like to calculate the average percent.mito for each cluster and place the value on the FeaturePlot/TSNE plots. I have no idea how to manipulate a Seurat object to calculate these values but I can probably figure out how to place those values on the graph -- probably with `geom_text`.

My Hail Mary approach is below:

```
p3 <- FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)
p3 <- lapply(X = p3, FUN = function(x) mean(p3))
CombinePlots(plots = p3)

```

My attempts yield the error below; I don't know if my approach is close to correct and therefore, I have spent little time working through the error message.

```
In mean.default(p3) : argument is not numeric or logical: returning NA
Execution halted

```

and

```
p3 <- FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)
p3 <- lapply(X = p3, FUN = function(x) mean(seurat_object$meta.data$seurat_clusters))
CombinePlots(plots = p3)

```

```
Error in `[[.Seurat`(x, i, drop = TRUE) :
  Cannot find 'meta.data' in this Seurat object
Calls: lapply ... lapply -> FUN -> mean -> $ -> $.Seurat -> [[ -> [[.Seurat

```

and

```
p3 <- FeaturePlot(seurat_object, features=c('percent.mito'), cols = c("gold", "firebrick4"), pt.size=0.5, label=TRUE, repel=TRUE)
p3 <- lapply(X = seurat_object, FUN = function(x) mean(seurat_object$meta.data$seurat_clusters))
CombinePlots(plots = p3)

```

Error in as.list.default(X) :
  no method for coercing this S4 class to a vector
Calls: lapply -> lapply -> as.list -> as.list.default```
</code></pre>

If anyone can provide some guidance -- you will be inducted into my personal hall of fame. :)

Thank you for your time.

P.S. I am using Seurat 3.1

Update:

```
md = surate_object@meta.data
md_mean <- ddply(md, .(seurat_clusters), summarize, mean_mito_per=mean(percent.mito))

```

This was not the approach I was hoping for but I have the values -- which when expanded to the length of the original df, I can insert it into a free slot in the seurat\_object.

If anyone has a better approach I would love to hear it.

UPDATE: Using @PPK's approach

```
# First the cluster annotation and the tsne embeddings are merged
label.df <- cbind(as.data.frame(seurat_object@meta.data), as.data.frame(seurat_object@reductions$tsne@cell.embeddings))

# using dplyr across to calculate the mean mitochondrial percentage and
# the median tsne values per cluster
label.df <- label.df %>% dplyr::group_by(seurat_clusters) %>% dplyr::summarise(dplyr::across(percent.mito, ~ mean(.), .names = "mean_{.col}"), dplyr::across(contains("tSNE"), ~ median(.)))

# running FeaturePlot and adding the label layer
FeaturePlot(seurat_object, reduction = "tsne", features = "percent.mito") +
  geom_text(data = label.df, aes(x = tSNE_1, y = tSNE_2, label = mean_percent.mito))

```

Error:

```
* object '.col' is a function.
ℹ Input `..1` is `dplyr::across(percent.mito, ~mean(.), .names = "mean_{.col}")`.
ℹ The error occured in group 1: seurat_clusters = "0".
Backtrace:
     █
  1. └─global::initial_clustering(BM0H, "BM0H")
  2.   └─`%>%`(...)
  3.     ├─base::withVisible(eval(quote(`_fseq`(`_lhs`)), env, env))
  4.     └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
  5.       └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
  6.         └─`_fseq`(`_lhs`)
  7.           └─magrittr::freduce(value, `_function_list`)
  8.             ├─base::withVisible(function_list[[k]](value))
  9.             └─function_list[[k]](value)
 10.               ├─dplyr::summarise(...)
 11.               └

```

# Answer
You can extract the necessary values and add them directly the plot as a second layer using `plot + geom_text()`. This is very similar to the inner workings of the `DimPlot` function with `label = TRUE` but allows you to use anything as label.

```
# I use dplyr v1.0.2 for all data frame manipulations
# First the cluster annotation and the tsne embeddings are merged
label.df <- cbind(as.data.frame(seurat_object@meta.data), as.data.frame(seurat_object@reductions$tsne@cell.embeddings))

# using dplyr across to calculate the mean mitochondrial percentage and
# the median tsne values per cluster
label.df <- label.df %>% dplyr::group_by(seurat_clusters) %>% dplyr::summarise(dplyr::across(percent.mito, ~ mean(.), .names = "mean_{.col}"), dplyr::across(contains("tSNE"), ~ median(.)))

# running FeaturePlot and adding the label layer
FeaturePlot(seurat_object, reduction = "tsne", features = "percent.mito") +
  geom_text(data = label.df, aes(x = tSNE_1, y = tSNE_2, label = mean_percent.mito))

```

The labels are positioned at the median position of the cells of a cluster. Typing `LabelClusters` (without brackets) in the R console will show you the labeling function used by Seurat that contains some additional ways of positioning the labels.
If you want to avoid overlapping labels you can use `geom_text_repel` from the `ggrepel` package.
> 2 votes

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


# Question
Title: Table error when combining counts columns
I have these files that I want to make table of (https://github.com/learnseq/RNAseqfiles01.git), the table that I want (since all the files have the same ID column),, I want merge the ID column with the count column from each file:

```
ID    File 1 count column     File 2 count column     File 3 count column   File 4 count column

```

I cleaned up the files from extra lines using:

```
grep -vE '^{' GSM2580324_counts.txt | grep -vE '^\\' | grep -v '^ *$' | sed '/^{/d'  | tr '\\' ' ' | sed '/}/d' | grep -vE '^__' > GSM2580324_counts04.txt

```

But when I try to make a table in R using this:

```
%R files = list.files(path=".", pattern="GSM")
%R file1 = read.table("GSM2580319_counts01.txt" , col.names=c("gene_id","GSM2580319_counts01.txt" ))
%R file2 = read.table("GSM2580320_counts02.txt" , col.names=c("gene_id","GSM2580320_counts02.txt" ))
%R file3 = read.table("GSM2580323_counts03.txt", col.names=c("gene_id","GSM2580323_counts03.txt" ))
%R file4 = read.table("GSM2580324_counts04.txt", col.names=c("gene_id","GSM2580324_counts04.txt"))
%R out.file = merge (file1, file2, by=c("gene_id"))
%R out.file2 <- merge(out.file, file3, by=c(\"gene_id\"))
%R out.file3 <- merge(out.file2, file4, by=c(\"gene_id\"))
%R outtable01 <- write.table(out.file, file = "out.file3",sep="\t", row.names = FALSE)
%R outtable01

```

I get this error:

```
RParsingError                             Traceback (most recent call last)
<ipython-input-179-2c471cf65113> in <module>
      7 get_ipython().run_line_magic('R', 'file4 = read.table("GSM2580324_counts04.txt", col.names=c("gene_id","GSM2580324_counts04.txt"))')
      8 get_ipython().run_line_magic('R', 'out.file = merge (file1, file2, by=c("gene_id"))')
----> 9 get_ipython().run_line_magic('R', 'out.file2 <- merge(out.file, file3, by=c(\\"gene_id\\"))')
     10 get_ipython().run_line_magic('R', 'out.file3 <- merge(out.file2, file4, by=c(\\"gene_id\\"))')
     11 get_ipython().run_line_magic('R', 'outtable01 <- write.table(out.file, file = "out.file3",sep="\\t", row.names = FALSE)')

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line, _stack_depth)
   2315                 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals
   2316             with self.builtin_trap:
-> 2317                 result = fn(*args, **kwargs)
   2318             return result
   2319

<decorator-gen-160> in R(self, line, cell, local_ns)

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
    185     # but it's overkill for just that one bit of state.
    186     def magic_deco(arg):
--> 187         call = lambda f, *a, **k: f(*a, **k)
    188
    189         if callable(arg):

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/ipython/rmagic.py in R(self, line, cell, local_ns)
    729             if line_mode:
    730                 for line in code.split(';'):
--> 731                     text_result, result, visible = self.eval(line)
    732                     text_output += text_result
    733                 if text_result:

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/ipython/rmagic.py in eval(self, code)
    262             try:
    263                 # Need the newline in case the last line in code is a comment.
--> 264                 value, visible = ro.r("withVisible({%s\n})" % code)
    265             except (ri.embedded.RRuntimeError, ValueError) as exception:
    266                 # Otherwise next return seems to have copy of error.

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/robjects/__init__.py in __call__(self, string)
    385
    386     def __call__(self, string):
--> 387         p = rinterface.parse(string)
    388         res = self.eval(p)
    389         return conversion.rpy2py(res)

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/rinterface_lib/conversion.py in _(*args, **kwargs)
     38 def _cdata_res_to_rinterface(function):
     39     def _(*args, **kwargs):
---> 40         cdata = function(*args, **kwargs)
     41         # TODO: test cdata is of the expected CType
     42         return _cdata_to_rinterface(cdata)

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/rinterface.py in parse(text, num)
     44         raise TypeError('text must be a string.')
     45     robj = StrSexpVector([text])
---> 46     return _rinterface._parse(robj.__sexp__._cdata, num)
     47
     48

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/rinterface_lib/_rinterface_capi.py in _parse(cdata, num)
    556         if status[0] != openrlib.rlib.PARSE_OK:
    557             raise RParsingError('Parsing status not OK',
--> 558                                 PARSING_STATUS(status[0]))
    559     return res

RParsingError: Parsing status not OK

```

What should I do?

# Answer
> 1 votes
Regarding your `grep`, why not just grep for what you want, *i.e.* `^ENSG.*`? (These last few lines that htseq-count outputs - which is the software that produces your files - are really important though. You should absolutely look into them to see if something strange is going on with one or more samples.)

Your error comes from the escaped quotes `\"` you use here:

```
out.file2 <- merge(out.file, file3, by=c(\"gene_id\"))

```

Before, you used:

```
out.file = merge (file1, file2, by=c("gene_id"))

```

which is the correct way to quote here. On a sidenote, you don't need `c()` if it's only a vector with one element. You can just write `by = "gene_id"`.

Here's a different way you could read the files:

```
files <- list.files(path=".", pattern="GSM")
count_files <- lapply(files, read.table, col.names = c("GeneID", "Count"))
count_matrix <- do.call(cbind, lapply(count_files, "[[", "Count"))
dimnames(count_matrix) <- list(count_files[[1]]["GeneID"], files)

```

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


# Question
Title: Mining documents for nucleotide and peptide sequences
The text mining literature has an emphasis on identifying and normalizing gene names, mutations, pathways, concepts, and so on. I haven't been able to find much, however, on methods for extracting nucleotide and peptide sequences from documents. I did find this tidbit from Aerts *et al.* (2008), emphasis mine:

> Text was split into words and words greater than 10 characters in length with greater than 40% of characters from the capitalized DNA alphabet \[ACGT\] were extracted using regular expressions to isolate putative DNA sequences. All putative DNA sequences extracted from each paper were concatenated in the order they appeared in the text into a single fasta sequence and labeled with the corresponding PMID. Concatenation of sequences was performed to merge sequences split by line breaks in the text conversion, and **because we reasoned that inappropriate joins would be reconciled at the genome level by local alignment procedures**. Extracted, concatenated sequences were used as queries to BLAST RepeatMasked versions of genome sequences downloaded from the UCSC genome database...

I'm struck by how simple the author's method for extracting DNA sequences is. But I guess with the limited complexity of the DNA alphabet, that makes sense. This approach wouldn't work with the much more complex protein alphabet.

Are simple regular expressions the text mining community's state of the art for finding DNA/RNA sequences in documents? What methods are used to identify protein sequences?

# Answer
> 2 votes
I am not sure there is a state of the art for that.

Some possible strategies I can think of:

* Applying a dictionary to identify language words and discard these.
* Using some context information like identifying some words in the same sentence than the sequence.
* I guess the format in papers for sequences will not be the same as for text. Maybe you can find a pattern.

I would use an incremental strategy. Use the simplest algorithm then see what it recovers. If there is garbage, find a strategy to remove some garbage with a common pattern. Repeat until you are happy with the result.

Disclaimer: I have some experience in Natural Language Processing.

---
Tags: proteins, phylogenetics, text-mining
---


title: How do I create Protein Ligand images

text:
# Question
Title: How do I create Protein Ligand images
I was wondering how I can create images showing the bonds and the residue codes involved in the bound Protein Ligand complex obtained after Molecular Docking as depicted in the image below.

# Answer
> 3 votes
Pymol, in either its commercial or open source edition, will be sufficient to make figures such as your examples. Pymol recognises "polar contacts" for the intermolecular interactions and has labels that can be used for residue codes. It is also scriptable and is widely used.

# Answer
> 3 votes
Yes, PyMol is one of the best programs for rendering IMO. Usually, you would create folder with your PDB file and a file that contains all the commands to create your image like shown here. You can create more than one file. Because PyMol is based on Python you can use all python syntax here. E.g. you can create multiple images and with for loops you can create series of images that you can combine to a movie later on.

Other programs to check out would be VMD. I use VMD for quick videos and images. The PyMol images I like better most of the time.

Or newer software like Samson modeling platform. I am not sure about the licensing conditions. Have not tried it yet.

Biopython does not support the generation of 3D molecular images. It is more a library for sequence data.

# Answer
> 2 votes
The other contributors have given great suggestions. I thought I ought to mention Michelaɴɢʟo, which is an online PDB file/PDB code/PyMOL-PSE/gene-name visualiser, where one can create an interactive page, add a description with special protein links (*prolinks*) that change the interactive protein view and then share the URL securely.

It has a modal to make the selection easier, but uses the NGL selection language and to show a ligand, one would add a residue type prolink selecting `ligand`. If the neighbourhood-showing radius setting is greater than 1.0 Å (default is 3.0 Å) one will likely see hydrogen bonds with the neighbours.

One can try it without registering (free), but cannot keep the page for longer than a month or edit the description.

**Disclaimer**: I am the author.

---
Tags: python
---


# Question
Title: How to fix RepeatMasker.lib.nsq missing error when running RepeatModelor?
I have run RepeatModelor and I don't see the classified output file. I got this error.

```
Missing /home/mudith/.conda/envs/braker/share/RepeatMasker/Libraries/RepeatMasker.lib.nsq!
Please rerun the configure program in the RepeatModeler directory before running this script.

```

I think I have to run configure file. I found it. But it doesn't run when I just try `perl ./configure` Can someone please explain me how to run this? I'm new to this. Also I found this in github. But I don't understand.

```
cd $(dirname $(which RepeatMasker))/../share/RepeatMasker
# ./configure downloads required databases
echo -e "\n2\n$(dirname $(which rmblastn))\n\n5\n" > tmp && ./configure < tmp

```

I don't understand the third line.

# Answer
**What that code snippet is doing:**

That snippet you found on github (where is that from exactly?) looks like it tries to run RepeatMasker's configure script, feeding it some options that are needed by RepeatModeler. The cryptic third line generates and dumps a block of text to a file named tmp:

```

2
path/to/bin/rmblastn

5

```

And then passes that in to RepeatMasker's own configure script as though it were entered on the keyboard. It's the equivalent of hitting enter to accept a default, entering "2" to select the RMBlast option, typing in the path to rmblastn, hitting enter to set RMBlast as the default, and then typing "5" to save and quit.

That all looks fine, but I think this might be a problem with your (and mine, see below) RepeatMasker install rather than RepeatModeler itself. These tools are new to me but I'm interesting in seeing this work too, so I gave a shot at troubleshooting it here.

**The way I was able to get it working (note I assume you're working with conda packages here) was:**

1. `pushd $CONDA_PREFIX/share/RepeatMasker`
2. `./util/buildRMLibFromEMBL.pl Libraries/RepeatMaskerLib.embl > Libraries/RepeatMasker.lib` to "automatically create the old RepeatMasker.lib FASTA style database" as the script's comments put it. **NOTE: you might want to first delete or move the old empty file rather than alter it because of how conda shares data between environments.** If you update it in place (which `>` will still do), you'll alter the file for all conda environments using the package, and conda may complaint about unexpected file sizes when checking packages.
3. `perl configure` to do RepeatMasker's setup including some makeblastdb commands using RepeatMasker.lib. (Note this is RepeatMasker's configure script and not RepeatModeler's.) You can type in the options manually or you could probably use that code snippet you found.
4. `popd` and run RepeatModeler/RepeatClassifier and hopefully it finds the BLAST library files now

**All the gritty detail**:

I see a conda path in your logs. Are you using the bioconda packages named repeatmodeler and repeatmasker? It looks to me like the repeatmasker package in particular includes an empty RepeatMasker.lib file and that's the cause of the trouble. I don't see any way that RepeatMasker.lib is set up automatically via the configure script, and the BLAST database files like RepeatMasker.lib.nsq that RepeatModeler expects are built from that (empty) file.

The `Please rerun the configure program in the RepeatModeler directory` message is a red herring, I think, in this case because the `$CONDA_PREFIX/share/RepeatMasker/Libraries/RepeatMasker.lib` file *does* exist, just not the other files for a database it assumes should exist. If I instead run `configure` for RepeatMasker *and* update its script to not hide error messages, I see this coming from makeblastdb:

```
Building a new DB, current time: 10/26/2020 16:10:07
New DB name:   /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib
New DB title:  /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib
Sequence type: Nucleotide
Keep MBits: T
Maximum file size: 1000000000B
BLAST options error: File /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib is empty

```

And it's true, that's an empty file provided by the conda package (the metadata in `$CONDA_PREFIX/conda-meta/repeatmasker-4.1.0-pl526_0.json` shows `size_in_bytes: 0` for that one so that's apparently "correct" for the package.)

I saw a reference to that filename in a buildRMLibFromEMBL.pl script in the RepeatMasker share directory, and running that produces a FASTA file:

```
./util/buildRMLibFromEMBL.pl Libraries/RepeatMaskerLib.embl > Libraries/RepeatMasker.lib

```

Then `./configure` sees a non-empty .lib file and builds the database, and `RepeatClassifier` finally sees the resulting database files.

**But why is this needed?**

Shouldn't RepeatMasker come with a non-empty FASTA file if it's needed by its own configure script? I can see bioconda's build.sh doesn't do anything in particular to handle that file, but, maybe it should? This is probably where somebody who knows these tools better can chime in.)
> 2 votes

---
Tags: genome, sequence-annotation, repeatmasker
---


# Question
Title: Are codons in RNA layered? Are we misinterpreting RNA codons?
I am analyzing nucleotide base-pair patterns in RNA and DNA, and had a thought about RNA and DNA (Let me first state though, I am not a biologist; I am an algorithmatician, employing a sort of recursive sequential neural network).

In any event, I was wondering if we are misinterpreting the way that we look at codons in RNA. Are we? Are we missing something when we interpret them the way that we do? A picture is worth a thousand words, so I made a powerpoint slide of my question and screen captured it. It highlights the specific question I am asking, which is: are codons layered?

Why is this important to me? When I look at patterns, I know for a fact that unless it is the end of the sequence (a stop codon) that if I have a U-nucleotide, then the next two nucleotides cannot be AA, AG, nor GA. Additionally, I know that unless there is a U-nucleotide, the sequence will definitely not end in two or less nucleotides. This is what got me to wondering if there indeed is layered information in the combination of nucleotides. Please see the picture.

So if we have (as in the picture) an RNA sequence: AUGGUCAGUCCAUAA

We might interpret that as 5 codons {Methionine, Valine, Serine, Proline, STOP}

But what if there are actually 13 codons: {Methionine, Tryptophan, Glycine, Valine, Serine, Glutamine, Serine, Valine, Serine, Proline, Histidine, Isoleucine, STOP}

# Answer
This relates to the concept of the ***reading frame***, which is the point at which you start dividing 3-letter sequences into codons. As you point out, there are multiple ways to divide a string of amino acids into triplets - in fact, there are six unique ways, found by shifting the start letter by one three times, and doing the same but reading in reverse.

Typically, there is one commonly used frame which can be transcribed and turned into protein, called the open reading frame. But it is not unheard of for a single piece of DNA to encode multiple, overlapping genes, depending on where you start reading it. So, the sequence you show cannot encode 13 amino acids in a single string, but it can encode three different strings, one of which has 5 amino acids and the other two which have 4 amino acids each.

See Multiple Reading Frames at Wikipedia, which shows the following example of a single piece of human DNA which encodes two different genes, MT-ATP6 and MT-ATP8.
> 3 votes

# Answer
To add complexity to this discussion, there are examples in some mitochondrial genomes where the stop point of a gene is based on the polyA-attached transcript, rather than the genomic sequence. I worked on a paper about *Nippostrongylus brasiliensis* mtDNA, and needed polyA completion to explain the transcript expression I observed:

https://doi.org/10.12688/f1000research.10545.1

There are plenty of other examples of gene transcription happening in a situation other than one-exon-per-fragment: transcription that's on both sides of the DNA, exon sharing for different genes present on the same strand (they're usually in the same frame, but there's the MT-ATP6/MT-ATP8 example that @Nuclear\_Hoagie pointed out), micropeptide expression from smaller chunks of genes. Some also have some degree of stop-codon read-through, so the stop codons are not a strict stop in all cases. Some viruses are also particularly good at exploiting transcription / translation complexity.
> 1 votes

---
Tags: sequence-alignment, sequence-analysis, rna, computational-biochemistry, codon
---


title: Why use "robust" FPKMs?

text:
# Question
Title: Why use "robust" FPKMs?
Both DESeq2 and edgeR have an FPKM/RPKM function that by default uses normalized library sizes ("robust" option in DESeq2). FPKMs have their own issues, but I thought the main benefit was to have comparable units across independent experiments. If you are adding an experiment-specific normalization factor, then why even use FPKMs?

# Answer
> 5 votes
FPKM are inherently experiment specific and can not be used to compare across samples. Let's consider the following two sequencing runs. Let $E1$ and $E2$ be the true, underlying expression in two samples of genes 1-6. Let $S1$ and $S2$ be the observed expression in our sequencing.

$$ \begin{matrix} Gene & E1 & S1 & E2 & S2 \\ G1 & 100 & 10 & 100 & 20 \\ G2 & 100 & 10 & 100 & 20 \\ G3 & 100 & 10 & 100 & 20 \\ G4 & 100 & 10 & 0 & 0 \\ G5 & 100 & 10 & 0 & 0 \\ G6 & 100 & 10 & 0 & 0 \\ \end{matrix} $$

Our totals are 60 counts for $S1$ and 60 for $S2$. We sequenced both libraries to the same depth, but as there are fewer genes expressed in $E2$, we just naturally capture the genes that are expressed more often. For simplicity's sake, let's assume the genes all have the same length of 1, so we ignore the "K" in FPKM. Let's also forget the scaling by a million, since that's just to make nice numbers:

$$ \begin{matrix} Gene & E1 & S1 & E2 & S2 \\ G1 & 100 & 0.167 & 100 & 0.333 \\ G2 & 100 & 0.167 & 100 & 0.333 \\ G3 & 100 & 0.167 & 100 & 0.333 \\ G4 & 100 & 0.167 & 0 & 0 \\ G5 & 100 & 0.167 & 0 & 0 \\ G6 & 100 & 0.167 & 0 & 0 \\ \end{matrix} $$

Now let's calculate size factors. In DESeq2, we create a pseudo-sample which is the geometric mean of all counts of each gene, so $\sqrt{10 \cdot 20} = 14.14$ for genes 1-3 and genes 3-6 are ignored as they have a member with 0-valued counts. The size factor for each library is the median ratio for each library compared to that pseudo-sample:

$$ med(\[\frac{10}{14.14}, \frac{10}{14.14}, \frac{10}{14.14}) = 0.71 $$

$$ med(\[\frac{20}{14.14}, \frac{20}{14.14}, \frac{20}{14.14}\]) = 1.41 $$

So where do we end up if we divide by the size factors?

$$ \begin{matrix} Gene & E1 & S1 & E2 & S2 \\ G1 & 100 & 141 & 100 & 142 \\ G2 & 100 & 141 & 100 & 142 \\ G3 & 100 & 141 & 100 & 142 \\ G4 & 100 & 141 & 0 & 0 \\ G5 & 100 & 141 & 0 & 0 \\ G6 & 100 & 141 & 0 & 0 \\ \end{matrix} $$

So in this simple example, size factors produce a better estimate than FPKM or TPM to compare values *across experiments*, hence the "robustness". Please note that this is still not great for any statistical testing, and in general should be used with caution anywhere.

Finally, I'll leave some more reading material concerning FPKMs and cross-sample comparisons here:

---
Tags: rna-seq, deseq2, edger, fpkm
---


# Question
Title: Splitting data table up based on separate file of gene names
I am very new so bare with me. I have a file called data (.csv that I turned into a data table) with 2 columns (geneid, normalized counts) and about 9000 rows. I have a second file called Xgenes of genes that I am interested in pulling out of the first table (1 column- geneid, about 50 rows).

Is there a way to make 2 new files- I want one with all the geneid/normalized counts minus the ones in the second file (so minus about 50 genes). Then I want to make another file that only has the geneid/normalized counts of the 50 genes in the second file.

Currently I am doing this manually which does indeed work but if there is a way to do it through R it would make my future analysis much easier.

Thanks!

# Answer
There are plenty of efficient and smarter solutions, but if you have so small tables, you can just use the very basics of R to do what you need:

```
# you need to adjust the following two lines to fit to your file formating and file names
genes_filt_table <- read.table(..., col.names = c('geneid'), stringsAsFactors = F)
norm_counts <- read.table(..., col.names = c('geneid', 'norm_counts'), sep = ',', stringsAsFactors = F)

genes_to_filt <- norm_counts$geneid %in% genes_filt_table$geneid
filt_genes <- norm_counts[genes_to_filt, ]
non_filt_genes <- norm_counts[!genes_to_filt, ]

write.table(filt_genes, 'filtered_genes.csv', quote = F, row.names = F, sep = ',')
write.table(non_filt_genes, 'non-filtered_genes.csv', quote = F, row.names = F, sep = ',')

```

This being said, there are like milions of more efficient solutions out there. I provided one that is using as basic R functions as possible, as it seemed that you would like to see a simple R solution.
> 3 votes

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


# Question
Title: Number of edges in Graphical Analysis using Jupyter with Go Programming language
I'm trying to learn graphical analysis in a biological setting.

The online course suggested using the Go Programming language, so I'm new to it and help to understand why **the number of edges** that I read from the following code is different than from the site I'm downloading it from.

```
import (
    "fmt"
    "log"

    "github.com/kortschak/graphprac"
)

infile := "YeastL.dot"
g, err := graphprac.NewGraph(infile)
if err != nil {
    log.Fatalf("failed to read graph: %v", err)
}
{
    fmt.Printf("G has %d nodes and %d edges.\n", g.Nodes().Len(), g.Edges().Len())
}

```

The output I'm getting is: G has 2361 nodes and 6646 edges.

But from the website http://vlado.fmf.uni-lj.si/pub/networks/data/bio/Yeast/Yeast.htm it's suppose to read 7182 edges.

Thanks!

# Answer
the dot file in the Git repo has 6646 edges, the pajek file in your link has 7182.

Chances are the latter has been updated since.
> 2 votes

---
Tags: go, jupyter
---


# Question
Title: Any suggestions for cultivar identification using SSR (simple sequence repeat) markers
This is the data I have now: 30 simple sequence repeat (SSR) markers for 80 cultivars of cucumber. 10 of the 80 cultivars belong to one cultivar (let's say A).

My goal is to classify an unknown cultivar into A or "not A" using tge 30 SSR markers.

I'm think of considering it as a classification problem (A vs non-A) and use machine learning method to build a model using the SSR markers as features. But the problem is that A cultivars don't have enough number of samples.

Do you have any suggestions which statistical method(s) I can try to solve this problem? Thanks in advance.

# Answer
> 2 votes
Use Fst (migration statistic) such as the Fstat program to generate a distance matrix and solve it via clustering, such as UPGMA or neighbor-joining. This analysis is particularly useful for diploids. I presume you have heterozygotes in your analysis. You can switch Fis (inbreeding) if you have excessive homozygosity.

The other way in is a haplotype network. When I get a minute I'll post the link to the programs. Machine learning isn't needed in this instance, although you could do a PCA style analysis, i.e. unsupervised learning. Unsupervised learning doesn't need a large sample size.

---

The other solution is a stepwise mutation phylogeny, but this is for microsatellites (one to four repeats). This is called Delta Mu squared and describes the method of repeat unit mutation. It is a distance method.

Again when I have a moment I will supply weblinks.

# Answer
> 2 votes
Not to discount @Michael's answer (unsupervised learning is very handy for descriptive analysis of this kind of problem), but classification based on microsatellite data should not be too difficult. For papers with examples, you can see e.g. this and this. They seem to apply various algorithms, you might be able to learn something from those papers.

More generally, I bet that you could get pretty far with just (multinomial) logistic regression for your outcome variables.

Of course, this depends to a certain extent on how you code the data and what the input data look like. Microsatellite allele coding and evolutionary/error models are a field in themselves.

---
Tags: python, phylogenetics, phylogeny, machine-learning
---


title: Lipidomics normalization (peak area)

text:
# Question
Title: Lipidomics normalization (peak area)
I have data from a lipidomic experiment comprising: 2 conditions and 5 samples in each condition; and for each lipid the value is a peak area. I am doing a lipidomics analysis using the web app lipid ontology, and they recommend to normalize the data prior to doing the analysis,

> "There are many ways to normalize data, so LION/web cannot check whether data is properly normalized. If you don't know how to normalize the data, we recommend normalization by expressing all lipid species as fraction of the sum. A popular website that supports several normalization algorithms is MetaboAnalyst."

I don't know how to normalize my data, I cant find any paper related to lipidomics normalization to know to best way to do that, and also I can't find normalization function in the MetaboAnalyst.

Thank you in advance!

# Answer
Basically and briefly,

* Normalisation = transforming a data set to conform to a normal distribution. Converting to a frequency is NEVER a normalised transformation
* Standardisation = ensuring each point of the data set is comparable to the others, either directly or relatively depending on the calculation

The information you have received is incorrect, moreover you only convert to a frequency for a final assessment of the data and in my experience NEVER as an intermediate step. Normally data point between 0 and 1 are bad news and often require a transformation to get them away from fractions.

I don't understand the behaviour of the variance involved in metabolism and particularly lipids, however you should refer this issue to "cross-validation" Stackexchange. I did have some understanding of enzyme kinetics but in a different area.

Cross-Validation Stackexchance statisticians will understand less of the biology than I do (which is a bit worrying of course), so you will need to give them quite detailed information on the data set. They will almost certain support my concern of the advice you are given and provide a better method of standardisation. One method is to take the difference in mean and divide it by the standard deviation, you can then normalise the standardised data. They will want to know the goals of your investigation to understand who to process the data further.

What I can say is variation between runs is often an issue (in your apparatus perhaps not, I dunno) and standardising against a common control is one approach to overcome this, e.g. the control same common to all runs would be the mean about which all other data is centred.

What I suspect is that lipids are less popular than other fields, so standard protocol may not be observed.
> 1 votes

---
Tags: lipid-bilayer
---


# Question
Title: Reading fasta sequence regions with faidx
Im coding in c++ and reading in different reference genomes to examine regions across the chromosomes a few hundred basepairs at a time.

To do it in c++ i use the library htslib and the command faidx\_fetch\_seq() which is similar to read specific regions using samtools faidx.

The problem is when im reading in the reference genomes, there are of course multiple regions consisting only of NNNNN (and so on). Since i only need the actual DNA sequence i wish to read in the entire file but removing the regions with NNNN.

I have tried reading in the regions in c++ and then just creating an if-statement to check the actual bases consist of N. This works but is incredibly slow.

So my question is: do anyone know how to only read in the actual DNA sequences, removing the NNN regions? either with samtools or with the c++ funciton faidx\_fetch\_seq().

# Answer
> 1 votes
If you can't read in the entire sequence for a given reference genome:

(1) Use a script to generate a BED file of masked regions of low complexity (i.e. where you find `N`s in your FASTA file). (2) Use `bedops --complement` with this mask file and the output of UCSC `fetchChromSizes` to generate a BED file containing regions that are *unmasked*. (3) Repeat for other reference genomes, as needed.

You can now loop over the unmasked regions in this file, making calls to `samtools faidx` or `faidx_fetch_seq()`. These calls will return unmasked sequence.

For performance or convenience, you could read chunks of sequence into a `std::deque` or `std::vector` buffer and then read through that however many hundreds of bases, refreshing the buffer, as needed.

---
Tags: samtools, phylogenetics, c++
---


# Question
Title: How to sort multiple FASTA files based on their content?
So I have around 10,000 FASTA files of Influenza A virus. These files contains sequences of each of the 8 segments of the viral genome and I want to separate these files into different locations based on the content of these FASTA files. In each FASTA file, first line has the segment number. Eg.

```
>KM368312.1 Influenza A virus (A/swine/Shandong/01/2009(H1N1)) segment 3 polymerase PA (PA) and PA-X protein (PA-X) genes, complete cds

```

Okay imagine you have 10,000 text files in one folder and there are 8 types of these text files based on the content of them. Eg. file1.txt has segment 1 written in it, file2.txt has segment 2 written in it upto segment 8. So I want to group the files containing segment 1 together in one folder, segment 2 in other folder together and so on.

# Answer
> 2 votes
That's quite easy with `cmd`:

```
@echo off
for /l %%i in (1,1,8) do (
  for /f "delims=" %%f in ('findstr  /mc:" segment %%i " *.txt') do (
    ECHO move "%%f" "segment %%i\"
  )
)

```

`findstr /mc:"string" *.txt` searches each matching file for the string and outputs the filename, if the string was found. The `for /f` loop around catches those filenames, so you are able to move them. The `for /l` loop just iterates from 1 to 8.

See `for /?` and `findstr /?` for more information about those commands.

You probably want to adapt the filemask `*.txt` to your needs.

Note: I "disarmed" the `move` command for troubleshooting. Remove the `ECHO` when you are sure it does what you want.

# Answer
> 1 votes
The general idea is:

1. iterate over all files
    1. read the first line of each file, extract the `segment` bit
    2. move the file based on that

In Python, this can be done in many ways. Here’s one way:

```
import os
import re

def find_segment(id):
    return re.search(r'segment (\d)', id).group(1)

for segment in range(8):
    os.mkdir(f'segment-{segment + 1}')

for entry in os.scandir():
    if not entry.is_file():
        continue

    with open(entry.path, 'r') as file:
        segment = find_segment(next(file))
    os.rename(entry.path, f'segment-{segment}/{entry.path}')

```

Note the absence of error handling — this code assumes that every operation succeeds; empty files, or files without “`segment X`” in their FASTA ID, will trigger an error, as will usual IO failures (missing permissions, etc.).

---

But many people wouldn’t use Python to solve this but a shell scripting language. For example, here’s something equivalent in Bash (untested):

```
mkdir segment-{1..8}

for file in *.fa *.fasta; do
    segment=$(sed 's/.*segment \([[:digit:]]\).*/\1/; q' "$file")
    if [[ "$segment" != '' ]]; then
        mv "$file" "segment-$segment"
    done
done

```

---
Tags: python, covid-19, virus
---


# Question
Title: How to download sequence alignments against UniProtKB sequence database for a particular Pfam Id using R?
I am trying to download sequence alignments for families against the UniProtKB sequence database as shown in the example below using R code. Can someone kindly suggest me how to?

```
https://pfam.xfam.org/family/PF02171/alignment/uniprot/format?format=pfam&alnType=uniprot&order=t&case=l&gaps=default&download=1

```

# Answer
The following code served my purpose.

```
library(readr)
pfam_url <-
  'https://pfam.xfam.org/family/PF02171/alignment/uniprot/format?format=pfam&alnType=uniprot&order=t&case=l&gaps=default&download=1'
seqs <-
  read_table(
    pfam_url,
    col_names = FALSE,
    col_types = cols(X1 = col_character(), X2 = col_character())
  )

```
> 2 votes

---
Tags: r
---