title: Showing name of genes in boxplot

text:
# Question
Title: Showing name of genes in boxplot
I have a long format data frame like this

```
> head(df)
     gene

1:  CCL28   2.030674       2.023333 1.921508    1.952007
2:  CXCL2   2.204847       2.275296 2.175790    2.241554
3:  CXCR4   2.121186       2.147181 1.950065    1.982639
4: IFNGR1   2.311390       2.324111 2.220727    2.255556

```

And by this code

```
ggplot(df, aes(variable,value)) +
  stat_boxplot(geom="errorbar", width=.5)+
  geom_boxplot(aes(fill=variable))+
  theme_bw()+
  theme(axis.title.x=element_blank(), axis.title.y=element_blank())+
  stat_summary(fun.y=median, colour="red", geom="line", aes(group = 1))+
  geom_jitter(position = position_jitter(0.2))

```

I have this box

I have added this line to my code where I got

```
geom_text(data = df_melt, aes(x = variable, y = value, label = gene), size = 5)

```

\[!\[enter image description here\]\[2\]\]\[2\]

I want to show the name of genes in box plot but I don't how, cam you help?

Thanks

# Answer
> 1 votes
In addition, to have several errors in your code for plotting (for example, you are using `df` for plotting and then `df_melt` for `geom_text` but you didnot mention what is it). Based on your graph, it seems that you have issues with the melting part.

However, these are quite secondary issues that are probably be fixed with the code below. with your example, there is one trick because you are using `jitter` to plot your points (that is not addressed by the link that @ATpoint provided in comments). This will randomly position the point at the right y value but with different x position in order to not overlap points.

So, when you will add your text label, the labeling won't fit to the position of points. But you can make the jitter reproducible by adding a `seed` (see this github discussion: https://github.com/tidyverse/ggplot2/pull/1996). However, `geom_jitter` function does not have the argument `seed`, so you can use `geom_point` function and the argument `position_jitter` instead.

So, using the few genes that you provided as an example. I have to reshape it (a step you probably did but you did not provide in your question), here I'm using `pivot_longer` function from `tidyr` package:

```
library(dplyr)
library(tidyr)
DF <- df %>% pivot_longer(., -gene, names_to = "variable", values_to = "value")

# A tibble: 16 x 3
   gene   variable       value
   <chr>  <chr>          <dbl>
 1 CCL28  Res_immune      2.03
 2 CCL28  Non_res_immune  2.02
 3 CCL28  Res_OBP         1.92
 4 CCL28  Non_res_OBP     1.95
 5 CXCL2  Res_immune      2.20
 6 CXCL2  Non_res_immune  2.28
 7 CXCL2  Res_OBP         2.18
 8 CXCL2  Non_res_OBP     2.24
 9 CXCR4  Res_immune      2.12
10 CXCR4  Non_res_immune  2.15
11 CXCR4  Res_OBP         1.95
12 CXCR4  Non_res_OBP     1.98
13 IFNGR1 Res_immune      2.31
14 IFNGR1 Non_res_immune  2.32
15 IFNGR1 Res_OBP         2.22
16 IFNGR1 Non_res_OBP     2.26

```

And then, for the plotting part, you can use `geom_text_repel` function of the `ggrepel` package to add labeling without overlapping your points:

```
library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)

```

**Warning**

If you have a long list of genes, it will be pointless to show all of them on the boxplot. It will makes things really illisible. If you want to show some particular genes of interest, make a subset of the original dataframe. A lot of posts on SO and here have addressed this subsetting.

**Reproducible data example**

```
structure(list(gene = c("CCL28", "CXCL2", "CXCR4", "IFNGR1"),
    Res_immune = c(2.030674, 2.204847, 2.121186, 2.31139), Non_res_immune = c(2.023333,
    2.275296, 2.147181, 2.324111), Res_OBP = c(1.921508, 2.17579,
    1.950065, 2.220727), Non_res_OBP = c(1.952007, 2.241554,
    1.982639, 2.255556)), row.names = c(NA, -4L), class = c("data.table",
"data.frame"), .internal.selfref = <pointer: 0x563d91b65350>)

```

---
Tags: r, ggplot2
---


# Question
Title: Is the quality score of fastq used somewhere besides trimming/fastqc?
In the fastq format every 4k+4-th line contains the positionwise qualityscore (ascii encoded):

```
@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65
...

```

**which program uses this information besides trimming tools (cutadapt) or quality control (fastqc)?**

I heard that segemehl (mapping) abandons this quality string in favor of running time (maybe other mapping tools too?) But what about abyss, trinity, ...?

UPDATE:

So i atrifically generated a fastq file containing 2 identical reads that perfectly match a sequence in some genome (C.fna). The two reads have the most extreme phred scores (only '!!'=highest sanger phread and only 'II'=lowest sanger phred)

Now using segemehl the mapping quality is identical (60)

```
@HD VN:1.0
@SQ SN:C_fna    LN:15440
@RG ID:A1   SM:sample1  LB:library1 PU:unit1    PL:illumina
@PG ID:segemehl VN:0.3.4    CL:segemehl -i C2.fna.idx -d C2.fna -q test.fq
A   0   C_fna   261 60  220=    *   0   0   ACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAGGTAACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCGATATTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCACCTGGTG    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    HI:i:0  NH:i:1  NM:i:0  MD:Z:220    RG:Z:A1 YZ:Z:0
B   0   C_fna   261 60  220=    *   0   0   ACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAGGTAACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCGATATTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCACCTGGTG    IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII    HI:i:0  NH:i:1  NM:i:0  MD:Z:220    RG:Z:A1 YZ:Z:0
[SEGEMEHL] Thu Feb  6 08:35:15 2020: matching w/ suffixarray has taken 0.000000 seconds.

```

So segemehl does NOT use the read quality phred (basically only using the read sequence for mapping)

# Answer
I think almost all alignment programs would use base quality as a factor to calculate mapping quality (See here). These would also be used in Variant Calling to pick reads where high quality bases confidently align as either the reference base or one of the alternate bases.
> 4 votes

# Answer
LoFreq is a variant calling tool that specifically leverages per-base quality scores.

