# Question
Title: Output when compute GC content is different between using user input and using string
I am learning bioinformatics and i was solving a problem on rosalind.info. Input is a string consist of multiple DNA strand names and the DNA sequences, and Output is the name of the sequence that has the highest GC percentage and the GC percentage. this is the input:

```
>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT

```

this is the correct output:

```
Rosalind_0808
60.919540

```

This is my code (sorry for this filthy mess):

```
DNA = input()
wholeRosalind_DNA = DNA.splitlines()
Rosalind = []
Rosalind_DNA = []
for i in wholeRosalind_DNA:
    if "Rosalind" in i:
        Rosalind.append(i)
Countstring = []
maxGCcontent = []
GCcontent = []
def eachRosalind_DNA(DNA):
    for i in range(0,len(Rosalind)):
        if Rosalind[i] in wholeRosalind_DNA:
            a = wholeRosalind_DNA.index(Rosalind[i])
            if a == wholeRosalind_DNA.index(Rosalind[-1]):
                b = len(wholeRosalind_DNA)+1
            else:
                b = wholeRosalind_DNA.index(Rosalind[i+1])
            Rosalind_DNA.append(wholeRosalind_DNA[a:b])
    for i in Rosalind_DNA:
        Countstring.append(''.join(i))
    for i in Countstring:
        bases = len(i)-14
        CG = i.count("C") + i.count("G")
        GCcontent.append([i[1:14],(100*(CG/bases))])
    for i in GCcontent:
        maxGCcontent.append(i[1])
    result = max(maxGCcontent)
    for i in GCcontent:
        if i[1] == result:
            return (i[0]),(round(i[1],6))
print(eachRosalind_DNA(DNA)[0])
print(eachRosalind_DNA(DNA)[1])

```

but i get this instead:

```
  File "GC_Content.py", line 26, in eachRosalind_DNA
    GCcontent.append([i[1:14],(100*(CG/bases))])
ZeroDivisionError: division by zero

```

i will get the desired output if i just paste the input into the variabe "DNA" without using user input. Please help me fix the code and distinguish those two ways of getting input

# Answer
Is it cheating to use Biopython(?), will a solution might be helpful anyway. This assumes you have your sequences stored in the seq.fa FASTA file:

```
from Bio import SeqIO
from Bio.SeqUtils import GC

print(max(((r.id, GC(r.seq)) for r in SeqIO.parse("seq.fa", "fasta")),
          key=lambda x: x[1]))

```

Output:

```
('Rosalind_0808', 60.91954022988506)

```
> 3 votes

# Answer
The error is that `bases = 0`. Briefly looking at the code carriage returns in your original code will trigger this, as will confusing a faster header with sequence. The latter is because the fasta headers are 14 characters long and subtracting 14 from bases which is of course zero. The former will likely result in 'nan' because it would be 14 subtracted from ''.

I've a simpler solution, download Biopython (using `conda` or last resort `pip`), you might want to use SeqIO (Biopython), which would be cooler because it allows you to call Rosalind\_0808 directly. Then:

```
from Bio.SeqUtils import GC
seq = 'CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGACTGGGAACCTGCGGGCAGTAGGTGGAAT'
print (GC(seq))

```

Done. If you want to splice off the last 14 characters its

```
GC(seq[:-14])

```
> 0 votes

# Answer
I think the hard thing in this problem is not the calculation of GC content, it is the processing of DNA strings because they are separated on different lines. I would suggest reading the file line by line to generate complete DNA strings for each Rosalind entry. (I am not able to spot what is causing the error in the code).

```
seqDic = {}
with open("YouInFile", "r") as inFile:
    for ln in inFile:
        ln = ln.rstrip()
        if ln.startswith(">"):
            ln.name = ln[1:]
            seq = inFile.readline().rstrip() + inFile.readline().rstrip()
            seqDic[ln.name] = seq
        else
            continue

```

Basically, whenever you encounter a "\>" in your file, you write down the Rosalind name, and then you read the next two lines and combine them together. You put the name and the sequence in a dictionary, and iterate over the dictionary to calculate GC content.

Note that I did not debug my code. Let me know if it is not working.
> 0 votes

---
Tags: python
---


title: BigWig circular plot representation

text:
# Question
Title: BigWig circular plot representation
I would like to ask if anyone knows or has experience in representing BigWig signals in a circular plot? I mean like IGV but in circular plot not linear. Or any library that can be used in plotting BigWig signals especially in Javascript?

I appreciate any suggestion, Thank you.

# Answer
> 2 votes
circosJS should probably work well, based on your requirements: https://github.com/nicgirault/circosJS

I posted a short Observable notebook that shows what calls could be used to bring data from remote bigWig files into a circosJS plot:

https://observablehq.com/@alexpreynolds/rendering-bigwig-file-data-in-a-circular-representation

It is less intended as an explicit how-to in Javascript, than as a starting point for the end user with some of the function calls and data scale issues involved. Hopefully it can be useful to the reader for exploratory purposes.

---
Tags: bigwig
---


# Question
Title: How to analyze co-occurrence of multiple SNPs?
I am interested in 20 different SNPs that all are either As or Gs, and they all occur on the same chromosome. How can I assess the co-occurrence of these SNPs? In other words, I want to know, if SNP1 is "G", does that make SNP2 to be more likely to be a "G" as well?

I know this is broad and I'm not looking for a script or anything so specific, just a point in the right direction. I'm working with a non-model organism, so any organism-specific application won't be useful without being able to learn the methodology.

Googling thus far has only brought up GWAS and similar methods of correlating a genotype and a phenotype, but not within a genotype.

# Answer
> 1 votes
I found the LDheatmap package in R very helpful at calculating linkage disequilibrium. You can make plots such as this:

---
Tags: snp, gwas, correlation
---


title: How I reformat this data

text:
# Question
Title: How I reformat this data
I have this data all is one line even gene from first sample joined the second one while I want this data column wise each field separated like I have provided

```
#Uploaded variation   #Uploaded_variation 1 883477 G T 1_883477_G/TPRDM16 1 3125744 A G 1_3125744_A/GPRDM16

```

But I want

```
         V1     V2  V3      V4  V5  V6
1_883477_G/T    1   883477  G   T   PRDM16
1_3125744_A/G   1   3125744 A   G   PRDM16

```

# Answer
> 1 votes
I am really curious what software or piece of code gave you this one-line data. I would think twice before using a tool giving such output. Anyway, I think regular expressions is the way to go in this case (this solution uses `R`):

```
library(stringr)

text <- "#Uploaded variation   #Uploaded_variation 1 883477 G T 1_883477_G/TPRDM16 1 3125744 A G 1_3125744_A/GPRDM16"

my_pattern <- "[0-9]+\\s[0-9]+\\s[ATGC]+\\s[ATGC]+\\s[0-9]+_[0-9]+_[ATGC]+\\/[0-9A-Z]+"

info_as_list <- str_extract_all(text,
                                pattern = my_pattern) %>% unlist()

info_as_list <- sub("/",
                    " ",
                    info_as_list)

info_as_list <- str_split(info_as_list, " ")

as.data.frame(Reduce(rbind, info_as_list))

```

This is the output, I guess the rowname "init" comes from `Reduce()`:

```
     V1      V2 V3 V4          V5      V6
init  1  883477  G  T  1_883477_G TPRDM16
      1 3125744  A  G 1_3125744_A GPRDM16

```

Here is what the regular expression does, you might need to alter this based on your data:

```
[0-9]+\\s       :   start position
[0-9]+\\s       :   end position
[ATGC]+\\s      :   ref nucleotide(s)
[ATGC]+\\s      :   alternate nucleotide(s)
[0-9]+          :   start
_               :   underscore
[0-9]+          :   end
_               :   underscore
[ATGC]+         :   ref
\\/             :   slash
[0-9A-Z]+       :   gene symbol

```

`+` means one or more of the preceding and `R` uses two backslashes, `\\`, as escape characters.

---
Tags: awk
---


# Question
Title: What is the meaning of these misaligned reads in a sequencing run?
I am analyzing some SARS-CoV-2 sequencing runs abd often find read alignments like the one in the image.

```
GACCCGGTGGGTTTTAC---ACTTAAAAACACAGTCTGTAC
GACCCTGTGGGTTTTAC---ACTTAAAAACACAGTCTGTAC
          GTTTCAC---ACTTAAAAACACAGTCTGTAC
          GGTTTAC---ACTTAAAAGCACAGTCTGTAC
CCACACTCTCCTAGCAC---ACTTAAAAACACAGTCTGTAC
CCACACTCTCCTAGCACCATACTTAAAAACACAGTCTGTAC

```

It basically says there are reads that share part of a sequence and part of another one. The first 4 reads are the same sequence and it's a sequence that you can effectively find in the official SARS-CoV-2 assembly but the other two reads are a mix of this one with another sequence (and both are the same mix) but with only a small (3bps) indel of difference.

The sequence of the last two reads does not appear anywhere in the official SARS-CoV-2 sequence.

* What does it mean?
* Are these sequencing errors?
* If that is the case how is it that the same mistake is repeated with similar but not exactly the same data?
* If these are real RNA sequences present in the probe, why are these not used anywhere in the final assembly and how did they arrive there?

Is there any other option that I am missing?