For long read technologies (PacBio SMRT, ONT), I think quality score usage is mostly limited to qc steps as well, such as filtering out poor reads either by mean quality score of the whole read (as in ONT's basecalling softwares and Filtlong) or by the lowest quality score within a window in a read (Filtlong). I think most long read mappers ignore these values (for example minimap2 and ngmlr do and can take fasta for read inputs).
> 1 votes

---
Tags: rna-seq, fastq
---


title: how to get FlyBase ID conversion table

text:
# Question
Title: how to get FlyBase ID conversion table
FlyBase uses systematic naming of genes, transcripts and proteins in a format `FBtypeNUMBER`. However, the numbers are not the same for genes and their products. For example gene `FBgn0259101` transcribes to `FBtr0299513` and translates to `FBpp0288787`.

On this webpage users can conver individual IDs, but what If I need to convert practically whole transcriptome? Where can I get the full translation table?

# Answer
Indeed the web interface does not seem to link the full translation table, but it is available on their ftp server at `releases -> release_ID -> precomputed_files -> genes -> fbgn_fbtr_fbpp_fb_2019_06.tsv.gz`

For example, this is the latest (`FB2019_06`) translation table available.
> 0 votes

---
Tags: identifiers, conversion
---


# Question
Title: Can I submit a R package to Bioconductor or CRAN if I have already published it a journal?
I have written a bioinformatics package in R that I want to publish in a bioinformatics Journal. Presently, I am maintaining a local repo of that package and I want to put in the Bioconductor repository (after publication or at least after submission to the journal), is it possible? or I must submit it to Bioconductor first? Thanks

# Answer
> 5 votes
You can submit a package to a repo after the publication, but I would say that you really should do it before.

Bioconductor accept only packages that are not published on CRAN, however, academic publications are fine.

Conversely, Bioinformatics does not enforce any particular platform for sharing the code, but you have to make it available somehow (e.g. GitHub repo, university webpage, bioconductor). From the instructions for authors

> If the manuscript describes new software tools or the implementation of novel algorithms the software must be freely available to non-commercial users at the time of submission, and appropriate test data should be made available.

However, you as a developer should make all the effort to make your software available and easily accessible. Also, sooner you make it available, sooner people start to cite you, which is a great advantage for you. Finally, you really need users to make sure that your package really works (i.e. debugging all the corner cases).

I could not emphasise more, the correct order is 1. make it available 2. write a publication

---
Tags: bioconductor
---


# Question
Title: Why does the FASTA sequence for coronavirus look like DNA, not RNA?
I'm looking at a genome sequence for 2019-nCoV on NCBI. The FASTA sequence looks like this:

```
>MN988713.1 Wuhan seafood market pneumonia virus isolate 2019-nCoV/USA-IL1/2020, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC
TAATTACTGTCGTTGACAGGACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTG
...
...
TTAATCAGTGTGTAACATTAGGGAGGACTTGAAAGAGCCACCACATTTTCACCGAGGCCACGCGGAGTAC
GATCGAGTGTACAGTGAACAATGCTAGGGAGAGCTGCCTATATGGAAGAGCCCTAATGTGTAAAATTAAT
TTTAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGGAGAATGACAAAAAAAAAAAA

```

Coronavirus is an RNA virus, so I was expecting the sequence to consist of `AUGC` characters. But the letters here are `ATGC`, which looks like DNA!

I found a possible answer, that this is the sequence of a "complementary DNA". I read that

> The term *cDNA* is also used, typically in a bioinformatics context, to refer to an mRNA transcript's sequence, expressed as DNA bases (GCAT) rather than RNA bases (GCAU).

However, I don't believe this theory that I'm looking at a cDNA. If this were true, the end of the true mRNA sequence would be `...UCUUACUGUUUUUUUUUUUU`, or a "poly(U)" tail. But I believe the coronavirus has a poly(A) tail.

I also found that the start of all highlighted genes begin with the sequence `ATG`. This is the DNA equivalent of the RNA start codon `AUG`.

So, I believe what I'm looking at is the true mRNA, in 5'→3' direction, but with all `U` converted to `T`.

**So, is this really what I'm looking at? Is this some formatting/representation issue? Or does 2019-nCoV really contain DNA, rather than RNA?**

# Answer
> 49 votes
That is the correct sequence for 2019-nCov. Coronavirus is of course an RNA virus and in fact, to my knowledge, every RNA virus in Genbank is present as cDNA (AGCT, i.e. thydmine) and not RNA (AGCU, i.e. uracil).

The reason is simple, we never sequence directly from RNA because RNA is too unstable and easily degraded by RNase. Instead the genome is reverse transcribed, either by targeted reverse transcription or random amplification and thus converted to cDNA. cDNA is stable and is essentially reverse transcribed RNA.

The cDNA is either sequenced directly or further amplified by PCR and then sequenced. Hence the sequence we observe is the cDNA rather than RNA, thus we observe thymine rather than uracil and that is how it is reported.

# Answer
> 17 votes
Most sequencing experiments, be it Illumina-based next-generation-sequencing or Sanger sequencing uses DNA as template, not RNA. Even if this virus is RNA-based it would be reverse-transcribed prior to any sequencing experiment. Therefore the output is DNA and this is what NCBI provides here.

# Answer
> 5 votes
> If this were \[cDNA\], the end of the true mRNA sequence would be ...UCUUACUGUUUUUUUUUUUU, or a "poly(U)" tail.

A cDNA sequence, maybe confusingly, refers to the *coding strand* of the cDNA (despite being called “complementary”). So while cDNA is the result of reverse transcribing RNA into DNA, by convention it has the same strandedness as the original RNA. That’s why what you’re seeing is read in 5′→3′ direction and contains a visible poly(A) tail. Having a single conventional reading direction for all archived sequences vastly simplifies data handling, and reduces errors.

In fact, since cDNA is double-stranded, there is no a priori reason why a computer-stored cDNA sequence should refer to the template strand (i.e. the opposite strand, which is synthesised from the RNA during reverse transcription).

The whole (simplified) synthesis process of cDNA is as follows:

1. A primer hybridises to the template RNA molecule.
2. The RNA template is reverse transcribed into DNA using reverse transcriptase.
3. The RNA template is removed.
4. A complementary strand is transcribed along the (currently) single-stranded cDNA, resulting in a double-stranded cDNA product.

# Answer
> 3 votes
It's not common to sequence directly from RNA because most sequencing platforms don't have that as an option. Nanopore sequencers do allow this, but I'm not aware yet of any 2019-nCov preprints involving nanopore RNA sequencing. I expect that will change in the next month or so.

Commercial kits exist; there are no insurmountable technical issues with it. Direct RNA sequencing can be done locally, on-site near the point of discovery without sample transfer or culture on a USB-powered device that fits in a pocket (RNA preparation takes about 2 hours). Flow cells that have potentially-infectious RNA inside them can be disposed as biohazard waste. However, the ease at which RNA can be quickly converted to more stable cDNA then amplified to create a much higher concentration DNA sample (which is quicker / more efficient to get results from) means that cDNA is generally preferred for sequencing unless the native RNA is needed (e.g. for looking at RNA base modifications that are destroyed when converting to cDNA).

There is a paper on coronavirus direct RNA sequencing with nanopore here; I would expect that 2019-nCoV would have a similar difficulty. The zika virus has an extremely low viral load in human blood, but has also been sequenced via direct RNA sequencing of \[carefully\] cultured cells (see here).

Regardless of whether or not RNA sequencing has actually been carried out, most genetic data analysis programs will only work with A/C/G/T sequences, so it's conventional to replace any U parts of an RNA sequence with T for data storage. There's no loss of information by doing this, as T replaces all Us in the RNA sequence.

---
Tags: fasta, phylogenetics, ngs, phylogeny, covid-19
---


title: Griding ggplot object in r

text:
# Question
Title: Griding ggplot object in r
I have this data frame

> head(d)

```
# A tibble: 6 x 3
  gene  variable       value
1 CCT5  Res_immune     10.6
2 CCT5  Non_res_immune 10.5
3 CDK1  Res_immune      9.25
4 CXCL5 Non_res_immune  9.82
5 GBP1  Res_immune      9.04
6 GBP1  Non_res_immune  9.43
>

```

And this code

```
library(ggplot2)
library(ggrepel)
ggplot(d, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)

```

The genes in first column are involved in cytokine signaling pathway in two groups of patients

This is my plot

I have another data frame for genes in IL17 signaling in the sample patients like

head(d)

```
# A tibble: 6 x 3
  gene  variable       value
  <chr> <chr>          <dbl>
1 CDK6  Res_immune     17.6
2 IL6  Non_res_immune 14.5
3 CXCL8  Res_immune      1.25

```

How I can put plots from both data in a grade like below where the name of pathways shown instead of DEL, INS, SNP and total like below picture?

# Answer
What you are looking for is `facetting`. Using this key word on any search engines you will find dozen of answers describing its use: https://ggplot2.tidyverse.org/reference/facet\_grid.html and https://ggplot2.tidyverse.org/reference/facet\_wrap.html

To prepare the faceting in `ggplot2`, a possible way is to bind both dataframes together by first specifying which group they belong such as:

```
cytokine$Group = "Cytokine"
il17$Group = "IL17"
DF <- rbind(cytokine,il17)
DF

    gene       variable value    Group
1:  CCT5     Res_immune 10.60 Cytokine
2:  CCT5 Non_res_immune 10.50 Cytokine
3:  CDK1     Res_immune  9.25 Cytokine
4: CXCL5 Non_res_immune  9.82 Cytokine
5:  GBP1     Res_immune  9.04 Cytokine
6:  GBP1 Non_res_immune  9.43 Cytokine
7:  CDK6     Res_immune 17.60     IL17
8:   IL6 Non_res_immune 14.50     IL17
9: CXCL8     Res_immune  1.25     IL17

```

Then, use `facet_grid` function of `ggplot2`:

```
library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)+
  facet_grid(Group~., scales = "free")

```

Alternatively, you can use `facet_grid` (or `facet_wrap`) to have your plot in a single row:

```
library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)+
  facet_wrap(.~Group, scales = "free")

```

**Reproducible data example**

```
structure(list(gene = c("CCT5", "CCT5", "CDK1", "CXCL5", "GBP1",
"GBP1"), variable = c("Res_immune", "Non_res_immune", "Res_immune",
"Non_res_immune", "Res_immune", "Non_res_immune"), value = c(10.6,
10.5, 9.25, 9.82, 9.04, 9.43), Group = c("Cytokine", "Cytokine",
"Cytokine", "Cytokine", "Cytokine", "Cytokine")), row.names = c(NA,
-6L), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x55bfc0db9350>)

```

and

```
structure(list(gene = c("CDK6", "IL6", "CXCL8"), variable = c("Res_immune",
"Non_res_immune", "Res_immune"), value = c(17.6, 14.5, 1.25),
    Group = c("IL17", "IL17", "IL17")), row.names = c(NA, -3L
), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x55bfc0db9350>)

```
> 4 votes

---
Tags: r, ggplot2
---


# Question
Title: No output when trying to obtain protein sequence from PDB file
I am trying to obtain the sequence from PDB file, but no output is given when I run the following code :

```
from Bio import SeqIO

with open("1a9l.pdb", "r") as handle:
    for record in SeqIO.parse(handle, "pdb-atom"):
        print(record.seq)

```

# Answer
If you are using a conda version of python you can use Pymol module for this, which is actually more robust in some aspects than `bio.PDB`, although errors may cause segmentation faults. The following:

```
import pymol2
## paralellisable version.
with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1a9l')
    print(pymol.cmd.get_fastastr('polymer.nucleic'))

```

Gives me:

```
>1a9l_A
GGGUGACUCCAGAGGUCGAGAGACCGGAGAUAUCACCC

```

Alternatively, you can use PDBe's API.

```
import requests

response = requests.get('https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/1a9l').json()['1a9l']

for entity in response:
    if entity['molecule_type'] == 'polyribonucleotide':
        print(entity['sequence'])

```

Do note that missing 'residues' (residues, small molecules and nucleotides) will cause issues. PDB files do not have missing 'residues', while mmCIF files have a dictionary that does. The API will give you the latter.
> 2 votes

# Answer
That PDB entry is an RNA molecule, so there is no protein sequence to extract.

If it extracted the RNA sequence by default here you can see there could be an ambiguity as to whether it was an RNA sequence or a protein sequence, as the alphabets share characters.
> 3 votes

---
Tags: python
---


title: Running SnakeMake on cluster

text:
# Question
Title: Running SnakeMake on cluster
I do not understand how to specify correctly the parameters on a `SLURM` cluster for `snakemake` to use them. I tried submitting the following `SLURM` file, but it does not work this way and the number of cores used is only `1`, not `20`:

```
#!/bin/bash

#SBATCH -p standard
#SBATCH -A overall
#SBATCH --time=12:00:00
#SBATCH --output=snakemake%A.out
#SBATCH --error=snakemake%A.err
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=20
#SBATCH --mem=120000

snakemake

```

Then, I tried following `snakemake` tutorial. And created `cluster.json` based on the `SLURM` parameters that I need:

```
{
    "__default__" :
    {
        "A" : "overall",
        "time" : "24:00:00",
        "nodes": 1,
        "ntasks": 1,
        "cpus" : 20,
        "p" : "standard",
        "mem": 120000,
        "output": "snakemake%A.out",
        "error": "snakemake%A.err"
    }
}

```

And ran `snakemake` inside a newly created `snakemake.sh` script:

```
#!/bin/bash

snakemake -j 999 --cluster-config cluster.json --cluster "sbatch -A
{cluster.A} -p {cluster.p} \
-t {cluster.time} --output {cluster.output} --error {cluster.error} --
nodes {cluster.nodes} \
--ntasks {cluster.ntasks} --cpus-per-task {cluster.cpus} --mem
{cluster.mem}"

```

And it is giving me an error now:

> sbatch: error: Unable to open file
> /bin/sh: line 1: -t: command not found Error submitting jobscript (exit code 127):

I am now completely lost at what I shoould actually do. I would prefer plain regular `.slurm` file submission, but how to make `snakemake` use them? Any suggestions would be greatly appreciated.

I removed `\` \- line separators in the `snakemake.sh` script:

```
#!/bin/bash

snakemake -j 10 --cluster-config cluster.json --cluster "sbatch -A {cluster.A} -p {cluster.p} -t {cluster.time} --output {cluster.output} --error {cluster.error} --nodes {cluster.nodes} --ntasks {cluster.ntasks} --cpus-per-task {cluster.cpus} --mem {cluster.mem}"

```

And it started to run. It is not convenient though for me. I would rather prefer submitting just one job using `.slurm` file passing all of the parameters from `#SBATCH`. Is it possible?

# Answer
You can omit `--nodes`, you need the following:

```
#!/bin/bash
#SBATCH --ntasks-per-node=1
#SBATCH -c threads_from_snakemake
#SBATCH -p some_partition

Commands go here

```

For slurm, you may want to modify my SlurmEasy script. We use that with snakemake all the time, usually of the form:

```
snakemake -c "SlurmEasy -t {threads} -n {rule} --mem-per-cpu {cluster.memory}" ...

```

You'll also want units with your memory requests, like `12G`. Also, in general it's best not to submit snakemake itself to the cluster, rather run it on an interactive node and have it submit jobs on your behalf (if the admins complain point out that snakemake is barely using any resources on the head node).
> 10 votes

# Answer
While I agree with not submitting snakemake to the cluster, it's not a great solution to make your pilot program run interactively - automation via batch script should not be dependent on an open shell.

How about allocating a whole node with many cores and working within the GNU parallel abilities of snakemake.
> 0 votes

---
Tags: snakemake, slurm
---


# Question
Title: VCF's unknown ('.') values causing problems for GATK's Liftover Vcf
I am having a tough time calling LiftoverVcf in GATK. I have concatenated all of the Vcfs of interest and removed rows with problematic values but now when I run LiftoverVCF I get the following error:

```
Intel; Inflater: Intel; Provider GCS is available; Picard version: Version:4.0.12.0
INFO 2020-02-10 14:36:03 LiftoverVcf Loading up the target reference genome.
INFO 2020-02-10 14:36:16 LiftoverVcf Lifting variants over and sorting (not yet writing the output file.)

[Mon Feb 10 14:59:31 EST 2020] picard.vcf.LiftoverVcf done. Elapsed time: 23.50 minutes.
Runtime.totalMemory()=7414480896
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
java.lang.IllegalStateException: Key . found in VariantContext field INFO at 1:121387974 but this key isn't defined in the VCFHeader. We require all VCFs to have complete VCF headers by default.
at htsjdk.variant.vcf.VCFEncoder.fieldIsMissingFromHeaderError(VCFEncoder.java:202)
at htsjdk.variant.vcf.VCFEncoder.write(VCFEncoder.java:141)
at htsjdk.variant.variantcontext.writer.VCFWriter.add(VCFWriter.java:248)
at picard.vcf.LiftoverVcf.rejectVariant(LiftoverVcf.java:456)
at picard.vcf.LiftoverVcf.doWork(LiftoverVcf.java:357)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:295)
at org.broadinstitute.hellbender.cmdline.PicardCommandLineProgramExecutor.instanceMain(PicardCommandLineProgramExecutor.java:25)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
at org.broadinstitute.hellbender.Main.main(Main.java:289)

```

Not sure why I am getting this error other than there may be a `.` in one of the info fields. Here is what the first several rows of my VCF looks like minus the header fields:

```
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT AltaiNea
1 10001 . T . 51.05 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=6.32;MQ0=242 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:21.05:0,21,265:1,2:0,1:1,0:116,129:2
1 10002 . A . 36.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.31;MQ0=248 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:6.01:0,6,69:171,184:0,9:0,0:0,2:0
1 10003 . A . 36.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.31;MQ0=248 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:6.02:0,6,72:222,201:1,0:0,0:0,8:0
1 10004 . C . 39.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.61;MQ0=247 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:9.02:0,9,97:0,1:268,211:0,0:0,1:0
1 10005 . C . 41.99 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=4.31;MQ0=246 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:12:0,12,119:0,0:341,287:0,0:0,0:0
1 10006 . C . 36.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.31;MQ0=248 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:6.02:0,6,72:0,1:408,347:0,0:2,0:0
1 10007 . T . 39.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=4.21;MQ0=246 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:9.02:0,9,97:0,0:0,0:0,0:506,451:0
1 10008 . A . 33.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=4.05;MQ0=247 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:3.01:0,3,40:547,507:0,0:0,0:0,0:0
1 10009 . A . 36.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.31;MQ0=248 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:6.02:0,6,70:600,534:0,1:0,0:0,1:0
1 10010 . C . 42.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=4.72;MQ0=244 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:12.02:0,12,125:0,0:649,547:0,1:1,0:0
1 10011 . C . 51.02 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=5.12;MQ0=242 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:21.03:0,21,209:0,2:732,612:0,1:0,0:0
1 10012 . C . 44.99 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.94;MQ0=243 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:14.99:0,15,141:0,3:804,657:0,0:6,1:0
1 10013 . T . 57.04 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.20;MQ0=238 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:27.04:0,27,285:0,0:1,0:0,0:904,760:0
1 10015 . A . 60.04 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.19;MQ0=238 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:30.04:0,30,301:985,850:0,0:1,0:0,0:0
1 10016 . C . 72.03 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.99;MQ0=235 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:42.04:0,42,427:2,1:1024,852:0,0:1,0:0
1 10017 . C . 84.05 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=8.37;MQ0=232 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:54.06:0,54,534:0,2:1100,912:0,0:0,0:0
1 10018 . C . 108.02 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=8.87;MQ0=224 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:78.03:0,78,757:0,0:1167,978:0,1:2,0:0
1 10019 . T . 96.05 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=10.24;MQ0=219 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:66.05:0,66,646:1,0:1,0:1,0:1258,1090:0
1 10020 . A . 65.97 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=9.74;MQ0=221 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:35.97:0,36,331:1309,1150:0,0:1,1:0,0:0
1 10021 . A . 65.34 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=8.14;MQ0=228 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:35.35:0,35,489:1339,1179:3,0:0,0:1,0:0
1 10022 . C . 92.98 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.79;MQ0=229 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:62.99:0,63,590:0,1:1389,1188:0,1:0,0:0
1 10023 . C . 98.93 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.69;MQ0=227 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:68.94:0,69,636:2,2:1453,1242:0,0:1,2:0
1 10024 . C . 98.95 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.46;MQ0=227 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:68.96:0,69,620:1,3:1497,1302:1,0:3,2:0

```

Is there a way for me to add `.` into the VCF header as a tag for all of the appropriate fields so that I can bypass this problem? I know from reading online that the `.` means that value is unknown. I would like to figure out how to instruct GATK to simply accept it regardless of the fact that it cannot be found in the VCF header.

Also, I am open to other ways I can liftover this genome to hg38. Below is the command I am using to do so using GATK.

`java -jar ~/work/progs/gatk-4.0.12.0/gatk-package-4.0.12.0-local.jar LiftoverVcf -I AltaiNea.hg19_1000g.concat.LowQualRemoved.vcf -O AltaiNea.hg38_1000g.concat.LowQualRemoved.vcf -C hg19ToHg38.over.chain -R hg38.fa --REJECT rejectedRecords.vcf`

# Answer
> 1 votes
My problem was solved using --ALLOW-MISSING-FIELDS-IN-HEADER.

---
Tags: vcf, vcftools
---


# Question
Title: How to calculate lipophilicity of proteins?
The question says it all. I want to calculate lipophilicity of proteins to screen the compounds out for antigen-antibody interactions.

# Answer
## Oligopeptide

The lipophilicity of an oligopeptide or a small molecule can be calculated with RDKit, the classic compChem python toolkit. The measure is logP or partition coefficient, namely the magnitude of the ratio of the concentration in octanol over than in water, so the higher it is, the more it partitions in the octanol phase. However, with a compound of several dalton it may be good to see what the contribution of each atom to this score are. This measure is the Crippen contribution to logP, which can be handily plotted on the molecule. Note that these are may be at first glance look similar to Gastaiger partial changes, but are not.

Okay, first the imports

```
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from rdkit.Chem.Draw import SimilarityMaps

```

Now let's make the molecules. Let's start with something we can simply use a SMILES string of the web: Bortezomib.

```
name = 'Bortezomib'
mol = Chem.MolFromSmiles('O=C(N[C@H](C(=O)N[C@H](B(O)O)CC(C)C)Cc1ccccc1)c2nccnc2')

```

Now we can do all the relevant calculations on the `mol` object:

```
print('logP', Descriptors.MolLogP(mol))
AllChem.Compute2DCoords(mol)
contribs = Chem.rdMolDescriptors._CalcCrippenContribs(mol)
fig = SimilarityMaps.GetSimilarityMapFromWeights(mol, [x for x,y in contribs], colorMap='BuPu', contourLines=10)
fig.savefig(name+'_crippen.png', bbox_inches='tight')

```

LogP=0.3

`colorMap='BuPu'` is the matplotlib colorscheme. I have no idea which work best, I have spent way too much time trying and I do not know. But with BuPu, white is hydrophilic, purple is hydrophobic.

Say the oligopeptide does not have a SMILES or we don't want to do it that way. Let's take β-Casomorphins 1–3:

```
name = 'β-Casomorphins 1–3'
mol = Chem.MolFromFASTA('YPF')
mol = AllChem.ReplaceSubstructs(mol,
                                Chem.MolFromSmiles('N'),
                                Chem.MolFromSmiles('[NH][OH]'),
                                replacementConnectionPoint=0)[0]

```

The result with the last block from before gives:

logP=1

## Proteins

With protein, it is a lot harder. Membrane affinity of protein is dependent on their structure. In PyMOL you can use the Adaptive Poisson-Boltzmann Solver (APBS) to show the surface by charge. For transmembrane protein there is clear ring of hydrophobicity —CHARMM-GUI membrane builder and a variety of tools for example can guess where it is. For a membrane associated protein only one side will be hydrophobic with a tell-tell charged ring around it (e.g. DHR1 domain) and there are no hard and tested methods for that.
> 2 votes

---
Tags: protein-protein-interaction
---


# Question
Title: Question: How to simulate 100k samples having 40 million SNPs in a proportion of case:control=30:70?
*Note: this question can also be found on Biostars*

I need to perform a stress test in a GWAS tool and the duty demands a dataset (plink format) having 100 thousand samples, having 40 million SNPs in a proportion of case:control=30:70. I'm performing the command:

> plink1.9 --simulate ds1.sim --make-bed --out ds1 --simulate-ncases 30000 --simulate-ncontrols 70000

But after about 3 hours execution, I'm facing `Error: File write failure.`, even though my environment has 256T of disk. (It's a container with a S3 (Object store) mounted through s3fs).

I was wondering if there's a clever way to do that? Perhaps breaking it into small chunks and merging up later somehow? It would prevent me to be waiting hours till an error happens and I loose either time and the sample generated so far.

Thanks in advance!

# Answer
> 1 votes
First comment: you should be using real data, rather than simulated data, as we are not very good at modelling genetic variation and recombination at a genome-wide scale. A way that this can be done with GCTA is mentioned in the PLINK manual.

I'll do some calculations to estimate the size of data you need to store:

Assume that this plink-formatted file is in tped format (one row per SNP), and the entire genome is diploid. That means that a single SNP/individual pair will take up three bytes. You want 30,000 cases and 70,000 controls, so each line from each sample will take up 300kB (plus a bit of change for the header). With 1 thousand SNPs, that's 300MB per sample; with 1 million SNPs, that's 300GB per sample; with 10 million SNPs, that's 3TB per sample; with 40 million SNPs, that's 12 TB for one sample.

Reading the plink manual, I see that the BED format is a two-bit compression for each SNP/individual pair, so one byte per 4 individuals, or about 12x compression (i.e. 1TB for one sample).

Okay, that's a lot. It's well within a 256TB storage limit, but might hit file size limits depending on the file system.

If there is a data generation issue with PLINK, it might be worth trying (as you suggest) splitting the input for simulation then merging afterwards. My guess is that doing the simulation on fewer SNPs (rather than fewer individuals) would work better. Plink has a `--bmerge` option for merging the compressed binary file formats.

I usually approach problems like these by working my way up from something that I know works. Try first merging two datasets of 1000 SNPs. If that works, try larger datasets. Rinse and repeat until the computer crashes, you get bored, or you reach your goal.

---
Tags: gwas, plink
---


title: Vegetable databases usage

text:
# Question
Title: Vegetable databases usage
Good morning, I'm aligning with blast some species, against all nt database to analyse the vegetables. However, in some cases we have species of animals mixed with vegetables and want only to detect the vegetable species. Normally, we would only expect to detect the vegetables only. But this is normal since we are using the all nt. So now, could anyone point me out, the way to detect only vegetables? Thanks for the help.

# Answer
> 1 votes
Using the blast command line tool, you can limit your search to a specific taxon using the `–taxids` option.

So, in your case `-taxids 3193` should limit your search to plants.

---

If you are using the ncbi blast webinterface, you can restrict your search by entering "plants" into the Organism field in the "Choose Search Set" section

---
Tags: database, nt
---


# Question
Title: How to visualize called narrowPeak files in UCSC Genome browser or IGV?
I have called peaks using MACS2. Then I got a narrowPeak file like this.

```
1   4491713 4491913 testONE_peak_1  37  .   2.86003 3.72580 1.66717 80
1   4492782 4493153 testONE_peak_2  27  .   2.22933 2.71107 0.82442 242
1   4571554 4572026 testONE_peak_3  72  .   4.16005 7.21000 4.85200 282
1   4739486 4739879 testONE_peak_4  64  .   3.63128 6.42948 4.14394 275
1   4784848 4785061 testONE_peak_5  31  .   2.24218 3.19160 1.21569 133

```

Then I tried to visualize it in UCSC Genome browser and I got following error.