# Answer
> 2 votes
The unaligned parts of the bottom two sequences are from elsewhere in the sars-cov-2 genome if you blast them. Likely causes of this are library prep. artifacts (not terribly common, but you'll find a few cases here in there once you have millions of reads) where fragments end up ligating together or (probably more likely) rearrangements of the viral sequence. My guess would be that the latter is what occurred, given how rapidly viruses reproduce. Whether that arrangement can lead to viable viruses I don't know (that's too far out of my field).

Assemblies represent at best a snapshot of a single example of a species at a single moment in time. Actual individuals in the wild can vary quite substantially, especially when they don't need to be viable.

---
Tags: sequence-alignment, ngs, assembly, covid-19, sars-cov-2
---


title: Hosting IGV web on a server

text:
# Question
Title: Hosting IGV web on a server
I have my analysis on a cloud server. The data is quite large, so copying it locally is not really an option. I want to be able to host IGV as a server on that server, so that it can access those files.

I installed IGV web (https://github.com/igvteam/igv-webapp) on my server, and ran it, but the Tracks → Local File option opens a file browser on *my* local machine, it doesn't let me select files on the server.

Is there any easy way to host IGV web on a cloud server?

# Answer
You can use X11 forwarding to run IGV on the headless server where the data is located and simply forward the GUI to your client machine using the X11 protocol. This is a solution I use consistently for my research when working with 20+ tracks that I want to keep on the server.

1. Download and install IGV for linux from the IGV site on to your server.
2. Connect to the server using ssh with X11 forwarding enabled

`ssh -X user@server`
3. Navigate to the recently downloaded and unzipped IGV folder. Execute the igv.sh script.

`./igv.sh`

The script launches the IGV session with 4gb of memory. Edit the script to launch with more, 32gb in this edited excerpt.

```
exec java -showversion --module-path="${prefix}/lib" -Xmx32g \
@"${prefix}/igv.args" \
-Dapple.laf.useScreenMenuBar=true \
-Djava.net.preferIPv4Stack=true \
--module=org.igv/org.broad.igv.ui.Main "$@"

```
> 7 votes

# Answer
According to their discription

> igv-webapp is a pure-client "genome browser" application based igv.js.

Pure-client, so it will be difficult to let it display tracks that are on the server itself. Perhaps you are looking for something more similar to the Apollo Genome browser. This, while more difficult to install, allows for the display of server-side data tracks and allows for manual curating of annotations as a nice bonus :)
> 1 votes

# Answer
One possibility might be to put the full IGV application on the server, and use the igvR package in R on the server to access the IGV API. You could then use R plotting facilities etc to show your data, accessed via the IGV API. There might be APIs for IGV in other languages, if you prefer, not sure.
Edit: I haven't tried this package, this is a speculative idea.
> 1 votes

# Answer
obviously I'm late but if this is still an open question open an issue in one of our github projects,
github.com/igvteam/igv (Java desktop), github.com/igvteam/igv.js (pure javascript component), github.com/igvteam/igv-webapp (html/javascript application, 100% client side) as appropriate.

As already mentioned the web igvs are 100% client, data is fetched from servers using range requests. It scales very well, as long as your data is indexed. If you really want to generate server side images a better choice might be the UCSC genome browser which is designed for that.
> 1 votes

---
Tags: sam, igv, cloud
---


# Question
Title: Is there a way to save a full pdb file in PYMOL after uploading biological unit?
pdb code: `2obt`

The pdb file has a coordinates of a monomer.

I have used following commands to generate an assembly, or in other words to upload one more biological unit so that the protein would be a dimer:

```
set assembly, ""
fetch 2obt, asu, async=0

set assembly, "1"
fetch 2obt, assembly1, async=0

```

After generating the functional dimer in pymol under name `assembly1`, I have an important issue that I cannot save a dimer structure in pdb file with dimer coordinates. Does anybody has an idea how to go around it or is it even possible to do that?

# Answer
> 4 votes
So on the GUI form of PyMOL on the right hand side you have a list of what is ambiguously called `objects` in PyMOL (models elsewhere), which can be enabled or disabled by clicking on their name. The name also has a `1/1` in green or similar. This indicates the state. In your case you have 2 states.

To split a multistate object use

```
split_state prot_A

```

To combine two objects one uses the command

```
create new_name, prot_A or prot_B

```

Now, this has two caveats which if violated result in a nasty piece of modern art:

* if the chains, segment ids and residues are the same
* if you had two overlapping macromolecules

To deal with the first, change the chain id on one of them.

alter prot\_B, chain='B' sort

Segment identifier is a weird thing that most comp-biochem tools dislike, so I like to remove that

alter prot\_B, segi='' sort

Now when you go to save your `new_name` object, in the popup of the GUI you get a dropdown menu defaulted to enabled (visible objects that have not be clicked to be disabled). There you can choose `new_name`.

## Aligning two

If someone that is in not the OP author has a template protein and needs to make a dimer of another based on that, the process is the same.

```
create twin, prot_B
align twin, prot_A and chain B
alter twin, chain='B'
sort

```

## SDF files

If someone else has a multicompound sdf where the state is the header name, I strongly advise not to use `split_state`. Yes, the python pymol module (not the app) handles states in a really weird way, but states can be selected with `state 1` selection algebra and `state` is an argument for most commands that do not accept `state` selector.

---
Tags: pymol
---


# Question
Title: Assembly reads having a copy of their beginning in their tail
I am analyzing the reads for the SARS-CoV-2 assembly with id SRR11140748. Apparently these reads were obtained with parallel sequencing by Illumina and Oxford Nanopore Technologies.

I have found these reads

```
GGTAAGTATGTACAAATACCTACAACTTGTGCTAATGACCCTGTGGGTTTTACACTTAAAAACACAGTCTGTACCGTCTGCGGTATGTGGAAAGGTTATGGCTGTAGTTGTGATCAACTCCGCGAACCCATGCTTCAGTCAGCTGATGCACAATCGTTTTTAAACGGGTTTGCGGTGTAAGTGCAGCCCGTCTTACACCGTGCGGCACAGGCACTAGTACTGATGTCGTATACAGGGCTTTTGACATCAGTATGTACAAATACCTACAACTTGTGCTAATGACCCTGTGGGTTTTACACTTAAAAACACAGTCTGTACCGTCTGCGGTATGTGGAAAGGTTATGGCTGTAGTTGTGATCAACTCCGCGAACCCATGCTTCAGTCAGCTGATGCACAATCGTTTTTAAACGAGTATGTACAAATACCTACAACTTGTGCT

```

As you can see, the second one has an exact copy at the tail of its beginning:

```
AGTATGTACAAATACCTACAACTTGTGCT

```

I have found other issues apparently caused by rearrangements of the virus genome, but this case seems to be different. It seems unlikely that the virus mutated in such a way and then a read was cut exactly this way. What is the most likely explanation for this read?

It's not the only one displaying the same or similar behavior. for example this other one has a copy at the tail but from bps not exactly at the beginning (CACCGGAAGCCAATATGGATC):

```
CAGTTACACCGGAAGCCAATATGGATCAAGAATCCTTTGGTGGTGCATCGTGTTGTCTGTACTGCCGTTGCCACATAGATCATCCAAATCCTAAAGGATTTTGTGACTTAAAAGGTAAGTATGTACAAATACCTACAACTTGTGCTAATGACCCTGTGGGTTTTACACTTAAAAACACAGTCTGTACCGTCTGCGGTATGTGGAAAGGTTATGGCTGTAGTTGTTACACCGGAAGCCAATATGGATC

```

# Answer
If you google AGTATGTACAAATACCTACAACTTGTGCT you'll find it's a primer sequence: https://artic.network/resources/ncov/ncov-amplicon-v3.pdf
> 2 votes

---
Tags: ngs, covid-19, illumina, genome-sequencing, sars-cov-2
---


# Question
Title: conda doesn't install latest version of snakemake
I would like to use the `--default-resources` parameter for profiles which is available in later versions of snakemake.

To install snakemake, I created a new conda environment:

`conda create -n snakemake -c conda-forge -c bioconda snakemake`

The currently installed version is listed as 5.14:

```
$ conda activate snakemake
(snakemake)$ conda list | grep snakemake
# packages in environment at /home/uvi/be/hde/.conda/envs/snakemake:
snakemake                 5.14.0                        0    bioconda
snakemake-minimal         5.14.0                     py_1    bioconda

```

However, command line tool says it is on version 5.1.4:

```
(snakemake)$ snakemake --version
5.1.4

```

This version does not know the new parameter `--default-resources` yet. So it gives an error when using it

> snakemake: error: unrecognized arguments: --default-resources partition = "thin-shared"

I already checked that I'm running the executable from inside the environment, rather than a different, external one:

```
(snakemake)$ which snakemake
~/.conda/envs/snakemake/bin/snakemake
(snakemake)$ $(which snakemake) --version
5.1.4

```

# Answer
> 4 votes
Python (whether installed with conda or not) will prefer packages installed in your home directory unless you tell it not to. This is a little known secret about using a virtualenv, actually, since internally those tell python to ignore everything in your home directory. Within conda, there's long been a debate about whether python (and also R) should continue looking for packages in the normal places or only within a conda environment. Since it would likely break many people's workflows, it's been decided to allow the normal language features of looking for packages not just within conda envs but also in the normal system and user locations. For python in particular you can disable this as follows:

```
PYTHONNOUSERSITE=True snakemake ...

```

Then only the packages installed within your conda environment will be used.

---
Tags: snakemake, conda
---


title: TEs vertical or horizontal transfers?

text:
# Question
Title: TEs vertical or horizontal transfers?
After predicting and annotating all the transposable elements. For four related species.

Now I would like to find out these predicted candidates found in my costume library, either vertical or horizontal transfers?

are there any methodologies to suggest?

# Answer
> 1 votes
You make a phylogenetic tree of the "host species" and compare it to the phylogeny of the TE's. Homology is not a good idea for evolutionary trees.

There are two snags:

1. Concerted evolution, the TE's will risk recombining with each other within a host genome and that can mess up the phylogeneic signal (the theory as why is quite complicated)
2. If host species phylogeny are higher eukaryotes thats a lot of genetic material, so you'd need a short-cut. You can either use someones species tree, or focus on TE's within a subsection of the host genome and produce a tree of this subsection.

With point number 2 "a subsection of the genome" you can then assess synteny, i.e. has the TE jumped between position between species. If the TE hasn't moved then its an orthogonal proof of putative vertical transmission.

What you are looking for is incongruence between host phylogeny and TE phylogeny as evidence of horizontal transmission and congruence to call vertical transmission. However, incongruence could occur due to the complexity of concerted evolution, so synteny would assist again as an orthogonal proof. Congruence is assessed as a "mirror tree", there is quite a lot of theory on this looking at host-parasite coevolution

For making trees RAxML or PhyML are flavours of the month particularly for TEs.

---

A state of the art-of-the-art review is here, but it goes into significant detail into precise proofs, which in my opinion arn't needed

---
Tags: transposable-elements
---


title: Extract end position from CIGAR

text:
# Question
Title: Extract end position from CIGAR
I was wondering how could I obtain the end position of an alignment using its start position and its CIGAR String.

If an alignment present some soft clipping, for example:

```

    20S55M

```

Then, would the start position indicated in the SAM file be the first soft clipped base, or the first correctly mapped base?

Thanks

# Answer
The start position is the first correctly mapped base. One way to do it, is to use GenomicAlignments in R to read in the alignments.

For example the sam file looks like this:

```
cat test.sam
@HD VN:1.3  SO:coordinate
@SQ SN:chrI LN:228861
D00689:415:CD4LUANXX:3:1316:7447:98083#CGGCTATGCAGGACGT 0   chrI    769 0   42M8S   *   0   0   CCATAGCCCGGTGGGGCCTGTGCTGCCACAGAGGGTGTGTCTAGATCGGA  B/B//FFFFFFBFFFF//7FF/<F/FFBFFFF//<F///<F/BFB//<F/  NM:i:1  MD:Z:21A20  AS:i:37 XS:i:37
D00689:415:CD4LUANXX:3:2308:9597:30445#CGGCTATGCAGGACGT 0   chrI    3965    0   3S47M   *   0   0   CCACACACCCACACACCACACCCACACCACACCCACACCACACCCACACC  BBBBBFFFFFFFFFFFFFFFFFFFFFFFFFFFFBFFFFFBF<FFFFF<FF  NM:i:0  MD:Z:47 AS:i:47 XS:i:47

```

We convert to bam:

```
samtools view -bS test.sam > test.bam

```

And read using R:

```
library(GenomicAlignments)
gr = readGAlignments("test.bam")
gr

GAlignments object with 2 alignments and 0 metadata columns:
      seqnames strand       cigar    qwidth     start       end     width
         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
  [1]     chrI      +       42M8S        50       769       810        42
  [2]     chrI      +       3S47M        50      3965      4011        47
          njunc
      <integer>
  [1]         0
  [2]         0
  -------
  seqinfo: 1 sequence from an unspecified genome

```

The start and the end can be retrieved using:

```
start(gr)
[1]  769 3965
end(gr)
[1]  810 4011

```
> 3 votes

# Answer
The match is starting with the first mapped base ignoring any softclipped parts of the read.

As an example this are two reads from a bam file of mine:

```
start   cigar   sequence
11710   17S83M  NNNNNNNNNNNNNNNNNGGCCTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCC
11713   100M    CTGGCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTTGAGAATGACTGCGCAAATTTGCCGGATTTCCTTTGCTGTTCCT

```
> 2 votes

---
Tags: r, sequence-alignment, sam, cigar
---


# Question
Title: SNPEFF genome build with specific codon usage
Failing to build my own Arabidopsis thaliana reference with Mt using 'codon.Mitochondrial' and Pt using 'codon.Bacterial\_and\_Plant\_Plastid' but the same commands works well with both lines with codon tables removed below.

The snpeff commands

```
java -Xms24g -Xmx24g -jar /opt/biotools/snpEff/snpEff.jar build -gff3 -v Arabidopsis_thaliana.TAIR10.46

```

returns:

```
java.lang.RuntimeException: Error parsing property 'Arabidopsis_thaliana.TAIR10.46.Pt.codonTable'. No such codon table 'codon.Bacterial_and_Plant_Plastid'
        at org.snpeff.snpEffect.Config.createCodonTables(Config.java:169)
        at org.snpeff.snpEffect.Config.readConfig(Config.java:650)
        at org.snpeff.snpEffect.Config.init(Config.java:480)
        at org.snpeff.snpEffect.Config.<init>(Config.java:117)
        at org.snpeff.SnpEff.loadConfig(SnpEff.java:451)
        at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:1000)
        at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:984)
        at org.snpeff.SnpEff.run(SnpEff.java:1183)
        at org.snpeff.SnpEff.main(SnpEff.java:162)

```

(same error for the Mt when removing Pt)

my fasta has the two chromosomes

```
>1 dna:chromosome chromosome:TAIR10:1:1:30427671:1 REF
>2 dna:chromosome chromosome:TAIR10:2:1:19698289:1 REF
>3 dna:chromosome chromosome:TAIR10:3:1:23459830:1 REF
>4 dna:chromosome chromosome:TAIR10:4:1:18585056:1 REF
>5 dna:chromosome chromosome:TAIR10:5:1:26975502:1 REF
>Mt dna:chromosome chromosome:TAIR10:Mt:1:366924:1 REF
>Pt dna:chromosome chromosome:TAIR10:Pt:1:154478:1 REF

```

my data folder is OK

```
data/Arabidopsis_thaliana.TAIR10.46/
drwxr-xr-x   3 u0002316 domain users 4.0K Apr 15 13:33 .
drwxr-xr-x 244 u0002316 domain users  16K Apr 10 09:39 ..
-rw-r--r--   1 u0002316 domain users 106M Apr 10 09:40 genes.gff
-rw-r--r--   1 u0002316 domain users 117M Apr 10 09:40 sequences.fa

```

my snpeff.config contains:

```
Arabidopsis_thaliana.TAIR10.46.genome : Arabidopsis_thaliana.TAIR10.46
Arabidopsis_thaliana.TAIR10.46.Mt.codonTable : codon.Mitochondrial
Arabidopsis_thaliana.TAIR10.46.Pt.codonTable : codon.Bacterial_and_Plant_Plastid

```

the two codon tables are present by default in the snpeff.config

```
codon.Mitochondrial                         : TTT/F, TTC/F, TTA/L+, TTG/L+, TCT/S, TCC/S, TCA/S, TCG/S, TAT/Y, TAC/Y, TAA/*, TAG/*, TGT/C, TGC/C, TGA/W, TGG/W, CTT/L, CTC/L, CTA/L, CTG/L+, CCT/P, CCC/P, CCA/P, CCG/P, CAT/H, CAC/H, CAA/Q, CAG/Q, CGT/R, CGC/R, CGA/R, CGG/R, ATT/I+, ATC/I+, ATA/I+, ATG/M+, ACT/T, ACC/T, ACA/T, ACG/T, AAT/N, AAC/N, AAA/K, AAG/K, AGT/S, AGC/S, AGA/R, AGG/R, GTT/V, GTC/V, GTA/V, GTG/V+, GCT/A, GCC/A, GCA/A, GCG/A, GAT/D, GAC/D, GAA/E, GAG/E, GGT/G, GGC/G, GGA/G, GGG/G

codon.Bacterial_and_Plant_Plastid           : TTT/F, TTC/F, TTA/L, TTG/L+, TCT/S, TCC/S, TCA/S, TCG/S, TAT/Y, TAC/Y, TAA/*, TAG/*, TGT/C, TGC/C, TGA/*, TGG/W, CTT/L, CTC/L, CTA/L, CTG/L+, CCT/P, CCC/P, CCA/P, CCG/P, CAT/H, CAC/H, CAA/Q, CAG/Q, CGT/R, CGC/R, CGA/R, CGG/R, ATT/I+, ATC/I+, ATA/I+, ATG/M+, ACT/T, ACC/T, ACA/T, ACG/T, AAT/N, AAC/N, AAA/K, AAG/K, AGT/S, AGC/S, AGA/R, AGG/R, GTT/V, GTC/V, GTA/V, GTG/V+, GCT/A, GCC/A, GCA/A, GCG/A, GAT/D, GAC/D, GAA/E, GAG/E, GGT/G, GGC/G, GGA/G, GGG/G

```

Anybody ever succeeded to build a custom genome with specific codon tables?

Thanks in advance

# Answer
**codon table issue**

I solved the issue with codon tables, the correct way to specify a table is the following:

```
Arabidopsis_thaliana.TAIR10.46.genome : Arabidopsis_thaliana.TAIR10.46
Arabidopsis_thaliana.TAIR10.46.Mt.codonTable : Mitochondrial
Arabidopsis_thaliana.TAIR10.46.Pt.codonTable : Bacterial_and_Plant_Plastid

```

*Note* that one should remove the leading prefix '**codon.**' from the name of the tables chosen from the config file

**Concerning accessory files:**

data to reproduce this issue:

```
baseurl="ftp://ftp.ensemblgenomes.org/pub/plants/release-46"
wget ${baseurl}/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
wget ${baseurl}/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
wget ${baseurl}/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
wget ${baseurl}/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.46.gff3.gz
wget ${baseurl}/variation/vcf/arabidopsis_thaliana/arabidopsis_thaliana_incl_consequences.vcf.gz

```

I suspect that the format used by Ensembl is not 100% OK for SnpEff; I therefore simplified the fasta headers to keep only the first argument (name) and remove the remaining of the line

eg.

```
>AT3G05780.1 pep chromosome:TAIR10:3:1714941:1719608:-1 gene:AT3G05780 transcript:AT3G05780.1 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:LON3 description:Lon protease homolog 3, mitochondrial [Source:UniProtKB/Swiss-Prot;Acc:Q9M9L8]

```

becomes

```
>AT3G05780.1

```

**Despite this modification, I still have an issue with**

```
java -jar /opt/biotools/snpEff/SnpEff.jar build -gff3 -v Arabidopsis_thaliana_TAIR10_46 2>&1 | tee -a log.txt

```
* a LOT of warning about not being able to add UTR while it seems that UTR's are well present in the GFF3 file

> Reading GFF3 data file : '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' WARNING: Cannot find transcript 'transcript:AT1G01010.1'. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 18 '1 araport11 five\_prime\_UTR 3631 3759 . + . Parent=transcript:AT1G01010.1' WARNING: Could not add UTR. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 18 '1 araport11 five\_prime\_UTR 3631 3759 . + . Parent=transcript:AT1G01010.1' WARNING: Cannot find transcript 'transcript:AT1G01020.3'. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 69 '1 araport11 three\_prime\_UTR 6788 6914 . - . Parent=transcript:AT1G01020.3' WARNING: Could not add UTR. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 69 '1 araport11 three\_prime\_UTR 6788 6914 . - . Parent=transcript:AT1G01020.3' WARNING: Cannot find transcript 'transcript:AT1G01020.5'. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 88 '1 araport11 three\_prime\_UTR 6788 6914 . - . Parent=transcript:AT1G01020.5' WARNING: Could not add UTR. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 88 '1 araport11 three\_prime\_UTR 6788 6914 . - . Parent=transcript:AT1G01020.5' WARNING: Cannot find transcript 'transcript:AT1G01020.4'. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 109 '1 araport11 three\_prime\_UTR 6788 6914 . - . Parent=transcript:AT1G01020.4' WARNING: Could not add UTR. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 109 '1 araport11 three\_prime\_UTR 6788 6914 . - . Parent=transcript:AT1G01020.4' WARNING: Cannot find transcript 'transcript:AT1G01020.1'. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 128 '1 araport11 three\_prime\_UTR 6788 6914 . - . Parent=transcript:AT1G01020.1' WARNING: Could not add UTR. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 128 '1 araport11 three\_prime\_UTR 6788 6914 . - . Parent=transcript:AT1G01020.1' WARNING: Cannot find transcript 'transcript:AT1G01030.1'. File '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/genes.gff' line 155 '1 araport11 three\_prime\_UTR 11649 11863 . - . Parent

* nothing good coming out of the cds.fa input although I treated it with the same name fix as the pep file

> 00:00:19 CDS check (FASTA file): '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/cds.fa'
>
> 00:00:19 Checking database using CDS sequences 00:00:19 Reading CDSs from file '/opt/biotools/snpEff\_v4.4/./data/Arabidopsis\_thaliana\_TAIR10\_46/cds.fa'... 00:00:20 done (75976 CDSs). 00:00:20 Comparing CDS... Labels: '+' : OK '.' : Missing '\*' : Error .................................................................................................... .................................................................................................... ....................................................................................................
>
> many more lines with dots .................................................................................................... .................................................................................................... .................................................................................................... ............................... CDS check: Arabidopsis\_thaliana\_TAIR10\_46 OK: 0 Warnings: 0 Not found: 104147 Errors: 0 Error percentage: NaN% Fatal error: No CDS checked. This is might be caused by differences in FASTA file transcript IDs respect to database's transcript's IDs. Transcript IDs from database (sample): 'transcript:AT3G12250.7' 'transcript:AT5G18570.1' 'AT3G27810.1' 'AT3G27810.2' 'transcript:AT2G16810.1' 'transcript:AT3G01170.1' 'AT4G37030.1' 'transcript:AT1G79940.4' 'transcript:AT2G44850.1' 'AT1G49540.1' 'AT1G49540.2' 'transcript:AT4G35290.5' 'transcript:AT4G15840.2' 'transcript:AT4G27500.1' 'AT4G37022.1' 'transcript:AT2G40008.1' 'AT3G27820.1' 'transcript:AT4G24265.1' 'transcript:AT2G20618.1' 'transcript:AT4G18375.2' 'AT4G37020.2' 'AT4G37020.1' Transcript IDs from database (fasta file): 'AT4G02541' 'AT1G12410.1' 'AT4G26510' 'AT4G02540' 'AT3G04500.1' 'AT3G04500.2' 'AT3G04500.3' 'AT2G36250' 'AT5G63350' 'AT2G36255' 'AT2G12280' 'AT1G03970.1' 'AT5G18390.1' 'AT3G04490.4' 'AT5G38386' 'AT3G04490.2' 'AT3G04490.3' 'AT3G04490.1' 'AT5G38380' 'AT4G02530' 'AT4G26520' 'AT1G69280.2'

**link to the full stdout log after running the command without cds.fa**

https://www.dropbox.com/s/gxa4e533p64wxm9/SnpEff\_build-log.txt?dl=1
> 2 votes

---
Tags: sequence-annotation, variants
---


# Question
Title: Chromosomal order of supplementary alignment in BAM file "SA" tag
I have several long-read sequences, which when aligned to a human reference genome, aligns to multiple chromosomal fragment as a result of chromothripsis shattering the linear DNA and randomly piecing them together. This creates for a single primary alignment multiple supplementary alignments.

With a single example of a read (only showing parts of all the bamfile information) being

`cdaf1fcb-27bf-4890-8ac7-a097dfa12c38 16 chr11 15151487 60`

With the primary alignment starting at chr11 at 15151487

Then the supplementary alignment can be seen as follows `SA:Z:chr13,49827969,+,1839S956M25D1584S,60,53;chr11,15152459,-,833S749M22D2797S,60,72;chrX,138324656,-,2540S404M2D1435S,60,20;chr11,15152066,+,278M11D4101S,60,20;`

Where the alignments are to chr13, 11, X and 11. With both regions on chr11 being close to the initial region from the primary alignment, making the regions on chromosome 11 the same (I know I haven't calculated the full chromosomal interval using the Cigar string for this example).

I assume the order of that supplementary information above is in the actual order. But the problem is that my DNA should be circular which doesn't fit with the order of the alignments since the circle would then be chr 11 (primary, start of circle) -\> chr 13 -\> chr11 (same as primary) -\> chrX -\> chr11 (back to the start of circle)

My question is: Q) So is the order of the supplementary information just pieces together randomly or is it true that I can assume the order is the true order of the sequencing read?

So thanks a lot for your help.

# Answer
> 0 votes
There is nothing in the SAMtags specification that dictates the order of the entries in an SA tag. It may be entirely dependent on the aligner that you used.

Looking at my own data that was aligned with bwa, I cannot see a definite pattern. Longer matches tend to be first, but not always.

Your aligner documentation may have more information on how it uses SA tags.

---
Tags: sequence-alignment, sam, long-reads
---


# Question
Title: Significance of k-mer length in COVID-19 sequence analysis?
I'm getting started in biology and bioinformatics with sequencing the SARS-Cov-2 or Coronavirus genome. I'm interested in this code which identifies k-mers in the genome:

```
from skbio import DNA
DNA(str(seq.seq)).kmer_frequencies(5)
kmers = DNA(str(seq1.seq)).kmer_frequencies(5)
kmers
{'ATTAA': 60, 'TTAAA': 95, 'TAAAG': 64, 'AAAGG': 45, 'AAGGT': 46, 'AGGTT': 51, 'GGTTT': 55,
'GTTTA': 69, 'TTTAT': 75,...}

```

The author of the code wants to use k-mer breakdowns to vectorize the data for machine learning purposes. What I'm curious about is why the k-mer in this case is five nucleotides? What factors should go into my selection of the length of the k-mer for analysis?

# Answer
> 0 votes
5bp produces a high resolution analysis, providing its sufficiently robust.

Its a sliding window of 5 bp wide and 1 bp step length which runs along the length of each genome. They are performing SNP analysis/classification. The original sequence they were examining was,

```
ATTAAAGGTTTAT

```

In any sliding window analysis the window size is arbitary towards the resolution and robustness of the analysis. The author believes that 5bp provides sufficiently robust information for a very high resolution classification/training analysis in their example.

The proof is the `accuracy` of the classification. I assume they have good results from 5bp. In ML the rationale is the general laws of ML (over/under training) and assessed on the test split in the data (probably 20%).

A key reason for the 1bp step is that this is the way ML works, you need a lot of data and by doing 5 bp sliding window at 1 bp step that will generate ~30 000 pieces of data for each sequence. So that's 3 million data points for 100 sequences.

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


# Question
Title: Database of RNA-Seq and SNP genotyping data
I need to obtain publicly available RNA-seq + SNP genotyping data for a cohort - ideally one where the RNA-seq data has been sequenced from brain tissue.

I have found CMC and GTEx to be sufficient for my needs, however the procedure to obtaining data from both these sources are extremely long and they request us to share our data with them as well; which we can not do because it is not ours to share.

I was wondering if anyone knows of another resource that can be used for obtaining such data?

Thanks

# Answer
> 1 votes
Human RNA-seq data, like genetic data is very sensitive data as it can be used to identify people (and their relatives), so people will likely not distribute it without entering a formal agreement that ensures that people can be held accountable if the data is misused, stolen or made public. Depending on the consent form signed by the patients (or their relatives), they may not be allowed to share it at all. One possibility would be to start a collaboration with a group that has access to such data and ask them to run your analyses on their data. You can then have a shared publication and it's a win for everybody.

---
Tags: rna-seq, database, snp, genotyping
---


title: Downloading SRA Files from AWS

text:
# Question
Title: Downloading SRA Files from AWS
I want to download the original BAM files that the authors had uploaded to SRA. Normally, I would just use `sam-dump`, but the files are having issues that seem related to this issue. Since according to the entry, AWS S3 also hosts the original BAM files, I thought I could download these directly.

## Example SRA Entry

NCBI documentation implies that I can't download this directly, but I can freely copy to other AWS locations within the region. To this end, I created my own S3 bucket (`mm-mneuron`) and am now trying to copy from the SRA bucket to mine. Here's what I try:

```
import boto3
import botocore

s3 = boto3.resource('s3')

bam_file = {
  'Bucket': 'sra-pub-src-6',
  'Key': 'SRR5253957/RPI25_0.bam'
}

my_bucket = s3.Bucket('mm-mneuron')

my_bucket.copy(bam_file, 'RPI25_0.bam')

```

This fails with:

> botocore.exceptions.ClientError: An error occurred (403) when calling the HeadObject operation: Forbidden

That is, it sounds like I can't access the SRA bucket. I've tested downloading and uploading to my bucket, so I know I have write permissions. Not sure what else to try here.

**How can I access the SRA data on S3?**

# Answer
A member of the SRA submission staff pointed out that using

```
prefetch --type all SRR5253957

```

will download the original files. In this case, it means running the above within an EC2 instance colocated with the S3 bucket (so, **us-east-1**) and having installed and configured SRA Toolkit to work from AWS (as per this documentation).

Unfortunately, the particular files I am concerned with are not currently accessible, but this should generically work in most situations.
> 2 votes

---
Tags: rna-seq, python, phylogenetics, aws
---


# Question
Title: How to deal with zeros in the denominator when calculating percent enrichment?
I have a `data.frame` with results from a CRISPR screen.

I want to calculate the % enrichment of a variable of interest `var` from baseline to treatment group.

Here is the distribution of `var` in linear:

I can do scale `var` with log transformation

But many rows with `var` = 0 will be lost. Moreover, I cannot compare % enrichment by dividing `treatmentgroup::var` to `Control/baseline::var` because I will run into the problem of dividing by zero.

I've seen people circumvent this problem by adding a small value ex)0.01 to the entire column of `var`. But I'm wondering if there's a better way of circumventing this problem. For example, a simple transformation that can deal with zeros and negative values.

Appreciate any input!

# Answer
> 1 votes
Cube root transformation is another way to deal with zeros without dropping them.

```
var^(1/3)

```

**EDIT** as Michael pointed out, square root transformation is another--perhaps more common--transformation that leaves zeros as zeros.

```
sqrt(var)

```

---
Tags: r, statistics
---


title: Java error when launching Pilon

text:
# Question
Title: Java error when launching Pilon
I am trying to use Pilon to improve a reference based on some Illumina data I have got. So, I aligned the Illumina reads to the reference using bwa. Then, I want to use Pilon to improve the reference using the alignment file in BAM format.

Information about Pilon is found here. I installed Pilon using conda.

The output of 'pilon' is:

```
Pilon version 1.22 Wed Mar 15 16:38:30 2017 -0400

    Usage: pilon --genome genome.fasta [--frags frags.bam] [--jumps jumps.bam] [--unpaired unpaired.bam]
                 [...other options...]
           pilon --help for option details

```

The output of:

```
pilon --genome ref.fa --frags readsSorted.bam --output pilo --outdir output/scaffold

```

is:

```
 Exception in thread "main" java.lang.reflect.InvocationTargetException
    at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
    at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
    at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
    at java.lang.reflect.Method.invoke(Method.java:498)
    at com.simontuffs.onejar.Boot.run(Boot.java:340)
    at com.simontuffs.onejar.Boot.main(Boot.java:166)
Caused by: java.lang.OutOfMemoryError: Java heap space
    at org.broadinstitute.pilon.GenomeRegion.<init>(GenomeRegion.scala:55)
    at org.broadinstitute.pilon.GenomeFile$$anonfun$contigRegions$1.apply(GenomeFile.scala:72)
 at org.broadinstitute.pilon.GenomeFile$$anonfun$contigRegions$1.apply(GenomeFile.scala:72)
    at scala.collection.TraversableLike$$anonfun$map$1.apply(TraversableLike.scala:234)
 at scala.collection.TraversableLike$$anonfun$map$1.apply(TraversableLike.scala:234)
    at scala.collection.immutable.Range.foreach(Range.scala:160)
    at scala.collection.TraversableLike$class.map(TraversableLike.scala:234)
 at scala.collection.AbstractTraversable.map(Traversable.scala:104)
 at org.broadinstitute.pilon.GenomeFile.contigRegions(GenomeFile.scala:72)
 at org.broadinstitute.pilon.GenomeFile$$anonfun$2.apply(GenomeFile.scala:52)
    at org.broadinstitute.pilon.GenomeFile$$anonfun$2.apply(GenomeFile.scala:52)
    at scala.collection.immutable.List.map(List.scala:277)
    at org.broadinstitute.pilon.GenomeFile.<init>(GenomeFile.scala:52)
    at org.broadinstitute.pilon.Pilon$.main(Pilon.scala:98)
    at org.broadinstitute.pilon.Pilon.main(Pilon.scala)
    ... 6 more

```

# Answer
Like gringer already commented, the cause of your error is nicely given by the `java` program. Your default `max heap` settings are probably too low for the program. If you have more memory on your machine you can give `max heap` a higher value. See this post and answers on Stack overflow.

You'll probably need to use the `-Xmxn` argument in `java`.
> 2 votes

# Answer
"If pilon was installed by condo, edit the file "/Users/user/anaconda/bin/pilon" by changing the maximum memory limit from 1G to higher numbers like "14G". It works. " check the thread: https://github.com/broadinstitute/pilon/issues/26
> 0 votes

---
Tags: assembly, reference-genome
---


title: Visualizing cell growth

text:
# Question
Title: Visualizing cell growth
I model the following events: a rod-shaped cell with center $x$ grows symmetrically until it reaches a maximum length/age, it divides into two identical cells, and the process continues. At division, the direction angle $\varphi$ of the daughter cells changes randomly.

So, the growth rate is $l=\frac{l\_{max}}{2}+\frac{a}{a\_{max}}\frac{l\_{max}}{2}$

Is there any way of naively visualizing this? I found different cell modellers on GitHub but they were very complicated biologically speaking. I have some basic Python knowledge, but I have never done simulations.

# Answer
> 1 votes
I would recommend looking into Python's plotting utilities, for example here.

I would recommend specifically looking at the examples that plot mathematical functions.

The basic idea is that you create a function that computes your equation based on the different parameters as inputs ($a, a\_{max}, l\_{max}$) and then you can plug in a range of values for each parameter and plot e.g. with $a\_{max}$ on the x axis, then different lines showing a range of values for $l\_{max}$, or whatever makes sense for you.

e.g.

```
import numpy as np
import matplotlib.pyplot as plt

def fn(a, amax, lmax):
  return lmax/2 + (a/amax) * (lmax/2)

plt.figure()
amax = np.arange(0.01, 1, .01)
plt.plot(amax, fn(.005, amax, .1), "ro") ## red points
plt.plot(amax, fn(.005, amax, .2), "go") ## green point
plt.show()

```

I have no idea if these parameter ranges are meaningful for you, you will have to figure that out for yourself.

---
Tags: python, modelling, bacteria
---


# Question
Title: Generate ligands candidates based on protein shape
Recent approaches to novel drug design using machine learning (ML) and deep learning, often involve generating hundreds of potential ligands which are later tested by docking with a target protein and recording the resulting binding affinity.

Is it possible to go the other way around by first constraining the ligand generation process to the shape of the protein thus limiting the search space, as opposed to generating ligands of random shapes?

# Answer
> 7 votes
TL;DR: docking is much slower than any ML approach, but the ML approach can be constrained by pharmacophores dictated by the active site.

## Side note: Scale

The scale for ligand space exploration is generally several orders of magnitude higher than "hundreds": Zinc DB lists 750 million enumerated compounds, GDB13 has enumerate all possible compounds up to 13 atoms (970 million), Enamine Real is 1.2 billion etc. In an actual experiment, one would use a more restrictive set. As an example of a real project, in Covid19 moonshot project there are 2,800 user suggested follow-ups.

## Side note: More than just binding

There are a few requirements for a drug to be effective, and binding is only one of them:

* Must bind...
* Must follow Lipinksi's rule of five, i.e. it must be membrane permeable and small.
* Must be synthesisable. Whereas the experience of someone learnèd is by far the best way to tell how an expansion can be done, there are several machine leaning approaches that can filter a huge dataset. IBM actually have a fun online tool —.
* Must not have side effects

So pre-filtering the list of compounds, regardless of docked poses, will happened, be it manually, by cut-offs or by advanced machine learning.

## CPU time: docking vs. machine learning

These various steps are generally followed in that order. Calculating properties (MW, logP etc.) for a given SMILES string to filter based on Lipinski's rule takes a fraction of a second, so that is generally always a first step.

Whether the compound can synthesised or bought is a general left last and done manually. Pharmacology is complicated so is ignored.

Docking a ligand takes time and different approaches have different time scales. A decent algorithm that uses implicit solvent model takes about a minute per core. A short MD simulation takes an 1 hour per core or more. There are docking programs for screening that can be faster but these dock only one conformer of the small molecule against a rigid protein and the results are not of much value. Also setting up docking properly is not easy (cf. steps required for Autodock 4, each with their own caveats). So to recap what needs to be considered that results in the slowness:

* Hundreds of conformers of the ligand
* Some flexibility of the protein, side-chains repacking only or even backbone changes.
* thousands of poses explored.
* Optionally. Implicit water is not as good as a TIP3 water model

Machine learning have a much smaller time requirement and is way more attractive for publications. Generally these work by using ligand only information, such as properties, splitting it into pharmacophores, colour etc. This requires a decent dataset, such as empirical results, say from a fragment screen, what compounds bound and which didn't. Often there is a lot of hype with machine learning and the result is something close to something already seen, as well put in this blog post by Pat Walters. This is a side effect of the starting data. The recent Stokes et al. 2020 Cell paper uses AI perfectly because of their staggering automated empirical data collection pipeline.

## ML after docking

However, ML can actually be used after docking to improve the quality of the score. The Boyles, Deane, Morris 2019 Bioinformatics paper investigated ligand-only machine learning to see why do ML methods that rely solely on ligand somehow often work. It is a good paper from a sceptical computational biochemistry viewpoint. The scorefunction used in docking by the sampler can be classified as forcefield-based, empirical or hybrid. In some cases empirical cases, many terms are ligand based properties (/features) that would normally go into a ML search (say one trained on failed hits and successful hits). For example Autodock Vina (not 4) has regression derived factors in the final weights of the score

## Structural features in machine learning: HotSpots

Information of the protein active site can be used in ranking compounds in a non-docking machine learning approach. Namely, finding out what kind of chemical groups bind preferentially to different parts of the protein's active site, *i.e.* making a "fragment hotspot map". Specifically, a hotspot map allows the determination of a pharmacophore, which is an abstract generalisation of a group of molecules with both both electronic and steric features —a coarse-grain model, if you will. This approach is also used to detect where the active site actually is in the extreme cases where it is not know (primarily a problem for pipelines, not humans). There are several implementations as described here, the first being FTMap (site, paper), which docks 16 different types of compounds. Waters can play a role and severely complicate matters, hence approaches can be taken to counter that. However, unlike docking protein rigidity is not a big nuisance: this is because a hotspot map is most often used for filtering/weighing by the "color" not the shape of the molecules. But shape based searching is possible.

### Shape based scoring

In terms of shape-based searching, this is actually a question of having pre-computed conformers, which can be problematic (but can and is often done on the fly), and shape based scoring —for which there are different metrics, such as the classical RMSD to more complexly weighted scores such as SuCOS.

---
Tags: protein-structure, machine-learning, docking
---


title: Protein misfolding

text:
# Question
Title: Protein misfolding
I am looking for literature on protein and RNA mis-folding. I am on the computational/biophysics side, and what interests me are the practical applications: diseases caused by misfolding, drug discovery, etc.

# Answer
Excellent question! I would recommend that you start with an NIH resource at the National Library of Medicine (NLM). PubMed (https://pubmed.ncbi.nlm.nih.gov/) provides search capabilities across the biomedical literature. You can also create an account and use it to manage your searches of the literature. An example search below for "protein misfolding disease" in PubMed retrieved almost 50,000 hits. You can filter to make the results more specific to your interests. You can view the citations, including the abstracts, and often links to the full text. Since you mention that you are on the "computational/biophysics side," you might be interested in the PubMed API (https://www.ncbi.nlm.nih.gov/pmc/tools/developers/) and python modules that use the API, such as Biopython.

E-utilities:

> https://www.ncbi.nlm.nih.gov/books/NBK25501/

Biopython:

> https://biopython.org/DIST/docs/api/Bio.Entrez-module.html

PubMed Search:

> https://pubmed.ncbi.nlm.nih.gov/?term=protein%20misfolding%20disease

Specific Articles:

> https://www.nature.com/scitable/topicpage/protein-misfolding-and-degenerative-diseases-14434929/ https://www.nature.com/articles/nrn1007
> 2 votes

---
Tags: proteins, protein-structure, rna, literature-search, rna-structure
---


# Question
Title: estimate genome size: kmer-based approach from PacBio reads
Can anyone suggest a software/method for kmer analysis using PacBio reads (RSII)?

Something similar to Jellyfish, that I saw in a nice tutorial \- but must be suitable for long, noisy reads. kmercounterexact from BBMapTools might also be an option, but again tricky with long, noisy reads.

The closest solution I could find was to estimate the genome size from overlap statistics (for the raw reads, the error corrected preassembled reads or by mapping the raw reads to the assembled contigs) by:

1. running an initial assembly using a small seed read length
2. then plot the pre-assembled read overlap histogram.

What do you think of this approach?

I would like to use this kmer analysis to estimate genome size, but for now I could not find much. I do not have companion Illumina reads - if so, I could try the approach proposed by Carvalho et al. 2016

I have reads & in-house genome assemblies (Canu) from different isolates of the same species of fungus to compare.

# Answer
> 5 votes
I don't think there is a method that would estimate a genome size using raw long reads.

The genome size estimates based on raw reads are done by fitting a model to kmer spectra (for instance Genomescope). The kmer spectra built from long reads are really messy due to the high error rate of long reads. That makes fitting of a model quite difficult. These methods assume small error rates. If the error rate is as high as 0.1 (which is normal for long reads), then the probability of having one kmer of 23 bases correct is ~8%. Meaning that if the coverage is 100x, just ~8x of homozygous kmers would be the real genomic kmers and it's very likely that kmer peaks of the sequencing errors, heterozygous kmers and homozygous kmers will just blend into one.

I can imagine, that if your coverage would crazy high (like 600x or more) the kmer spectra could separate errors and genomic kmers. You could also try to decrease kmer size to increase the proportion of correct kmers.

Another option would be to use a trick that they used to generate solid kmers for genome assembly using generalized de brujin graph using noisy reads. Maybe it would be better to take a look on the newer version of the assembler, that is called Flye.

But to be honest, I think for a long reads, a genome assembly is much more reliable way how to predict a genome size than any of the kmer tricks. I would simply assemble the genome, remap back reads, check the uniformity of coverage and conclude the genome size out it.

-- edit --

I tried to find some comprehensive evaluation of the assembly sizes compared to genome sizes, but I could not, so I just checked two examples, one big, one small and they both match quite well. *Ammopiptanthus nanus* was assembled to 823.74 Mb while kmer estimate was 889 Mb (ref). *Plasmodium falciparum* was assembled to 23.5 Mb while the original genome size estimate is 22.8 Mb. I know it's just two, but there is also a good reason to think that the genome size will be reasonable, feel free to add here more examples.

If you would be worried about contamination in the assembly, you can first run it though one of the methods to detect a contamination (like Blobology).

# Answer
> 1 votes
I have recently come across this submission from Wang et al. (2020), which is a method using the tool kmerfreq on error-corrected long reads (e.g. generated using Canu, FLAS, or LoRMA). It looks to be a promising approach.

---
Tags: genome, software-recommendation, assembly, k-mer, pacbio
---


# Question
Title: Grouping a set of genes in a ggplot object
I have a boolean matrix of mutated or non mutated genes across some samples

```
> head(muts)
                  TP53 CDKN2A TSHZ3 ARID1A EPHA3 NOTCH1 ACVR2A ARID1B ARID2
LP6008336-DNA_G01    1      1     1      0     0      0      0      0     0
LP6008460-DNA_D01    1      1     0      1     0      0      0      0     0
LP6005500-DNA_D03    1      1     0      0     0      0      0      0     0
LP6008334-DNA_B02    1      0     1      0     0      0      0      0     0
LP6008334-DNA_A04    1      0     1      0     0      0      0      0     0
LP6005334-DNA_H01    1      0     0      0     1      0      0      1     0

```

By this code

```
pw_rt = data.frame(
  gene = c('PIK3CA', 'EGFR', 'ERBB2', 'PTEN'),
  pw = c('RTK'),
  stringsAsFactors = F
)

pw_cc = data.frame(
  gene = c('CDKN2A', 'CCNE1'),
  pw = c('Cell Cycle'),
  stringsAsFactors = F
)

pw = rbind(pw_rt, pw_cc)

pw_col = RColorBrewer::brewer.pal(n = 3, name = 'Set1')[1:2]
names(pw_col) = unique(pw$pw)

require(dplyr)

lmuts = reshape2::melt(muts) %>% as_tibble()
colnames(lmuts) = c('sample', 'gene', 'value')

pw_n = pw$pw
names(pw_n) = pw$gene

Np = nrow(muts)

lmuts %>%
  mutate(gene = paste(gene), PW = pw_n[gene]) %>%
  group_by(gene, PW) %>%
  summarise(N = sum(value == 1)) %>%
  ungroup() %>%
  mutate(
    PW = ifelse(is.na(PW), "None", PW),
    N = paste0(N, ' (', round(N/Np * 100, 1), '%)')
    ) %>%
  ggplot(aes(x = PW, y = gene, fill = PW)) +
  geom_tile() +
  geom_text(aes(label = N)) +
  theme_light() +
  theme(legend.position = 'bottom') +
  guides(fill = guide_legend('Pathway')) +
  scale_fill_manual(values = c(pw_col, `None` = 'gainsboro')) +
  scale_x_discrete(limits = c(names(pw_col), 'None')) +
  labs(
    x = "",
    y = 'Gene',
    title = "Occurrence of mutations"
  )

```

I am able to get proportion of altered pathway by these genes like

However I want to highlight oncogenes and tumour suppressors within this plot

For instance TP53 is a tumour suppressor and CCNE1 is an oncogene

I don't know how to manipulate the above code to highlight or partition genes to tumour suppressor or oncogene

Can you help?

EDITED

I have add extra column to lmuts like

```
> head(lmuts)
             sample gene value gene_class
1 LP6005690-DNA_H02 TP53     1        TSG
2 LP2000333-DNA_A01 TP53     1        TSG
3 LP6005409-DNA_D03 TP53     1        TSG
4 LP6008141-DNA_H02 TP53     1        TSG
5 LP6008336-DNA_E02 TP53     1        TSG
6 LP6008269-DNA_B06 TP53     1        TSG
>

```

but

```
> lmuts %>%
+     mutate(gene = paste(gene), PW = pw_n[gene]) %>%
+     group_by(gene, PW) %>%
+     summarise(N = sum(value == 1)) %>%
+     ungroup() %>%
+     mutate(
+         PW = ifelse(is.na(PW), "None", PW),
+         N = paste0(N, ' (', round(N/Np * 100, 1), '%)')
+     ) %>%
+     ggplot(aes(x = PW, y = gene, fill = PW)) +
+     geom_tile() +
+     geom_text(aes(label = N)) +
+     theme_light() +
+     theme(legend.position = 'bottom') +
+     guides(fill = guide_legend('Pathway')) +
+     scale_fill_manual(values = c(pw_col, `None` = 'gainsboro')) +
+     scale_x_discrete(limits = c(names(pw_col), 'None')) +
+     labs(
+         x = "",
+         y = 'Gene',
+         title = "Occurrence of mutations"
+     )+facet_wrap(~gene_class, ncol=1)
Error: At least one layer must contain all faceting variables: `gene_class`.
* Plot is missing `gene_class`
* Layer 1 is missing `gene_class`
* Layer 2 is missing `gene_class`
Run `rlang::last_error()` to see where the error occurred.

```

# Answer
If you had the classification into oncogene/TSG for each gene in a column (let's call it `gene_class`) in `lmuts`, you could split both groups into facets within the plot by adding:

```
+ facet_wrap(~gene_class, ncol=1)

```
> 3 votes

# Answer
I can't reproduce your graph, because you don't make `muts` available, but it appears your problem is having `x = PW` in your ggplot call. Don't you want to that to be whatever your variable name for sample is (ie, the rownames in your `but` dataframe you have at the beginning).
> 0 votes

---
Tags: r, ggplot2
---


# Question
Title: What does it mean when PDB lists multiple organisms for the same structure?
I'm browsing through ribosomes which have been modeled with Cryoem on PDB and am pretty confused by the fact that for some, PDB, or rather the authors of the deposition probably, put multiple organisms from different domains of life... How does one interpret that?

So far as i understand *bos taurus* is a cow and *e.coli* is a bacterium. So eukarya and bacteria.. Thank you.

# Answer
Since this is a 40S ribosomal protein, and therefore eukaryotic-specific (as the description says) the listing of *E. coli* is not intuitive as you note. If you go down to "Macromolecules" and click on "Nucleic Acids/Hybrids" you can see the listed *E. coli* sequence is just 28bp, and all `U`s (poly(U) tail?). So, you can safely ignore the bacterial annotation I think.
> 3 votes

---
Tags: pdb
---


# Question
Title: How to calculate average minor allele frequency (MAF) difference between two populations (i.e. European, Hispanic)
The data is present in the PLINK format .bed (individual genotypes), .bim (genetic markers) and .fam (sample IDs and disease phenotype file). I have also separate text files containing information on subjects related to two populations. (i.e. European and Hispanic)

# Answer
> 1 votes
1. You need to report allele frequencies using `plink` for some European samples
2. Do the same for Hispanic samples
3. Use a software to subtract the two allele frequencies obtained above and obtain the result.

---
Tags: bioconductor, plink
---


title: No MQ tags in VCF files

text:
# Question
Title: No MQ tags in VCF files
To call minority variants in my Mtb sequences I'm using a pipeline of

```
fastq ->
    BWA MEM for mapping ->
    picard MarkDuplicates ->
    GATK3 for indel realignment and base recalibration ->
    lofreq for variant calling (https://csb5.github.io/lofreq/)

```

However the resulting VCF files don't have an MQ tag for me to filter on. Where in this process would mapping quality information be generated or lost so that it doesn't end up in the VCF file?

# Answer
> 2 votes
As @ATpoint explains in his comment, the mapping quality (`MQ`) applies to read mappings and, hence, is contained in BAM files, but not VCFs. The `MQ` values carry important information that ideally the variant caller should use when determining the variant quality (`QUAL` in the VCF) of individual variant loci.

In the LoFreq documentation it says that

> # Choice of mapper: Bowtie and BWA-SW
>
> BWA-SW assigns very low mapping qualities to mapped reads. **LoFreq incorporates mapping quality by default into its variant calling model**, which doesn’t help in this scenario and will reduce its sensitivity. Consider disabling the use of mapping quality for `lofreq call` with `-N`.

So LoFreq is already aware of the `MQ` values and takes them into consideration during variant calling.

Since you are using BWA-MEM (and not BWA-SW) you should probably not be concerned about the issue they mention and trust their quality filtering.

If you are sceptic about the quality of mappings, though, you could still filter them out directly after mapping, e.g. using `samtools view -q 20 ...`

---
Tags: ngs, variant-calling, bwa, gatk, bowtie2
---


# Question
Title: Programmatically retrieve graph features from Cytoscape
Is there a way to programmatically retrieve graph features (a networkx graph) from Cytoscape in Python 3?

EDIT: To add more details, I had a dictionary of proteins, each protein was associated with a list of proteins that it interacts with. I wanted to create a protein-protein interaction graph so that I could extract graph features from it, including things like k-core value of each node.

# Answer
I was able to extract the graph features that I wanted just using networkx, without going through the intermediate step of Cytoscape.

EDIT: Initially, I planned to import the data into Cytoscape and extract graph features from Cytoscape. However, I noticed that a much easier solution would be to just use the methods of the networkx library. I initially had a few graph features in mind that I wanted to extract like k-core value of each node, and I was able to do this using methods from networkX. To calculate k-core for instance, I used the k\_core method. For more information on networkX methods, check the docs
> 1 votes

---
Tags: python
---


# Question
Title: What's the meaning of having different alignments with the same qname in a bam with no duplicates?
I have removed duplicates from my bam file (not secondary, nor suplementary alignments), and then I have looked for some repeated `qnames` (to assess the failure of a program) and I've found a lot of them. All of them has `TLEN` equal to 0.

As an example: this is part of a bed file that a generate from the bam.

```
CHR  Start   End       QNAME                                 PNEXT  TLEN
12 86945793 86945944 A05325:240:SS348DSZX:2:5105:12246:3790 90044835 0
4 168461445 168461485 A05325:240:SS348DSZX:2:5105:12246:3790 90044835 0

```

And those are these reads in the original bam:

```
A05325:240:SS348DSZX:2:5105:12246:3790  2177    4       168461445       0       111H40M 11      90044835        0       ATAAAATTAAAAAAATAAAAAAATTAAAAAAAAAAAAAAA        :,F,F,,FF:FF,F,,F,FFFF,,F,F,:,:,FFFFFFFF     NM:i:2  MD:Z:28T8G2     MC:Z:151M       AS:i:32 XS:i:31 RG:Z:HM327DSXX.3        SA:Z:12,86945793,+,20S37M94S,1,1;
A05325:240:SS348DSZX:2:5105:12246:3790  129     12      86945793        1       20S37M94S       11      90044835        0       ATGGTGAAAAAACAAATTATTTAAAAAAATCAAAAAAAAAAAAAAAACCAAAAAACATAAAAAAAAGGGCAGTAATTTTATTCATGGGGTCAAACAAAAAAAAAGCATAATATAAAATTAAAAAAATAAAAAAATTAAAAAAAAAAAAAAA      FF,,,:,F,F,,,,F,,,,,,,F,FFFF:,:FFFF:FFF,FF::F,:,,FFFF,,,F,FF::,,F,::,,,,F,,F:FF,,,,,,,F::F,F,,,,F,FFFFF:,F,,F,,:,F,F,,FF:FF,F,,F,FFFF,,F,F,:,:,FFFFFFFF      NM:i:1  MD:Z:9G27       MC:Z:151M       AS:i:32 XS:i:29 RG:Z:HM327DSXX.3        SA:Z:4,168461445,+,111S40M,0,2;

```

This are the commands used to remove duplicates:

```
samtools view -F 1024 -bo filtered.bam original.bam

```

The commands used to generate this bed were:

```
samtools view sample.bam | awk '{ print $3, $4, $4+length($10), $1, $8, $9}' > sample.bed

```

Why do I present those alignments?

# Answer
I think that you are dealing with two different phenomena:

**1) TLEN == 0 reads**

From the BAM spec:

> TLEN: signed observed Template LENgth. If all segments are mapped to the same reference sequence, the absolute value of TLEN equals the distance between the mapped end of the template and the mapped start of the template, inclusively (i.e., end − start + 1).14 Note that mapped base is defined to be one that aligns to the reference as described by CIGAR, hence excludes soft-clipped bases. The TLEN field is positive for the leftmost segment of the template, negative for the rightmost, and the sign for any middle segment is undefined. If segments cover the same coordinates then the choice of which is leftmost and rightmost is arbitrary, but the two ends must still have differing signs. **It is set as 0 for a single-segment template or when the information is unavailable (e.g., when the first or last segment of a multi-segment template is unmapped or when the two are mapped to different reference sequences).**

Note that in each case the mate read appears to be mapped to a different chromosome/contig. So, the mapper can't estimate insert size in these cases, so it just says zero.

**2. Multiple reads with same qname**

There are various reasons why the same qname can appear multiple times. Again from the BAM spec, one such explanation is that it is a chimeric alignment:

> **Chimeric alignment:** An alignment of a read that cannot be represented as a linear alignment. A chimeric alignment is represented as a set of linear alignments that do not have large overlaps. Typically, one of the linear alignments in a chimeric alignment is considered the “representative” alignment, and the others are called “supplementary” and are distinguished by the supplementary alignment flag. All the SAM records in a chimeric alignment have the same QNAME and the same values for 0x40 and 0x80 flags (see Section 1.4). The decision regarding which linear alignment is representative is arbitrary.

You may have to do some digging and looking at the different flags to decide what exactly is going on, but my bet is that they are supplementary alignments.
> 2 votes

---
Tags: sequence-alignment
---


title: Available phenotype data from GTEx

text:
# Question
Title: Available phenotype data from GTEx
I am using gene expression data from the The Genotype-Tissue Expression (GTEx) project, available here: http://gtexportal.org. I would also like to use phenotype data for each sample (especially age, sex, race, body mass index, smoking history).

Access to such information is controlled by the NIH, so you need to apply for access.

Before I go through that probably lengthy process, I would like to know if it is worth it. **Are any of the phenotype data age, sex, race, body mass index, smoking history available on a per-sample level, if you have access?** I could not find this information anywhere.

# Answer
> 3 votes
If helpful for anyone else, I found that all and much more are available. See below 2 links.

https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/dataset.cgi?study\_id=phs000424.v8.p2&phv=169091&phd=3910&pha=&pht=2742&phvf=&phdf=&phaf=&phtf=&dssp=1&consent=&temp=1

https://ftp.ncbi.nlm.nih.gov/dbgap/studies/phs000424/phs000424.v8.p2/pheno\_variable\_summaries/phs000424.v8.pht002742.v8.p2.GTEx\_Subject\_Phenotypes.var\_report.xml

---
Tags: gene, public-databases, gene-expression, public-dataset, gtex
---


# Question
Title: How to run a one-time 'setup' step for amrfinder in snakemake using conda
I'm trying to add amrfinder running within conda to a snakemake workflow, working on a slurm HPC.

The first time you run `amrfinder`, you have to run `amrfinder -u` to update/download the database of genes/proteins which it searches for. Since I'm using the `--use-conda` flag of `snakemake`, the first time I run the workflow, it creates a new env for the `conda` package, so this command has to be run before that copy of `amrfinder` will work.

My question is: how do I run `amrfinder -u` within `snakemake`? If I do it as part of the standard `amrfinder` command, it runs as many times as there are samples, and as all the jobs are looking at the same database, when it updates for one job, it breaks the other jobs.

I can't figure out how to run this command in a 'snakemake-y' way, as it has no input, it's just literally `amrfinder -u`.

Of course, I can make it work with a hack/workaround, but I'm interested in a snakemake way of doing it.

This is the main snakemake rule for running amr\_finder\_plus, the new rule, to run `amrfinder -u` needs to run before this

```
rule amr_finder_plus:
    input:
        assembly = 'my_assembly.fasta'
    output:
        amr_finder_plus_results = 'my_amr_finder_plus_results.tsv'
    conda:
        '../../envs/amrfinderplus.yaml'
    shell:
        'amrfinder -n {input.assembly} -O Salmonella --output {output.amr_finder_plus_results} --threads 4'

```

# Answer
There are two approaches you can take:

1. There's a `--conda-create-envs-only` option that you can use to exit after creating the envs. You should get information about where the envs are so you can run `amrfinder -u`.
2. Add a new rule called "amr\_setup" or something like that that runs:

.

```
amrfinder -u
touch amrfinder.setup

```

Then put `amrfinder.setup` into your snakemake workflow as an input file into the amr\_finder\_plus rule. Either of these will work
> 1 votes

---
Tags: snakemake, conda, hpc
---


# Question
Title: How to concatename file names in Snakemake?
I want to pass a series of files (e.g., `a`, `b`) to a program `prog` like: `prog a,b`.

Is there a way to do this in Snakemake?

My attempt is:

```
rule run_prog:
  input: "a", "b"
  output: "output"
  run: "prog {input}"

```

But the actual command Snakemake generated is `prog a b` instead of `prog a,b`.

Any solution to pass a comma concatenated file names to Snakemake?

# Answer
> 1 votes
You can use join to convert the list of input file into a string. I think this would work as expected.:

```
rule run_prog:
    input: "a", "b"
    output: "output"
    run: f"prog {','.join(input)}"

```

---
Tags: snakemake
---


# Question
Title: Adding certain barcodes to metadata in Seurat
I want to upload an excel file sheet that has certain barcodes that I would like to show on my umap. How do I go about adding the file and linking it to the metadata?

Below is my following code.

```
# set up the working directory
wd = "/home/PTX_AAC656
setwd(wd)

# load counts
# SD-52-JZ33 patient sample but aligned to both mouse and human
SD_52_JZ33.data <- Read10X(data.dir = "../PTX_AAC656/metadata/patient_sample/SD-52-JZ33/outs/filtered_feature_bc_matrix/")
BNM_human_only.data <- Read10X(data.dir = "../PTX_AAC656/metadata/human_only_alignment/SD-52-JZ33_BNM_5K/outs/filtered_feature_bc_matrix/")
BNL_human_only.data <- Read10X(data.dir = "../PTX_AAC656/metadata/human_only_alignment/SD-52-JZ33_BNL_5K/outs/filtered_feature_bc_matrix/")
BNL_mouse_only.data <- Read10X(data.dir = "../PTX_AAC656/metadata/mouse_only_alignment/SD-52-JZ33_BNL_5K/outs/filtered_feature_bc_matrix/")
BNM_mouse_only.data <- Read10X(data.dir = "../PTX_AAC656/metadata/mouse_only_alignment/SD-52-JZ33_BNM_5K/outs/filtered_feature_bc_matrix/")

   BNL_both.data <- Read10X(data.dir = "../AAC656/SD-52-JZ33_BNL_5K/outs/filtered_feature_bc_matrix/")
    BNM_both.data <- Read10X(data.dir = "../AAC656/SD-52-JZ33_BNM_5K/outs/filtered_feature_bc_matrix/")

setwd(wd)

# create the Seurat objects
SD_52_JZ33 <- CreateSeuratObject(counts = SD_52_JZ33.data, min.cells = 0, project = "human")
BNL_human_only <- CreateSeuratObject(counts = BNL_human_only.data, min.cells = 0, project = "human_alignment_ptx")
BNM_human_only <- CreateSeuratObject(counts = BNM_human_only.data, min.cells = 0, project = "human_alignment_ptx")
BNM_both_ptx <- CreateSeuratObject(counts = BNM_both.data, min.cells = 0, project = "human_mouse_alignment_ptx")
BNL_both_ptx <- CreateSeuratObject(counts = BNL_both.data, min.cells = 0, project = "human_mouse_alignment_ptx")
BNL_both_ptx$model = "mouse_human_alignment_ptx"
BNM_both_ptx$model = "mouse_human_alignment_ptx"
SD_52_JZ33$model = "patient_aligned_mouse_human"
BNM_human_only$model = "human_alignment_only_ptx"
BNL_human_only$model = "human_alignment_only_ptx"
BNL_both_ptx$model = "mouse_human_alignment_ptx"
BNM_both_ptx$model = "mouse_human_alignment_ptx"
BNM_mouse_only <- CreateSeuratObject(counts = BNM_mouse_only.data, min.cells = 0, project = "mouse_alignment_ptx")
BNL_mouse_only <- CreateSeuratObject(counts = BNL_mouse_only.data, min.cells = 0, project = "mouse_alignment_ptx")
BNL_mouse_only$model = "mouse_alignment_only_ptx"
BNM_mouse_only$model = "mouse_alignment_only_ptx"

# merge samples
mouse_human_patient <- merge(BNL_mouse_only, y = c(BNM_mouse_only, SD_52_JZ33), add.cell.ids = c("A", "B", "C"), project = "HTB2876")
human_ptx_human_patient <- merge(BNL_human_only, y = c(BNM_human_only, SD_52_JZ33), add.cell.ids = c("D", "E", "F"), project = "HTB2876")
ptx_human_patient.1 <- merge(BNL_both_ptx, y = c(BNM_both_ptx, SD_52_JZ33), add.cell.ids = c("G", "H", "I"), project = "HTB2876")

# mouse ptx model with patient
# mark the mito genes
mito.genes <- grep(pattern = "^MT-", x = rownames(x = mouse_human_patient), value = TRUE)
mouse_human_patient[["percent.mt"]] <- PercentageFeatureSet(mouse_human_patient, pattern = "^MT-")
mouse_human_patient[["log_nCount_RNA"]] <- log2(mouse_human_patient[["nCount_RNA"]]+1)

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
mouse_human_patient <- subset(mouse_human_patient, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 30)
mouse_human_patient <- NormalizeData(mouse_human_patient, normalization.method = "LogNormalize", scale.factor = 10000)
mouse_human_patient <- FindVariableFeatures(mouse_human_patient, selection.method = "vst", nfeatures = 2000)
all.genes.mouse <- rownames(mouse_human_patient)
mouse_human_patient <- ScaleData(mouse_human_patient, features = VariableFeatures(object = mouse_human_patient), vars.to.regress = c("nCount_RNA"))
mouse_human_patient <- RunPCA(mouse_human_patient, features = VariableFeatures(object = mouse_human_patient))
mouse_human_patient <- FindNeighbors(mouse_human_patient, dims = 1:20)
mouse_human_patient <- FindClusters(mouse_human_patient, resolution = 0.4)
mouse_human_patient <- RunUMAP(mouse_human_patient, dims = 1:20)
DimPlot(mouse_human_patient, reduction = "umap", group.by = "model")

# human only alignment pix with patient
# mark the mito genes
mito.genes <- grep(pattern = "^MT-", x = rownames(x = human_ptx_human_patient), value = TRUE)
human_ptx_human_patient[["percent.mt"]] <- PercentageFeatureSet(human_ptx_human_patient, pattern = "^MT-")
human_ptx_human_patient[["log_nCount_RNA"]] <- log2(human_ptx_human_patient[["nCount_RNA"]]+1)
human_ptx_human_patient <- subset(human_ptx_human_patient, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 30)

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
human_ptx_human_patient <- NormalizeData(human_ptx_human_patient, normalization.method = "LogNormalize", scale.factor = 10000)
human_ptx_human_patient <- FindVariableFeatures(human_ptx_human_patient, selection.method = "vst", nfeatures = 2000)
all.genes.human.ptx <- rownames(human_ptx_human_patient)
human_ptx_human_patient <- ScaleData(human_ptx_human_patient, features = VariableFeatures(object = human_ptx_human_patient), vars.to.regress = c("nCount_RNA"))
human_ptx_human_patient <- RunPCA(human_ptx_human_patient, features = VariableFeatures(object = human_ptx_human_patient))
DimPlot(human_ptx_human_patient, reduction = "pca")
human_ptx_human_patient <- FindNeighbors(human_ptx_human_patient, dims = 1:20)
human_ptx_human_patient <- FindClusters(human_ptx_human_patient, resolution = 0.4)
human_ptx_human_patient <- RunUMAP(human_ptx_human_patient, dims = 1:20)
DimPlot(human_ptx_human_patient, reduction = "umap", group.by = "model")

# ptx and human
# mark the mito genes
mito.genes.ptx_patient <- grep(pattern = "^MT-", x = rownames(x = ptx_human_patient.1), value = TRUE)
ptx_human_patient.1[["percent.mt"]] <- PercentageFeatureSet(ptx_human_patient.1, pattern = "^MT-")
ptx_human_patient.1[["log_nCount_RNA"]] <- log2(ptx_human_patient.1[["nCount_RNA"]]+1)

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
ptx_human_patient.1 <- subset(ptx_human_patient.1, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 30)
ptx_human_patient.1 <- NormalizeData(ptx_human_patient.1, normalization.method = "LogNormalize", scale.factor = 10000)
ptx_human_patient.1 <- FindVariableFeatures(ptx_human_patient.1, selection.method = "vst", nfeatures = 2000)
all.genes.both.patient.1 <- rownames(ptx_human_patient.1)
ptx_human_patient.1 <- ScaleData(ptx_human_patient.1, features = VariableFeatures(object = ptx_human_patient.1), vars.to.regress = c("nCount_RNA"))
ptx_human_patient.1 <- RunPCA(ptx_human_patient.1, features = VariableFeatures(object = ptx_human_patient.1))
ptx_human_patient.1 <- FindNeighbors(ptx_human_patient.1, dims = 1:20)
ptx_human_patient.1 <- FindClusters(ptx_human_patient.1, resolution = 0.4)
ptx_human_patient.1 <- RunUMAP(ptx_human_patient.1, dims = 1:20)

#run plots
DimPlot(ptx_human_patient.1, reduction = "umap", group.by = "model")
DimPlot(ptx_human_patient.1, reduction = "umap", label = "TRUE")
saveRDS(ptx_human_patient.1, file = "derived_data/PTX_human_2876_res04_6000g.rds")
DimPlot(ptx_human_patient.1, reduction = "umap", group.by = "model")
DimPlot(ptx_human_patient.1, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
ptx_human_patient.1$clustering.idents <- Idents(object = ptx_human_patient.1)
Idents(ptx_human_patient.1) <- "clustering.idents"
new.cluster.ids <- c("Human PTX Tumor", "PTX Fibroblast", "Human Fibroblast", "PTX Macrophage", "Human PTX Tumor", "Human Mitochondrial", "Human Macrophage", "Human Endothelial", "Human T-cells", "Human PTX Tumor","Human Fibroblast", "PTX Fibroblast", "Human Tumor Ductal", "Human Acinar", "PTX Macrophage", "PTX Stellate", "PTX Endothelial")
names(new.cluster.ids) <- levels(ptx_human_patient.1)
ptx_human_patient.1 <- RenameIdents(ptx_human_patient.1, new.cluster.ids)
saveRDS(ptx_human_patient.1, file = "derived_data/PTX_human_2876_res04_6000g.rds")
p1 = DimPlot(ptx_human_patient.1, reduction = "umap", label = TRUE, pt.size = 0.5,  group.by = "clustering.idents") + NoLegend()
p2 = DimPlot(ptx_human_patient.1, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
print(CombinePlots(plots = list(p1, p2)))
markers2 <- FindAllMarkers(ptx_human_patient.1, min.pct = 0.25, logfc.threshold = 0.25)
markers2 %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
saveRDS(markers, file="derived_data/markers_HTB2876_PC20_res04_6000g.rds")

```

Part of my excel file

```
barcode muTaTionsTaTus
GAGATGGCACTAGGCC-1  G12V
CTGAGCGAGACTCAAA-1  G12V
TACCGGGGTCCGAAGA-1  G12V
TATACCTCATAGGCGA-1  G12V
CTGCCATGTATTGGCT-1  G12V
AGTGTTGTCAGGAGAC-1  G12V
AATGCCAAGCAGTACG-1  G12V
TTTGTTGAGGCTTCCG-1  G12V
ACAGAAATCCATTCGC-1  G12V
GTCATTTCAACTAGAA-1  G12V
GAGTGTTGTGAGATTA-1  G12V
TCTTGCGCATAACGGG-1  G12V
ACATCGATCGGCATCG-1  G12V
CATCCCACACACAGCC-1  G12V
GACAGCCGTAGCTGTT-1  G12V
CACAGATGTTAGGCCC-1  G12V
AACGGGATCACAACCA-1  G12V

```

# Answer
If you want to add metadata to a Seurat object the new metadata needs to have the same length as the old metadata (== number of cells).

For an example, if you have a 8 cells `c("a", "b", "c", "d", "e", "f", "g")` and you want to be able to highlight the 4 cells with the "barcodes" `c("a", "c", "f", "g")` in a umap plot you have to add a new grouping column:

```
# The Seurat object is called "ExampleData" and columns can be directly adressed using "$"
# The barcodes of all cells are stored as the rownames of the meta.data slot
ExampleData$cells_of_interest <- ifelse(rownames(ExampleData@meta.data) %in% c("a", "c", "f", "g"), "interesting", "boring")

```

From this you will get a new column in the meta data with the values:

```
interesing, boring, interesting, boring, boring, interesting, interesting

```

This way you can use this column in `DimPlot`:

```
DimPlot(ExampleData, reduction = "umap", group.by = "cells_of_interest")

```

With two colors in the plot one called "interesting" and one "boring"
> 2 votes

# Answer
From my understanding, you are just trying to add some new information to the metadata of your Seurat object. To do so you can just add the column to `meta.data` slots in your Seurat object.

(https://github.com/satijalab/seurat/wiki/Seurat#slots)

Something like this:

```
myBarcode = rownames(seurat_obj@meta.data) #get barcode from seurat
mutation.df = mutation.df[match(myBarcode, mutation.df$barcode), ] #match the order of seurat barcode with you data
seurat_obj$mutation = mutation.df$muTaTionsTaTus #put mutation into metadata

```
> 0 votes

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


title: Economist article on coronavirus

text:
# Question
Title: Economist article on coronavirus
I am wondering about an article in the Economist here: https://www.economist.com/briefing/2020/02/29/covid-19-is-now-in-50-countries-and-things-will-get-worse There is a graph there

The explanation is as follows The course of an epidemic is shaped by a variable called the reproductive rate, or R. It represents, in effect, the number of further cases each new case will give rise to. If R is high, the number of newly infected people climbs quickly to a peak before, for want of new people to infect, starting to fall back again (see chart 2). If R is low the curve rises and falls more slowly, never reaching the same heights. With sars-cov-2 now spread around the world, the aim of public-health policy, whether at the city, national or global scale, is to flatten the curve, spreading the infections out over time.

I am having trouble with this. There are 2 curves here. First of all, I am not sure of the different processes that determine them. It says one is without measures. The other is with measures. But it also says one is with high R and the other with low R. Is R a fixed constant for a disease? Alternatively, is there some variability in R and does that degree of variation affect transmission? This is not explained or stated in the article and I haven't heard R described in this way. At best it is not written clearly in my opinion.

In addition, does without measures necessarily imply high R? Is that what they are getting at? Also, it looks to me as though the blue graph goes to zero very quickly, basically when t is equal to the mean of the orange curve. That means that without measures the epidemic will end quickly (because it exhausts all the possible patients).

Finally the areas under the two curves is about the same. This means that the number of patients is about the same with or without measures.

Anyway, can anyone clear up any of this?

# Answer
R is a reference to R0, otherwise known as the basic reproductive rate, and means the number of new cases from a single patient infection. R0 has to be above 1 for a disease to persist. It is a classic equation which account for the probability of transmission against immunity. The y-axis is the number of infecteds, very poorly presented.

The idea is that R0 will be lower if there is containment (less people to infect because you can't move around), but will last longer, possibly because the number of naive patients (not encountered the virus) will be larger over a longer period of time.

That particular graph appears to be saying that everyone at risk of infection will get infected with or without control measures (the surface area under the graphs is approximately the same\*) and the only thing that will halt the virus is acquired immunity. If this is mainstream economic thinking then I beg to differ, albeit I note the source is CDC (again see note\*).

For example, all other respiratory viral infections are seasonal. If COVID-19 is seasonal containment is the right choice, because the window of transmission is limited.

\*, Given this is of CDC origin I think either its a mistake, or at the very least was taken out of context, because the idea is that overall there is a lower number of people infected as a result of "social distancing" and that does not appear to be what the two curves show.
> 5 votes

# Answer
R, the reproductive number, relates to the average number of (new) people that will get sick (infected) per person that is already sick.

For instance, if $R=2$ and you start with a single infected person, then the next generation will be 2 people, those 2 people will make 4 people sick, those 4 people will make 8 people sick, and so on.

---

> Is R a fixed constant for a disease?

This reproduction number, how many people the disease is passed on to, depends a lot on the circumstances. For instance, **how frequent** are contacts and how likely is **transmission per contact**. The probability of transmission per contact will depend on the type of disease (and from that stems those comparisons of different disease), but can also be influenced by measures like reducing frequency and intensity of contacts (social distancing).

The coronavirus spreads in droplets from breath, coughs and sneezes, and will reach other people mostly during close contact (Although there may be events of long-distance spreading if the virus is spread in aerosols, which is claimed to be not the case but it may still exist to a small extend, like with SARS in 2003 explained here and here). So, the reproduction number, R, will likely be influenced by measures like preventing close contact, what they also call social distancing, as well as proper hygiene (washing hands to prevent transmission from virus picked up from fomites). A clear example is the case of a sexually transmittable disease, where condom use, or not having sex at all, will reduce the spreading of the disease.

> Alternatively, is there some variability in R and does that degree of variation affect transmission?

There is a *lot* of variability in transmission rate.

The spread of the covid-19 virus has likely been accelerated by celebration of New Year's in China and skiing holidays and carnaval in Europe (as well as large scale meetings during soccer matches and church services). (this is why large events are bad)

Also, the disease will at some point satiate the surrounding. Once people have been infected they are likely immune for some while and less likely to be infected again. So once the entire household, classroom, or workplace is infected then new infections will be less likely. (this is why long distance travelling is bad because it brings the virus quickly into untouched area's)

It is wrong to speak about *the* reproduction number. The reproduction occurs heterogeneously and it is wrong to speak of a single number.

The average reproduction might be 1 or less, but there can *still* be effectively exponential growth within a sub-population with relatively higher transmission rates.

Especially when this sub-population is strongly connected network on a large scale. (the virus is like hopping from cluster to cluster in a large network, the transmission occurs mostly in dense clusters of the network and some area's remain relatively little affected, this is why Sweden does not need to take though measures, their network is relatively sparse and safe)

---

> Finally the areas under the two curves is about the same. This means that the number of patients is about the same with or without measures.

That sketch in the economist article is just a simplistic sketch. The 'flattening the curve' will reduce the rate of growth, but it should *also* reduce the total number of cases.

An example is the image below from this SE post:

https://stats.stackexchange.com/questions/461319/using-different-nodes-in-a-networked-compartmental-model-sir-for-different-reg

In the image you see the effect of different social distancing measures (more details about it in that post). On the left the cumulative cases are plotted (this relates to the area under the curves) and you see a reduction of the total cases. The curves are not *only* flattened, but also reduced in volume. (Of course, this does require that the measures are continuously followed for the entire period. That image is an idealistic situation)

> Also, it looks to me as though the blue graph goes to zero very quickly

The endemic will eventually halt when a lot of people have already gotten sick (it may also continue for a long time, and become endemic, because the population is not static, and births introduce new people that have not been sick before, also immunity may become lost over time).

The *basic* reproduction value, often denoted $R\_0$, is the number of people that the infection is passed on to *when nobody in the surrounding is immune*. But if a lot of people have gotten sick already, or when there is a vaccine, then a lot of people are immune, and less new people will acquire the disease. Effectively the reproduction will be less.

When the disease spreads faster, then this point of satiation/immunity will also be reached more quickly. Currently, strategies are being developed to make the disease spread as fast as possible (to get over with it quickly) but not too fast (to not overload the intensive care units and health care system). This may take a long time (there is still a lot of uncertainty about it) and possibly not economically viable.

If you see that the spreading might be over more quickly when we do nothing about it, then it could be better to just let the disease continue without measures, accept the high amount of deaths, and be done with it quickly instead of letting it linger for a long time. (But not many would dare to suggest these things seriously, because of the complicated ethical issues, and a stigmatization of an utilitarian approach to health and preventing deaths. )
> 4 votes

---
Tags: statistics, covid-19, sars-cov-2, epidemiology
---


# Question
Title: A .gtf file for HG38 for HISAT2 indexing
I am trying to build a HISAT2 index using Pertea et al. procedure (https://www.ncbi.nlm.nih.gov/pubmed/27560171) using their `extract_splice_sites.py` and `extract_exons.py` pieces of code. For that, I need a `.gtf` file of Human Genome HG38. I have obtained this file using recommendations described here http://www.metagenomics.wiki/tools/fastq/ncbi-ftp-genome-download

The file I have downloaded is `GCF_000001405.39_GRCh38.p13_genomic.gff` (1,1 Gb) and when I use on it this command `python3 extract_splice_sites.py GCF_000001405.39_GRCh38.p13_genomic.gff > test.ss` my `test.ss` file is 0 Kb! I wonder if there's another `.gff` file that I must use and where I can get it? This is my first time working with the human genome so I am not sure if I am doing the procedure the right way.

# Answer
You can download a GTF file from NCBI now. For hg38 (also known as GRCh38), just go to NCBI Assembly and search for 'GRCh38'. Click on the first result to go to this page and use the blue 'Download Assemblies' button to download data. Here, after you choose 'RefSeq' as your source of data, you will see a GTF file to download from the drop-down list, as shown in the screenshot below:
> 0 votes

# Answer
I had the same problem and I was able to solve it with that https://github.com/NBISweden/GAAS/blob/master/annotation/knowledge/gff\_to\_gtf.md
> 1 votes

---
Tags: python, gtf, human-genome
---


title: Filtering genes from cuffdiff results

text:
# Question
Title: Filtering genes from cuffdiff results
*This question has also been asked on Biostars*

I have run ***cuffdiff*** (with statistics turned ON) to compare two groups of samples: *Control group* and *Late AD group*.

This is the command I ran to be precise:

```
cuffdiff -L Control,AD_Late_Braak -p 8 --total-hits-norm --frag-bias-correct ../References/ensembl.GRCh38.99.fa --multi-read-correct --library-norm-method quartile ../References/Homo_sapiens.GRCh38.99.chr.gtf Early_Braak_Control_1,Early_Braak_Control_2,Early_Braak_Control_3 Late_Braak_Sample_1,Late_Braak_Sample_2,Late_Braak_Sample_3

```

I'm looking at the output from *cuffdiff* and I see a ***gene\_exp.diff*** file which contains the results of the differential expression testing. I want to know what is the best way to filter the results of this gene\_exp.diff file so as to restrict the number of genes that are up-regulated and down-regulated to a list between ~50-850.

P.S My thoughts are to tweak the **p-values** and **log2 fold change** values, but seems like a trail-and-error method, so I was wondering if there was a more "formal" method/approach?

# Answer
> 1 votes
To filter your gene\_exp.diff either you export it to excel and then filter because it is columns separated by tabs or use the command (awk and head -n) to filter and have the number of lines you want if you master Linux

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