```
Error File 'testONE_peaks.narrowPeak' - Unrecognized format line 1 of file: 1 4491713 4491913 testONE_peak_1  37  .   2.86003 3.72580 1.66717 80 (note: chrom names are case sensitive, e.g.: correct: 'chr1', incorrect: 'Chr1', incorrect: '1')

```

So I added chr to the chromosome column instead of 1 (Chromosome number).

```
chr1    4491713 4491913 testONE_peak_1  37  .   2.86003 3.7258  1.66717 80
chr1    4492782 4493153 testONE_peak_2  27  .   2.22933 2.71107 0.82442 242
chr1    4571554 4572026 testONE_peak_3  72  .   4.16005 7.21    4.852   282
chr1    4739486 4739879 testONE_peak_4  64  .   3.63128 6.42948 4.14394 275
chr1    4784848 4785061 testONE_peak_5  31  .   2.24218 3.1916  1.21569 133

```

When I tried to load the file again I got following error.

```
Error File 'testONE_peaks_NEW2.narrowPeak' - Error line 1 of custom track: thickStart out of range (chromStart to chromEnd, or 0 if no CDS)

```

Then I added following line to the top and I could load the file to UCSC Genome browser.

```
track type=narrowPeak name="Somite narrowPeak"
chr1    4491713 4491913 testONE_peak_1  37  .   2.86003 3.7258  1.66717 80
chr1    4492782 4493153 testONE_peak_2  27  .   2.22933 2.71107 0.82442 242
chr1    4571554 4572026 testONE_peak_3  72  .   4.16005 7.21    4.852   282
chr1    4739486 4739879 testONE_peak_4  64  .   3.63128 6.42948 4.14394 275
chr1    4784848 4785061 testONE_peak_5  31  .   2.24218 3.1916  1.21569 133

```

But the problem is I can not see anything in the Genome browser. Only the name is there. But nothing else. Do you know what I should do?

# Answer
> 2 votes
I think you might have changed the separator (or at least have some kind of inconsistency from the required format) for your file. Note that peak output files from MACS2 are variants of BED files.

It seems you **need to have a *tab* separator for this type of file.**

I copied your example file and then ran `awk -v OFS="\t" '$1=$1' your_peaks.narrowPeaks > fixed_your_peaks.narrowPeaks` on it. This produced a result that I was able to see in both the UCSC genome browser and IGV.

---
Tags: rna-seq, chip-seq, ucsc, igv, macs2
---


title: Running tophat dockerfile in background

text:
# Question
Title: Running tophat dockerfile in background
I need to use tophat2 to align some sequences, and I wish to use the docker container. These are large sequences so it will take a long time, plus I'm working on a university server so the chances of being kicked off the server are very high. How can I use nohup and & to run tophat2 through docker, if I need some user interaction?

# Answer
> 2 votes
Assuming your user is in the group permitted to use docker:

```
docker pull quay.io/biocontainers/tophat:2.1.1--py27_3
nohup docker run -v /some/path:/data quay.io/biocontainers/tophat:2.1.1--py27_3 tophat2 ... &

```

The only thing to remember is to mount everything you'll need into the docker container (`/data` in the example above). If you need to debug things then you can add the `-it` flag to `docker run` and use `bash` as the command.

As an aside, you might find using `screen` or `tmux` preferable to `nohup ... &`, since you can reconnect to sessions with them.

---
Tags: sequence-alignment, bash, docker
---


# Question
Title: How to exclude the repetition of gene-gene correlation calculation in python?
I am trying to make a paired matrix of gene-gene correlation. Considering that I have a huge matrix (13000 genes and 900 samples) and for some reasons I don't want to decrease the number of my genes, my gene-correlation matrix would be 13000\*13000 and my paired matrix will become 169 million \*4 (Column 1: Gene 1; Column 2 : Gene 2; Column 3: Correlations; Column 4: P-values) . In this case, I have to exclude unnecessary calculations as much as I can. I have excluded the situation that Gene 1 = Gene 2. But I couldn't find a way to exclude the condition that "Column 1: Gene 1 ; Column 2: Gene 2 = Column 1: Gene 2; Column 2: Gene 1 ". To make a long story short, correlation between G1 and G2 is equal to G2 and G1. It is like calculating just lower section of diag in a symmetric matrix. I would be grateful if anybody help me in this case. I have enclosed my python codes here:

```
import pandas as pd
import numpy as np
import scipy
import math
import openpyxl
from openpyxl import Workbook
from scipy.stats import spearmanr

din=pd.read_csv('m_test.csv', index_col=0)

out=pd.DataFrame()
outdf=pd.DataFrame()

for g1 in din.index:
    for g2 in din.index:
        temp=din.loc[[g1, g2]]
        if g1==g2:
            next
        else:
            spR, spP=spearmanr(temp.loc[g1], temp.loc[g2])
            frame={'g1':[g1], 'g2':[g2], 'spR':[spR], 'spP':[spP]}

            out=pd.DataFrame(frame)

            outdf=pd.concat([outdf, out])

```

# Answer
> 1 votes
You don't need `g2` to go from `0` to `din.index`, just from `0` to `g1 - 1`. This way, you'll end up calculating just for the "lower triangle":

```
g1-> 0   1   2   3
g2
0                    #(0,0 to 0,1-1 = nothing)
1    X               #(1,0 to 1,1-1 = 1,0)
2    X   X           #(2,0 to 2,2-1 = 2,0; 2,1)
3    X   X   X       #(3,0 to 3,3-1 = 3,0; 3,1 and 3,2)
4    X   X   X   X   #(4,0 to 4,4-1 = 4,0; 4,1; 4,2 and 4,3)

```

# Answer
> 1 votes
Actually I found an easy and efficient way.

```
import pandas as pd
import numpy as np
import scipy
import math
import openpyxl
from openpyxl import Workbook
from scipy.stats import pearsonr
from scipy.stats import spearmanr
from scipy.stats import kendalltau

pro=pd.read_csv('m_test.csv', index_col=0)

xInds = []
yInds = []
for i in range(len(pro.index)):
    for j in range(i+1, len(pro.index)):
        xInds.append(i)
        yInds.append(j)

z=0
Rlist_sp = []
Plist_sp = []
Rlist_pe = []
Plist_pe = []
while z < len(xInds):
    b=xInds[z]
    c=yInds[z]

    spR, spP = spearmanr(pro.iloc[b].values, pro.iloc[c].values)
    peR, peP = pearsonr(pro.iloc[b].values, pro.iloc[c].values)

    Rlist_sp.append(spR)
    Plist_sp.append(spP)
    Rlist_pe.append(peR)
    Plist_pe.append(peP)
    z=z+1

G1=pd.DataFrame(xInds)
G2=pd.DataFrame(yInds)
R_sp=pd.DataFrame(Rlist_sp)
P_sp=pd.DataFrame(Plist_sp)
R_pe=pd.DataFrame(Rlist_pe)
P_pe=pd.DataFrame(Plist_pe)

Final=pd.concat([G1, G2, R_sp, P_sp, R_pe, P_pe], axis=1)
np.savetxt('Final_gg_A549.txt', Final.values, fmt='%s', delimiter='\t')

```

---
Tags: python, statistics, gene-expression, correlation
---


# Question
Title: How to find a homolog that has a PDB structure available
I tried to find a homolog that has a PDB structure available so I can use this PDB file for comparative modeling. I ran a BLAST search but none of the search results seemed to have a known structure. Is there any way to find the results I want?

# Answer
If you are using NCBI Blast and blastp against the pdb database gives you nothing, you can do a psiblast against the normal refseq database, do a few rounds, save the PSSM matrix and submit a new job with the PSSM matrix but against the PDB database.

You might need to limit the species of you get lots of too close homologues —say blasting E. coli protein for which you'd want to exclude Enterobacteriales.
> 1 votes

# Answer
You could use HMM-based comparison, for example with HHpred. This is generally the best way to find structural templates because structure is more conserved than sequence, and comparing HMMs will turn up hits of the same structure but minimal sequence overlap.

If HHpred turns up nothing you are in the region of de novo modelling, which is a different problem. However, it has advanced in recent years to the point where it may be useful to you. For example, you could try our method DMPfold.
> 2 votes

# Answer
To expand the first answer by @MatteoFurla: re-run the P-blast and under database option, which is set by default to non-redundant database (nr) select Protein databank (PDB), the third option from the bottom of the menu. All the Blast results will be from PDB and thus all results will have associated structures. Also note that the Genbank codes will be PDB codes.

If you have zero hits then consider PSI-Blast, how I suspect you are running P-Blast via the nr database.
> 1 votes

---
Tags: blast, pdb
---


# Question
Title: Make a fasta out of mapped reads without taking into account the reference
I've aligned reads onto an haploid reference genome. I'd like to have a consensus sequence of all my aligned reads that doesn't take into account the reference genome I used (i.e. I just want the nucleotides from my reads at each base (taking into account frequency and quality) but if there is no read a one position I rather have N/? than the reference nucleotide. Is there a way to do this ? Thanks

# Answer
> 2 votes
vSNP outputs a VCF file that includes positions with no coverage, represented as "N". Then there is a tool to merge a VCF file into the FASTA reference `vsnp_merge_vcf_into_fasta.py`. It will merge those "N" regions into the FASTA reference. There are options for considering frequency and quality. It's limited to just SNPs and those regions along the reference with no coverage. It does not merge small indels or unmapped reads. Pilon is a better tool for applying indels, but I'm not aware of it letting regions with no coverage be represented as "N".

```
vSNP_step1.py -r1 *_R1*fastq.gz -r2 *_R2*fastq.gz -r ref.fasta

```

Use `*zc.vcf` output to merge into FASTA

```
vsnp_merge_vcf_into_fasta.py -f ref.fasta -v *_zc.vcf

```

# Answer
> 1 votes
This is similar to this question, which links to this solution.

The difference here is that you want to go through and replace low coverage (maybe zero-coverage?) positions with Ns rather than using the reference. You can likely accomplish this using `samtools depth` to measure coverage and then writing a script that will go through each position in your reference, check its coverage, and string-replace with an N. Probably a bash one-liner could do this if someone very devoted to bash one-liners came up with one (I don't care for them, but other people are more creative than I). Very quick and dirty pseudocode might look like this:

```
### Someone could come up with something better if they spent >3min on it.
### should be easy to get it working if you use biopython
# read_and_preprocess() stands for whatever you need to do to get data in the right shape
coverage_data = read_and_preprocess('samtools_coverage.txt')
# reference_with_variants.fasta comes from the solution linked above
reference_fasta = read_and_preprocess('reference_with_variants.fasta')

threshold = 0
fasta_string = ''
for position in reference_fasta:
  if coverage_data[position] <= threshold:
    fasta_string += "N"
  else:
    fasta_string += reference_fasta[position]

# nonexistent fn that writes a fasta.
write_as_fasta(fasta_string)

```

I am a little skeptical of this approach because it seems like **if you do not trust the reference sequence to stand in for zero-coverage areas, then you probably shouldn't trust it to align to and call variants against**. Read mapping can be a rather noisy process, especially when you are using a more distantly related reference. I think this is related to the @swbarnes2 comment.

I would suggest instead indicating the regions of low/zero coverage by putting them in lowercase `acgt` with your covered regions uppercase `ACGT`, which preserves the coverage information that you care about but still allows you access to the reference information. It's not a perfect solution, but I think it makes the most out of the marginal data you're working with, and it is trivial to change the above example code to do it.

Hope that helps-

# Answer
> 0 votes
Have you googled *de novo assembly*?

---
Tags: sequence-alignment, sam
---


title: Index a BAM file using pysam

text:
# Question
Title: Index a BAM file using pysam
(How) can you index a BAM file using `pysam`?

When I tried the intuitive `pysam.index` I got:

```
import pysam
my_bam = pysam.AlignmentFile("regular_bwamem_mapping.bam", "rb")
pysam.index(my_bam)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-5-a0b3a04ecb2c> in <module>
----> 1 pysam.index(L_bam)

~/miniconda3/lib/python3.7/site-packages/pysam/utils.py in __call__(self, *args, **kwargs)
     59             args,
     60             catch_stdout=kwargs.get("catch_stdout", True),
---> 61             save_stdout=kwargs.get("save_stdout", None))
     62
     63         if kwargs.get("split_lines", False):

pysam/libcutils.pyx in pysam.libcutils._pysam_dispatch()

TypeError: object of type 'pysam.libcalignmentfile.AlignmentFile' has no len()

```

# Answer
Oh you silly sausage, `pysam.index` takes a bam file name, not a python object.

```
import pysam
pysam.index("regular_bwamem_mapping.bam")

```

will index your `.bam` file.
> 11 votes

---
Tags: sam, indexing
---


# Question
Title: Variant calling for a subset of genes using whole gneome sequencing data
I have few 100 raw fastq files from whole-genome sequencing data and I would like to map these files to a set of genes only (and not whole genome) so as to find SNP's associated with them. Can anyone tell me what snp analysis pipeline is best for targeted variation and not whole genome? I have this in mind: Mapping the reads with bwa and then using varscan Hope to hear for more options

# Answer
Your approach is not recommended. Always align against the entire genome. If you align against a subset you might accumulate false mappings. The aligner will always try to find the best match for every read. If the true origin of the read is not in the reference then the aligner will still try to find an acceptable mapping position. Therefore, align against the full reference and then extract the regions you are interested in, e.g. with the region option of `samtools view`:

`samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]`

Check its manual for details.

While there is basically nothign wrong with `VarScan` it is quite old and not maintained anymore. I suggest you give `bcftools` a try (check its manual). Other maintained alternatives can be (among others) `freebayes`, strelka2, VarDict or the GATK from the Broad institute.
> 3 votes

# Answer
vSNP can be used to target specific regions of a genome. Pipeline is a 2 step process:

1. create VCF files (via Freebayes)
2. then in your case run all 100 VCF files to create SNP table

A filter file can be supplied to step 2. By filtering all position ranges other than the genes of interest, step 2 will output a SNP table (.xlsx) with only your genes of interest. Also, table will include annotation if GBK is provide.
> 0 votes

---
Tags: snp, variant-calling, bwa, wgs
---


# Question
Title: Why do these NCBI representative genomes for ape species have no Y chromosome?
I recently downloaded a genome from NCBI for chimpanzee and was surprised to see no Y chromosome. I added in one from a separate accession and carried on, but then today went to do the same for gorilla and saw the same thing, no Y.

Why would this be? I can see a separate entry for a Y chromosome for that same exact chimp used for the representative genome (Clint) here, but I don't see why it's not in the representative genome listing itself:

https://www.ncbi.nlm.nih.gov/nuccore/DP000054.3

# Answer
tl;dr: technical difficulties, or sex

For the gorilla, that genome is from a female gorilla (Kamilah the gorilla) so she doesn't have a Y chromosome.

For the chimp link, there is a visualization of the chromosomes that you can click on that takes you to the genome browser for the Y chromosome (see image). Notably, there is a more recent genome build that does have a Y chromosome sequence assembly that is directly linked, they just haven't updated the main reference sequence yet.

To directly address the question, the Y chromosome is very hard to assemble. The easy part of it to assemble is already on the X chromosome anyways (the pseudoautosomal region), and the hard part of it is often a mess of repeats, especially in mammals. Many groups will decide to just assemble the X as well as they can and call it good, due to the difficulties associated with assembling the Y. If they are working with a male individual the Y sequence is still there, it just may not be anchored to a chromosome sequence, so it's only in the unscaffolded/unanchored assembly ("master WGS" I believe, usually). It is also probably highly fragmented (it isn't a whole chromosome of sequence stitched together, but rather a lot of little pieces of sequence).

In some cases they will focus all sequencing efforts on female individuals specifically to make their assembly job easier, because in females the X is at the same copy number as the autosomes and doesn't have interference from homologous regions on the Y.
> 8 votes

---
Tags: reference-genome, phylogenetics
---


# Question
Title: Calculate molecular mass from chemical formula C9H15N4O8P with Python?
Can I calculate the molecular mass from a formula like this `C9H15N4O8P` with Python e.g. using `rdkit`? I think `rdkit.Chem.Descriptors.MolWt` needs a molecule object. Though, to generate that I would need e.g. a SMILES string. Is there any other python module that can do the trick?

I found this: http://verahill.blogspot.com/2012/10/molecular-weight-calculator-in-python.html

# Answer
I tried to adapt this approach within rdKit to calculate exact mass (for mass spectrometry), but was not successful.

This library (https://pypi.org/project/molmass/) appears to be able to do this:

```
from molmass import Formula
f = Formula('C7H11N3O2')
f.formula
f.isotope.mass
'169.08512660696'

```
> 3 votes

# Answer
As you mentioned, rdKit does not accept a chemical formula as input. The main reason for this are probably isomers. You could however use rdKit to initiate all the atoms, get their mass and add them together.

```
import re
from rdkit import Chem

def get_mass(formula):

    parts = re.findall("[A-Z][a-z]?|[0-9]+", formula)
    mass = 0

    for index in range(len(parts)):
        if parts[index].isnumeric():
            continue

        atom = Chem.Atom(parts[index])
        multiplier = int(parts[index + 1]) if len(parts) > index + 1 and parts[index + 1].isnumeric() else 1
        mass += atom.GetMass() * multiplier

    return mass

```

Then...

```
>>> get_mass("C9H15N4O8P")
338.21299999999997

```
> 3 votes

---
Tags: python
---


# Question
Title: How to calculate Centimorgan lengths from run lengths in one-to-one dna profile comparison
For easier management of a database of dna profiles, I'm creating a tool to perform one-to-one matching on two profiles with snp genotype information.

Basically, something that works similar to GEDMatch's one-to-one matching feature.

Determining the length of runs of half or fully matched snp's, including some allowance for a genotyping error rate, is the easy part.

My question is about converting an snp run length to centimorgans.

I figure that the length needs to be adjusted by the likelihood of a change in the particular region of the chromosome.

Since my dna profiles and the positions of snps are based on Genome Reference Consortium Human Build 37, I figure that his information could be in the reference build, directly or derivable in some way.

Looking at the GRCh37 download files and the definitions of their formats, I'm unable to find such information.

Is this information present in GRCh37, or can it be downloaded from another reference source, and from which field in the dataset can it be retrieved or derived?

# Answer
Try having a look at this link:

http://bochet.gcc.biostat.washington.edu/beagle/genetic\_maps/plink.GRCh37.map.zip

If you unzip the file, there is a file for each chromosome (1:22) and then 3 for the x chromosome. The first column of each file details which chromosome the file corresponds to, the 3rd column gives you the cumulative genetic distance in cM and the 4th column gives you the genomic coordinate (position) in bp.

You can then calculate the distance in cM between two positions from this data by subtracting the cM positions between the two positions you are interested in.
> 3 votes

---
Tags: snp, phylogenetics
---


title: Edit all the fasta headers using awk

text:
# Question
Title: Edit all the fasta headers using awk
kindly your comment highly appreciated

have more than 4k header sequences look like:

```
>LTR_retrotransposon100_Gypsy?
cgtgcaccccaatgttcatagcagcactatttacaatagccaagacatggaaacaaccta aatgtcctttgacaggtgactggataaagaagctgtggtatatttatacaatggaatgct actcagccataaaaaagaataagataatgccatttacagcaacatggatagacctggaga atgttattctaagtgaagtaagccagaaagaaaaagaaaaataccatatgctatcactta

>LTR_retrotransposon10011_ERVK
taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>LTR_retrotransposon10011_ERVL-MaLR
    taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>LTR_retrotransposon100_Copia
    cgtgcaccccaatgttcatagcagcactatttacaatagccaagacatggaaacaaccta aatgtcctttgacaggtgactggataaagaagctgtggtatatttatacaatggaatgct actcagccataaaaaagaataagataatgccatttacagcaacatggatagacctggaga atgttattctaagtgaagtaagccagaaagaaaaagaaaaataccatatgctatcactta

```

I would like to use awk or sed to remove everything between the underscore, including the underscore and replace it with / plus add sequential RNA followed by #.

and the target output could be:

```
>RNA1#LTR/Gypsy?
ctcagcagcactatttacaatagccaagacatggaaacaacctaaatgtcttatcaatag atgactggataaaggagctgtggtatatctatacaatggaataccatcagccataaaaaa gaataaaatattgccatttgcagcaacatggatggacctggagattatcattctaaggga agtaagccagaaagagaaagaaaaataccatatgatatcacttatatgtggaggtaaaaa aaaaaaaaaagacacaaatcaatttatttgcaaaacatacaTGGActttcagacatagaa

>RNA2#LTR/ERVL-MaLR
taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>RNA3#LTR/ERVK
    taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>RNA4#LTR/Copia
    cgtgcaccccaatgttcatagcagcactatttacaatagccaagacatggaaacaaccta aatgtcctttgacaggtgactggataaagaagctgtggtatatttatacaatggaatgct actcagccataaaaaagaataagataatgccatttacagcaacatggatagacctggaga atgttattctaagtgaagtaagccagaaagaaaaagaaaaataccatatgctatcactta

```

# Answer
**EDIT after OP's update:**

Try this Perl one-liner:

```
perl -pe 'BEGIN { $i = 1 } $i += s{>([^_]+)_.*_}{>RNA${i}#${1}/}' input_file > output_file

```

Here, the command line flags are:
`-e`: tells the perl interpreter to look for the code inline, rather than in a file with the script,
`-p`: loop over the input one line at a time, execute the code in the one-liner for each iteration of the loop, and `print` the result into `STDOUT`.
The code has regex replacement (see `perlre` for more details):
`s{>([^_]+)_.*_}{>RNA${i}#${1}/}`:
replace an underscore, followed by zero or more occurrences (`*` quantifier) of any character (`.`), followed by an underscore, with slash (`/`).
Keep `>` and the original header part before the underscore, captured by parens into `$1` variable.
Add string `RNA${i}`, where `$i` is the counter of the fasta headers. The counter is incremented by 1 when the substitution occurs, that is, when the next header occurs.
The variables are used like so in the substitution: `${i}`, `${1}` to prevent wrong interpolation.
I am also using `s{}{}` as delimiters for readability (as in most cases) instead of the default `s///`.

**EXAMPLE:**

```
echo ">LTR_retrotransposon100_Gypsy?\nACGT\n>LTR_retrotransposon10011_ERVK\nTCGA" | perl -pe 'BEGIN { $i = 1 } $i += s{>([^_]+)_.*_}{>RNA${i}#${1}/}'

```

Prints:

```
>RNA1#LTR/Gypsy?
ACGT
>RNA2#LTR/ERVK
TCGA

```
> 4 votes

# Answer
Your own answer doesn't seem to do what you ask as it doesn't add the sequential identifier to each sequence. Maybe try this?

```
awk -v rec=0 '{
    if($0 ~ "^>"){
        rec++
        sub(">", ">RNA" rec "#", $0)
        sub("_.*_", "/", $0);
    }
    print $0
}' yourfile.fa

```
> 3 votes

# Answer
```
sed -e 's:^\(.*\)_.*_\(.*\)$:\1/\2:' < input > output

```
> 2 votes

---
Tags: fasta, phylogenetics, filtering, data-preprocessing, awk
---


# Question
Title: Why my bim file doesn't match to my ped file as the Plink documentation suggests?
Plink documentation about .bim says the columns 5 and 6 are the Allele 1 and Allele 2 respectively. The documentation about the .ped in turn, says "...The first six fields are the same as those in a .fam file. The seventh and eighth fields are allele calls for the first variant in the .map file ('0' = no call); the 9th and 10th are allele calls for the second variant; and so on.."

I've simulated a dataset as per:

```
$ plink1.9 --simulate ds1.sim acgt --out ds1
...
$ head ds1.bim
1       null_0  0       1       C       T
1       null_1  0       2       A       C
1       null_2  0       3       C       A
1       null_3  0       4       T       A
1       null_4  0       5       G       C
1       null_5  0       6       T       C
1       null_6  0       7       T       C
1       null_7  0       8       G       C
1       null_8  0       9       C       G
1       null_9  0       10      T       A

```

Than I converted it to plain PED and MAP:

```
$ plink1.9 --bfile ds1 --recode --out ds1
...
$ head ds1.ped > ds1.head.ped
$ vim ds1.head.ped
per0 per0 0 0 2 2 T T A C A A T A C C T C T C C C....

```

Given the data above, why my .bim file shows "C and T" for the first variant while my .ped file shows "T and T"? Apparently they don't match as per described in the documentation (or more likely, I'm misinterpreting it)...

Also, how do those 2 files related to each other? E.g.: How do I know that my Sample 1 has the SNP's null\_3, null\_4 or null\_5, etc?

I need to clarify those doubt because I want to fit a LightGBM model with the Genomic Data (can be plink or vcf). Not sure yet how to do. I was trying to convert the plink files to a CSV-style somehow, to be able to use in my model. Not sure if this is the way to go though...

Thanks!

# Answer
> 2 votes
Responded in the plink2-users group (https://groups.google.com/forum/#!topic/plink2-users/zgJxdXxvdLo ).

---
Tags: machine-learning, plink
---


title: Using cellranger for non-10x data

text:
# Question
Title: Using cellranger for non-10x data
I'm trying to use cellranger to `mkfastq` and then `count` and `aggregate` single-cell data. This data was not generated using 10x Genomics RNASeq.

When I run `count`, cellranger fails to auto id the chemistry and crashes out. Is it possible to wrangle the data into the cellranger method or should I look elsewhere. I'd ideally like to use the same method for all my single-cell data but not all of it is 10x.

# Answer
Possible? Sure. Your biggest problem is that 10x is expecting barcodes from its whitelist. I used to alter the python scripts to accept a 'new' kind of chemistry where the barcodes and umi lengths matched what I had, but I found that my modifications stopped working when the version updated, and I couldn't figure it out, so I started modifying fastqs instead.
> 1 votes

---
Tags: single-cell, cellranger
---


# Question
Title: How to generate ultrametric phylogenetic trees?
*This question has also been asked on Biostars*

I am going to use ultrametric tree for downstream analysis. I would like to see use the ultrametric tree for Ranger-DTL for infer gene family evolution by speciation, gene duplication, horizontal gene transfer, and gene loss.

I have generated rooted species tree with 1000 bootstraps by RAxML. Further, I need to generate ultrametric tree for the same. Therefore, Please suggest me how to generate ultrametric tree from the rooted species tree.

Thank you.

# Answer
An ultrametric tree is simply a global molecular clock and can be easily generate in PAML using a single evolutionary rate. You take you RAxML tree in phylip format (PAML might do nexus format know), and click the "clock function". PAML will do the rest.

You can alternatively generate it de novo via UPGMA neighbor-joining tree which is found in the MEGAX package (which is probably easier and considered a technical definition of an ultrametric tree).

A global molecular clock (as opposed to a local molecular clock) is quite a specific request, but if you performed a likelihood ratio test, it would identify "clock-like" behaviour.
> 1 votes

---
Tags: python, phylogenetics, phylogeny, linux, perl
---


# Question
Title: select unique names while removing the duplicates from a column
I have a file with Ensembl gene and Transcript IDs :

>

I want to retain only the first entires from both the column so that the output looks like this:

```
ENSG00000003137 ENST00000001146
ENSG00000003402 ENST00000309955

```

and ignore the rest. Have tried awk `awk '!a[$1$2]++'` and versions of it. Nothing works. Kindly suggest some way to get the desired output.

regards.

# Answer
> 1 votes
Keep it simple:

```
awk '!a[$1]++' file.txt

```

If you have a different command line separator you need to specify with the `-F` option.

Space

```
awk -F' ' '!a[$1]++' text.txt

```

Tab

```
awk -F'\t' '!a[$1]++' text.txt

```

If you don't have the first column in order maybe should be a good idea sort first.

```
sort text.txt | awk -F' ' '!a[$1]++'

```

The name of the array in this case is `a`. There is a better explanation in this post.

# Answer
> 1 votes
That sound like a job for a perl one-liner:

```
perl -ane 'unless(exists $seen{$F[0]}){print $_; $seen{$F[0]}++}'

```

Explanation:
The input is read line-wise and printed back to the command line if the ID in the first column has not been seen before.

* `perl -ane '` reads input line-wise and splits fields to array @F
* `unless(exists $seen{$F[0]}){` checks if ID exists in hash %seen and continues if that is not the case
* `print $_;` prints the original input to STDOUT
* `$seen{$F[0]}++}'` adds the ID to the hash %seen

---
Tags: ensembl
---


# Question
Title: How to create a .phy file for constructing a phylogenetic tree?
I'm new here. I need to construct a phylogenetic tree from 16s Sequence data (merged 151 bp) and color the branches/nodes based on my experimental conditions. From this blog, it looks like I need a .phy file to make the tree. I want to use Biopython for this since I'm reasonably comfortable with python.

I have a csv file with the 16s sequences, the phylogenetic information, and the "name" of the sequence, and the experimental conditions. How can I make a .phy file from this?

I would appreciate any help. Thank you.

# Answer
> 1 votes
You simply use Biopython's SeqIO and AlignIO. The documentation is https://biopython.org/wiki/AlignIO etc ... The format you seek is phylip

# Answer
> 1 votes
Python works great for tying everything together and file manipulations. I find stepping out of Python, using system calls, best for alignment and tree building. For example:

```
os.system("mafft --auto {} > {}" .format("unaligned_16s_seq_in_file.fasta", "alignment_out_file.fasta"))

os.system("raxml -f a -s {} -p 12345 -x 12345 -# 100 -m GTRCAT -n {}" .format("alignment_out_file.fasta", "tree_out_file.tre"))

```

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


# Question
Title: Q: What analyses can I perform with a completely phased genome assembly?
Phasing a genome is the process of determining which variants lie on which copy of each chromosome. For example, position 5,430,500 might have an A on the paternal chromosome, but a T base on the maternal chromosome.

A haplotype block is a section of the genome for which we know which haplotype has what variants. A completely phased genome has phase blocks that are the full length of each chromosome. Every variant on every chromosome can be assigned specifically to the maternal to or the paternal haplotype.

**My question is**, what can we learn from having completely phased genomes that we cannot from an unphased genome, or from a phased genome with short haplotype blocks?

# Answer
> 2 votes
Some population genetic inference methods require the data to be phased because they look for long stretches of matching alleles between individuals, rather than just matching genotypes. Here, 'genotype' data is un-phased, and 'allele' data is phased.

This is because there is a degree of ambiguity when it comes to matching heterozygous positions - if the data is un-phased, you don't know if 2 individuals share mutations on the same chromosome (here, I mean different copies of the same chromosome - like either copy of chromosome 1, for example). This is important, because the fact that the mutation is shared on the same chromosome means that they are much more likely to share recent common ancestry with one another.

Let's say we have genotype data for two individuals (where the number corresponds to the number of non-reference alleles present):

```
Ind 1

0 ----- 1 ----- 2 ----- 0

Ind 2

0 ----- 1 ----- 2 ----- 0

```

If were looking at unphased genotype data only, we know the alleles match exactly in the 1st, 3rd and 4th positions, because they are homozgyous (either 0 or 2 copies of the non-referencee allele), but because the second position is heterozygous, we can't be sure which chromosome the non-reference allele is on. This is because the haplotype configurations for individual 1 could be either

```
0 ----- 0 ----- 1 ----- 0 Chromosome1a
0 ----- 1 ----- 1 ----- 0 Chromosome1b

```

or

```
0 ----- 1 ----- 1 ----- 0 Chromosome1a
0 ----- 0 ----- 1 ----- 0 Chromosome1b

```

i.e., we don't know which chromosome the non-reference allele at position 2 is on.

If we then phased the data, and found that the non-reference allele was on chromosome 1a for Individual 1, and chromosome 1b for Individual 2 there is less evidence of recent shared ancestry than otherwise. I think this is more or less the approach that 23andme takes when they are inferring your relationship to other individuals in the dataset. After writing this, I realised there is a much better description of what I just said on the 23andme website!

For example, the above applies to the ChromoPainter/fineSTRUCTURE algorithms. This means that they are able to detect more subtle genetic differences between individuals than if you were to use non-phased data.

This is a slightly niche example, and there occasions where you need phased data, but this is the one I am most familiar with :). Off the top of my head, methods such as Multiply Sequential Markovian Coalescent (MSMC), RELATE, and Extended Haplotype Homozygosity all require the data to be phased.

# Answer
> 0 votes
A situation where having phasing data is useful is with allele-specific protein binding. Some individuals will have a heterozygous mutation within the binding site of a particular protein. If you perform a standard chromatin immunoprecipitation (ChIP) sequencing assay, you'll find sequencing reads around the protein's binding site, but you won't know which allele they'll come from since sequencing is an ensemble measurement.

However, if you have phasing information, you can measure whether this specific mutation has an effect on that protein's binding, and get some insight into the function of that mutation in the context of that individual's cells.

---
Tags: vcf, snp, phasing, variants
---


# Question
Title: Snakemake: Migrating from deprecated cluster.json to new profiles.yaml
I am an avid user of Snakemake. Recently we have been refreshing our pipelines and I saw that a cluster.json file is no longer the recommended way to store the cluster configuration.

I used to start my pipelines like this:

```
snakemake --cluster-config cluster.json --cluster "qsub -W umask=002 -d cluster_outputs  -l nodes={cluster.nodes}:ppn={cluster.ppn} -N {rule}" --jobs 245 --use-conda --rerun-incomplete --keep-going --latency-wait 60

```

With the cluster.json looking like this:

```
{
    "__default__" :
    {
        "cluster"   : "qsub",
        "jobs"      : 70 ,
        "nodes"     : 1,
        "ppn"       : 1
    },
    "rsem_quantify" :
    {
        "ppn"       : 8
    },
    "star_map_human" :
    {
        "ppn"       : 4
    },
}

```

Now, I understand that some of this goes in to the new profiles.yaml like this:

```
cluster: qsub
jobs: 70

```

But how do I set the ppn values per rule? And how do I set the other configuration options into profile.yaml (like: snakemake "qsub -W umask=002 -d cluster\_outputs -l nodes={cluster.nodes}:ppn={cluster.ppn} -N {rule}", etc.)

Any examples are appreciated, https://github.com/snakemake-profiles is not really helping

# Answer
It looks like you can just bundle your cluster config in its own YAML file, and reference that from within a profile YAML. Here's a clumsy little example using a fake job submission script: https://github.com/ressy/example-snakemake-profiles

---

The current documentation makes it sound like cluster configuration support is going away and being replaced by profiles:

> While still being possible, cluster configuration has been deprecated by the introduction of Profiles.

...but from those example profiles, that might be misleading; here's an example template for a profile YAML that references a separate cluster-config YAML:

https://github.com/Snakemake-Profiles/generic/blob/master/%7B%7Bcookiecutter.profile\_name%7D%7D/config.yaml

```
...
cluster-config: "$((INSTALDIR))/overwrite_cluster_param.yaml" #abs path
...

```

That cluster-config YAML has the same default and per-rule entries just your JSON version does:

https://github.com/Snakemake-Profiles/generic/blob/master/%7B%7Bcookiecutter.profile\_name%7D%7D/overwrite\_cluster\_param.yaml

```
...
# specific parameters for certain rules, which need more time/memory

#run_assembler:
#  queue: bigmem
#  time: 100
#   threads and memory definen in config file

```

So I think you can just set up a pair of YAML files accordingly including options for your rsem\_quantify and star\_map\_human rules. From my own test it looks like it handles the cluster argument templating just as it did before when using a profile. (If this is correct we should nag them to make the new documentation more explicit though.)
> 1 votes

---
Tags: snakemake
---


# Question
Title: How can i get data set of illness and not illness people?
I'm reading a virology paper where the authors made a classification for the patients as to whether they had a specific disease outcome or not against the underlying sequence ...

> The research focused on two different class of DNA which are Hepatitis B Virus (HBV) Genotype B nucleotide sequence and Non-HBV nucleotide sequences\`

* How can I get the data set or how can search to find it ?
* How can data sets of sequences date classify between illness and asymptomatic patients?

# Answer
The data set will be available on Genbank, you simply search for the "authors surname and HBV" and their entire data set should be revealed.

The HBV genome is 5Kb, the smallest viral genome so its not alot of data.

I am confused by "Genotype B nucleotide sequence and non-HBV nucleotide sequences". What I suspect is the investigators looked at HBV genotype B in comparison to HBV genotypes A, C to I. Essentially they wanted to assess the relationship between genetic variation and pathogenesis, which is an attempt to assess virus-induced cancer. Obviously I don't know the paper but I'm aware of the basic research front in this area. HBV has a very wierd evolution and I personally wouldn't want to investigate it.
> 2 votes

---
Tags: phylogenetics, phylogeny, covid-19, nucleotide-models, sars-cov-2
---


# Question
Title: how to export the monocle data into Seurat
I have tried to export the monocle data into Seurat and got the error

```
"Error in Seurat::CreateSeuratObject(raw.data = data, normalization.method = "LogNormalize",  :
  unused arguments (raw.data = data, normalization.method = "LogNormalize", do.scale = TRUE, do.center = TRUE, is.expr = monocle_cds@lowerDetectionLimit)
Calls: exportCDS
In addition: Warning message:
In if (export_to == "Seurat") { :
  the condition has length > 1 and only the first element will be used"

```

Here is my code `seurat_SAN_cds_subset <- exportCDS(SAN_cds_subset)`

Is anyone familiar with this?

# Answer
There's no default for the `export_to` option, so you just need to specify `export_to="Seurat"`. Note that monocle2 doesn't seem to be compatible with the newest Seurat (3.1) releases.
> 1 votes

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


# Question
Title: How I change a function not to need an input
There is a function does need an input as selected genomic site across genomic regions. But my genomic region bed file already is a subset of interesting site (like promoters). I want this function to work not complaining about genomic site However when I leave genomic site says

```
> prepare_elements_regions_tf <- prepare_elements("regions.bed",
+                                                 window_size = 50000)
Successfully read elements
 Error in utils::read.delim(sites_file, stringsAsFactors = FALSE, header = FALSE) :
  argument "sites_file" is missing, with no default

```

The whole function is here

https://github.com/wenmm/ActiveDriverWGS/blob/master/R/ActiveDriverWGS.R

But part of function that nags is this

```
prepare_elements = function(elements_file, sites_file, window_size = 50000, mc.cores = 1, recovery_dir = NA) {

    try_error = "try-error"

    element_coords_filename = paste0(recovery_dir, "/ADWGS_elements_coords_recovery_file.rsav")

    element_coords_1 = NULL

    load_result = suppressWarnings(try(load(element_coords_filename), silent = TRUE))

    if (class(load_result) == try_error) {

        element_coords_1 = prepare_element_coords_from_BED(elements_file)

        if (!is.na(recovery_dir)) {

            save(element_coords_1, file = element_coords_filename)

        }

        cat("Successfully read elements\n")

    } else {

        cat("Recovered elements\n")

    }

    element_coords = coords_sanity(element_coords_1)

    element_3n_filename = paste0(recovery_dir, "/ADWGS_element_3n_recovery_file.rsav")

    element_3n = NULL

    load_result = suppressWarnings(try(load(element_3n_filename), silent = TRUE))

    if (class(load_result) == try_error) {

        element_3n = add_trinucs_by_elements(get_3n(element_coords))

        if (!is.na(recovery_dir)) {

            save(element_3n, file = element_3n_filename)

        }

    }

    window_coords_filename = paste0(recovery_dir, "/ADWGS_window_coords_recovery_file.rsav")

    window_coords = NULL

    load_result = suppressWarnings(try(load(window_coords_filename), silent = TRUE))

    if (class(load_result) == try_error) {

        window_coords = coords_sanity(get_window_coords(element_coords, window_size))

        if (!is.na(recovery_dir)) {

            save(window_coords, file = window_coords_filename)

        }

    }

    window_3n_filename = paste0(recovery_dir, "/ADWGS_window_3n_recovery_file.rsav")

    window_3n = NULL

    load_result = suppressWarnings(try(load(window_3n_filename), silent = TRUE))

    if (class(load_result) == try_error) {

        window_3n = add_trinucs_by_elements(get_3n(window_coords))

        if (!is.na(recovery_dir)) {

            save(window_3n, file = window_3n_filename)

        }

    }

    raw_site_coords_filename = paste0(recovery_dir, "/ADWGS_raw_site_coords_recovery_file.rsav")

    site_coords_filename = paste0(recovery_dir, "/ADWGS_site_coords_recovery_file.rsav")

    raw_site_coords = NULL

    site_coords = NULL

    sites_load_result = suppressWarnings(try(load(site_coords_filename), silent = TRUE))

    if (class(sites_load_result) == try_error) {

        load_result = suppressWarnings(try(load(raw_site_coords_filename), silent = TRUE))

        if (class(load_result) == try_error) {

            raw_site_coords = utils::read.delim(sites_file, stringsAsFactors = FALSE, header = FALSE)

            colnames(raw_site_coords) = c("chr", "starts", "ends", "site_info")

            raw_site_coords$chr = gsub("chr", "", raw_site_coords$chr, ignore.case = TRUE)

            raw_site_coords$chr = paste0("chr", raw_site_coords$chr)

            if (!is.na(recovery_dir)) {

                save(raw_site_coords, file = raw_site_coords_filename)

            }

            cat("Successfully read active sites\n")

        } else {

            cat("Recovered active sites\n")

        }

        if (!is.null(raw_site_coords[[1]])) {

            # first only keep sites that overlap with element coordinates

            site_gr = GenomicRanges::GRanges(raw_site_coords$chr,

                IRanges::IRanges(start=raw_site_coords$starts, end=raw_site_coords$ends))

            element_gr = GenomicRanges::GRanges(element_coords$chr,

                IRanges::IRanges(start=element_coords$starts, end=element_coords$ends))

            rm(raw_site_coords)

            overlaps_found = GenomicRanges::findOverlaps(element_gr, site_gr)

            site_here_gr = site_gr[S4Vectors::subjectHits(overlaps_found)]

            rm(site_gr, element_gr, overlaps_found)

            site_coords = do.call(rbind, parallel::mclapply(1:nrow(element_coords),

                get_sites_per_element, element_coords, site_here_gr, mc.cores=mc.cores))

            rm(site_here_gr)

            site_coords = coords_sanity(site_coords)

        }

        if (!is.na(recovery_dir)) {

            save(site_coords, file = site_coords_filename)

        }

    } else {

        cat("Recovered active sites\n")

    }

    rm(sites_load_result, element_coords)

    site_3n_filename = paste0(recovery_dir, "/ADWGS_site_3n_recovery_file.rsav")

    site_3n = NULL

    load_result = suppressWarnings(try(load(site_3n_filename), silent = TRUE))

    if (class(load_result) == try_error) {

        site_3n = get_3n(site_coords)

        if (!is.na(recovery_dir)) {

            save(site_3n, file = site_3n_filename)

        }

    }

    list(element_coords=element_coords_1, element_3n=element_3n, window_coords=window_coords, window_3n=window_3n, site_coords=site_coords, site_3n=site_3n)

}

```

How I get rid of nagging part?

# Answer
Think about it this way - your `elements_file` already has only the sites of interest and hence does not need to be subject to a subset operation. In other words, if you were to have a `sites_file`, it would have the same content as the `elements_file`.

Can you see where I am going with this?
> 1 votes

---
Tags: r
---


# Question
Title: Is there a "definitive" database for Mobile Genetic Elements?
This isn't quite my speciality, but I've been tasked with finding any evidence of Mobile Genetic Elements (like transposable elements, retrotransposons, bacteriophages, etc.) within some WGS samples.

Is there a "definitive" database used with sequences? How "species-specific" would this be?

Furthermore, are there standard algorithms used to find these MGEs? I suspect we don't have a comprehensive database of all possible MGEs, but it may be possible to query against the sequences we do indeed have (or detect any repetitive elements within the genome).

# Answer
Would advise doing Denovo identification method structure-based. Then use these public databases ( RepBase, Dfam ) to look for high similarity or homology exploration, this possibility if your species, the repeats not yet discovered, or deposited in the databases.

If the species well-investigated earlier and the mobile element recognised, I assume it will be a straightforward task using species-specific search via Repeatmasker tool with the -species flag to obtain the MobileElement annotation in the genome.

Let me know if you need any data from Repbase related to your species, i have an entree.
> 2 votes

# Answer
Back in the Day™ RepeatMasker\[1\] was the go-to tool everyone used for TE annotation in eukaryotes. At the time, the RepBase database from GIRI\[2\] (license required) was the *de facto* database recommended for RepeatMasker users. I know *de facto* isn't necessarily *definitive*, but it may have been the only public database supported for some time.

It looks like more recent efforts have focused on developing the Dfam database, mirrored after the design of databases like Pfam and Rfam.

---

\[1\]Smit AFA, Hubley R, Green P (2013-2015) RepeatMasker Open-4.0, http://www.repeatmasker.org.

\[2\]Jurka J, Kapitonov VV, Pavlicek A, Klonowski P, Kohany O, Walichiewicz J (2005) Repbase Update, a database of eukaryotic repetitive elements. *Cytogentic and Genome Research* 110:462-467, doi:10.1159/000084979.
> 1 votes

---
Tags: blast, repeat-elements, transposable-elements, bacteria
---


# Question
Title: Can I get longer contigs by changing MEGAHIT settings?
My current settings are `megahit -r input.fq --num-cpu-threads 32 --min-contig-len 300 --presets meta-large -o output`. I picked the 'presets meta-large' because I am dealing with a quite large and diverse metagenome out of which I only need some genes.

I've got a contig of interest ~19kb, with a depth of 16. Could I make this contig larger with more relaxed settings? Help would be appreciated!

# Answer
> 1 votes
You can use the option `--k-step 10` and `--kmin-1pass` for low-coverage datasets.

---
Tags: assembly, metagenome
---


title: STAR quantMode vs RSEM vs Kallisto

text:
# Question
Title: STAR quantMode vs RSEM vs Kallisto
I recently discovered this Snakemake pipeline for RNASeq that uses STAR's quantMode to quantify gene expression for DESeq2 differential expression analysis.

In the past I've always seen the workflow as:

1. STAR alignment
2. RSEM/HTSeqCount quantification
3. DESeq2 analysis

More recently I've also seen:

1. Kallisto/Salmon pseudo-alignment quantification
2. DESeq2/Sleuth analysis

I tried searching online for differences between STAR quantMode vs other quantification algorithms but couldn't find many details. What are the benefits and drawbacks of using STAR quantMode vs RSEM/Kallisto/Salmon?

# Answer
> 3 votes
STAR quantMode (GeneCounts) essentially provides the same output as HTSeq-Count would, ie. number of reads that cover a given gene. This is the most simple measure of expression you could get from RNA-seq data. Kallisto and Salmon utilize pseudo-alignment to determine expression measures of transcripts (as opposed to genes). RSEM uses some algorithm to determine isoform fractions, but still based on alignments. (Note that I don't have any personal experience using RSEM.)

It ultimately depends on the research question which method will be most appropriate for you.

# read counts

pros:

* relatively easy to understand concept-wise
* can be used to generate many different expression measures (FPKMs, TPM, CPM, etc.) and is, thus, more versatile in down-stream analyses

cons:

* you can't distinguish isoforms (from the counts themselves)
* requires alignment and, thus, more time to obtain

# salmon/kallisto:

pros:

* faster (due to lack of alignment step)
* quantifies (known) isoforms

cons:

* expression values are more abstract, less easy to explain
* gene-level analyses will require some extra work

# RSEM

Likely lies somewhere in between, still requiring alignments will make it slower. But it does give you expression values for isoforms.

# Answer
> 0 votes
A shorter version...STAR's is convenient, in that it gets done right alongside the alignment, but RSEM and Kallisto are smarter about dealing with ambiguous reads. So I'd say your results from either of those are superior.

Personally, I use STAR for alignment, then give the bam to RSEM. STAR's gene count is good for a reality check to make sure that I got the strandedness right, and RSEM's gene level expected counts are good for DESeq once they are rounded.

---
Tags: rna-seq, differential-expression, deseq2, rsem, kallisto
---


# Question
Title: Protein ligand docking: how to convert <protein>.pdb to <protein>.maps.fld?
Hello I'm helping to develop a cloud docking tool for screening compounds, similar to Swissdock but with mass throughput and GPU optimizations.

Specifically helping screen existing drugs against coronavirus proteins (couple structures solved, rest TASSER).

I'm planning to use GPU optimized AutoDock-GPU , which takes in .maps.fld

1. I know you can use autogrid to select the bounding box and generate the .maps.fld, but I've been unable to figure out the workflow.
2. for preliminary screening I want to search the whole protein without specifying a bounding box.
3. Is there a script for converting protein.pdb to .maps.fld?

Collaborators welcome.

Thanks!

# Answer
> 2 votes
There are a few steps in between `pdb` and `.maps.fld`.

Here is the list of several scripts that can do tasks for you that you downloaded with autodock MGLTools: http://autodock.scripps.edu/faqs-help/faq/where-can-i-find-the-python-scripts-for-preparing-and-analysing-autodock-dockings. Look at the `prepare_` files. Also note that the scripts come with a pre-bundled python, but you can install these scripts with conda. Here are the basic order of things:

```
prepare_receptor4.py -r protein.pdb
prepare_ligand4.py -l ligand.mol2
prepare_gpf4.py -l ligand.pdbqt -r protein.pdbqt -y
autogrid4 -p protein.gpf
prepare_dpf4.py -l ligand.pdbqt -r protein.pdbqt

```

## partial charge

Your PDB needs partial charges so first you convert it to `PDBQT` file. This is true for your protein and your ligand. That is no complicated parameterisation for your ligand —which is not necessarily a good thing! In this step it is essential to correct protonation if absent each have their own flags.

## grid

The you make a grid file (`gpf`), which contains your grid. To make a box that spans the whole protein simply use the unit cell dimensions form the PDB. For an mmCIF, if you open the header dictionary you need `_pdbx_struct_oper_list.matrix` using `MMCIF2Dict` from biopython's PDB submodule. For a PDB file, just search for the `CRYST1` line (cf. format).

Once you have have the grid you can run autogrid. This step will create the `map.fld`

### Why a big box is bad for a VS

However, all those who take this approach opt for a box that is huge: as there is no solvent having a huge bounding has little penalty. However, this is considered a very bad strategy for a virtual screen. This not because it uses up computer resources as these are dirty cheap. Structural knowledge is important and targeting the active site of an enzyme stops it, while binding to the surface does nothing unless it's an interface. So manually reading what the stuff are and choosing the box wisely will save you a lot of time downstream in analysis. There is after all the saying "one week in the lab will save you an hour in the library"...

## Protein and ligand

Lastly, merge the protein with the small molecule in a `dpf` file.

## Caveat against docking with models

It is not a good idea dock against models. Swissmodel (unless virtually identical), I-Tasser, Phyre or EVFold models. Docking is very sensitive to small structural difference that may have resulted from the modelling, so docking to models is highly discouraged. For Coronavirus protein, dock against SARS, which is highly similar. Or at the very least pay strong attention to the I-Tasser C-scores and discard the bottom 2/3 of their models.

## Note about coronavirus

Your objective of docking coronavirus may be a bit too late however. I did a wee interactive summary of the literature about the solved structure which is the protease. Proteases are very easy to rationally design drugs for. There are already:

* The Covid Moonshot project is an open collaboration between a fragment-screening X-ray facility (Diamond XChem) and many researchers from various disciplines (myself included).
* a paper that has an empirically validated coronavirus specific ligand
* many many manuscripts about virtual screens in bioarxiv/chemarxiv.
* licensed HIV protease inhibitors lopinavir and ritonavir that were initially suspected to be effective —a clinical trial

Also you may want to check out Galaxy project, a tool to run pipelines for genomics, but have started moving into biochemistry from genomics. Eg. https://covid19.galaxyproject.org/cheminformatics/.

---
Tags: protein-structure, covid-19, sars-cov-2, docking, structural-biology
---


# Question
Title: How to initiate RNA-Seq analysis of TCGA files?
I want to compare RNA-Seq datasets obtained from the TCGA to investigate how my gene of interest is implicated in different types of cancer.

I'm trying to download the data from the GDC Data Portal (https://portal.gdc.cancer.gov), but it is not user-friendly.

Can somebody give me some initial directions? I found a bioconductor package that facilitates access to the TCGA (https://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html), but first of all, I would like to understand how the files are organized, the nomenclature and how to inspect them.

thank you very much, Fabiano

# Answer
> 3 votes
The GDC API user guide (https://docs.gdc.cancer.gov/API/PDF/API\_UG.pdf) provides some guidance on the structure of GDC data. So for getting RNA-Seq data from GDC, you would provide this filter to get the 551 text files with workflow\_type set as "HTSeq - FPKM-UQ" and "TCGA-LUSC" for the project. There are other workflow types as well:

```
cases.project.program.name in ["TCGA"]
and cases.project.project_id in ["TCGA-LUSC"]
and files.access in ["open"]
and files.analysis.workflow_type in ["HTSeq - FPKM-UQ"]
and files.data_format in ["txt"]
and files.data_type in ["Gene Expression Quantification"]
and files.experimental_strategy in ["RNA-Seq"]

```

(You don’t need all those filters to get the 551 files; I included them to show what is available.) The txt files are returned as gzipped tar files that hold text files with rows of the form:

ENSG00000158486.12 5341.21601153

The above filters provide all the FPKM-UQ files for TCGA-LUSC. If you wanted to just get the data for normal samples as controls, that would require filtering on sample type as well, adding in the additional filter:

```
and cases.samples.sample_type in [“solid tissue normal”]

```

That will return 49 files. If you instead wished to get only the tumor files, the additional filter would instead be:

```
and cases.samples.sample_type in [“primary tumor”]

```

and that returns 502 files. Note that there are several different sample types for tumor and normal tissues; those just happen to be the ones used here. (Note that if you are using the GDC Repository portal, adding that additional filter requires you to select “Add a Case/Biospecimen Filter” on the “Cases” filter tab.)

Since GDC is part of the NCI Cancer Research Data Commons (CRDC), instead of downloading the data, you can actually work with it in the cloud using the various CRDC Cloud Resources such as the ISB-CGC: https://www.isb-cgc.org/. In fact, all the RNA-Seq data for TCGA is in an ISB-CGC Google BigQuery table: (https://console.cloud.google.com/bigquery?p=isb-cgc&d=TCGA\_hg38\_data\_v0&t=RNAseq\_Gene\_Expression&page=table) And the above count data can be pulled using the query:

```
SELECT HTSeq__FPKM_UQ
FROM `isb-cgc.TCGA_hg38_data_v0.RNAseq_Gene_Expression`
WHERE Ensembl_gene_id_v = "ENSG00000158486.12" AND project_short_name = "TCGA-LUSC"

```

If you want to discriminate between normal and tumor expression levels, the ISB-CGC BQ tables make it pretty easy to tag the expression values with the type of tissue involved. The sample type can be assigned by joining the previous query with the biospecimen table:

```
WITH bio as (
SELECT sample_barcode, sample_type_name
FROM `isb-cgc.TCGA_bioclin_v0.Biospecimen`
),
expr as (SELECT sample_barcode, HTSeq__FPKM_UQ
FROM `isb-cgc.TCGA_hg38_data_v0.RNAseq_Gene_Expression`
WHERE Ensembl_gene_id_v = "ENSG00000158486.12" AND project_short_name = "TCGA-LUSC")
SELECT bio.sample_type_name, expr.HTSeq__FPKM_UQ
FROM expr
JOIN bio
ON bio.sample_barcode = expr.sample_barcode

```

The GDC has listed some of the available sample types on this page: https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes Also, BigQuery has a user interface that makes exploring data sets and writing queries easy, but you need to log in with a Google ID, so the ISB-CGC has a BigQuery Table Search (https://isb-cgc.appspot.com/bq\_meta\_search/) that allows you to explore the available tables without logging in. If you want to run queries, Google offers a free tier that you can use to try things out and run small jobs (https://cloud.google.com/free/) and ISB-CGC offers free credits too. ISB-CGC also has examples of working in R with BigQuery within the Community Notebook Repository (https://isb-cancer-genomics-cloud.readthedocs.io/en/latest/sections/HowTos.html). A few useful notebooks would be 'How do I create cohorts of patients?', 'How do I plot a heatmap using data in BigQuery?', and 'How do I get started fast?'.

# Answer
> 0 votes
I find UCSC Xena to be well-organized and easy to use: https://xenabrowser.net/datapages/

Plain text files that are clearly labeled for all TCGA data (that is not protected).

---
Tags: rna-seq
---