# Question
Title: How can I classify the 3 clades(S, G, V) of the coronavirus without using protein data?
On GISAID they classified the coronavirus using 4 clades(S, G, V, Other).
I downloaded around 1,000 complete genomes of the coronavirus from GISAID and I would like to classify each one as belonging to one of the 4 clades(S, G, V, Other).

On the top left of the image below is the classification of the 3 clades(S, G, V).
More specifically, the three strains were classified based on three specific variants(This answer was taken from my previous question):

1. Strain S, variant ORF8-L84S: a variant in the gene "ORF8" which changes the leucine (L) residue at position 84 of the gene's protein product to a serine (S).
2. Strain G, variant S-D614G: a variant in the gene "S" which changes the aspartic acid (D) residue at position 614 of the gene's protein product to a glycine (G).
3. Strain V, variant NS3-G251V: a variant in the gene "NS3" which changes the glycine (G) residue at position 251 of the gene's protein product to a valine (V).

My problem is that using this method of classification for the coronavirus genome I am seeking information on the protein products of the genes in the genome.

So to be more specific here are my concise questions:

1. How can I classify a genome as "Strain S" using only the complete genome assembly, without any information about its proteins?
2. How can I classify a genome as "Strain G" using only the complete genome assembly, without any information about its proteins?
3. How can I classify a genome as "Strain V" using only the complete genome assembly, without any information about its proteins?

This picture was taken from this PDF on GISAID

# Answer
> 8 votes
It seems to be explained right there in the image you posted:

So, the three strains were classified based on three specific variants:

1. Strain S, variant ORF8-L84S: a variant in the gene "ORF8" which changes the leucine (L) residue at position 84 of the gene's protein product to a serine (S).
2. Strain G, variant S-D614G: a variant in the gene "S" which changes the aspartic acid (D) residue at position 614 of the gene's protein product to a glycine (G).
3. Strain V, variant NS3-G251V: a variant in the gene "NS3" which changes the glycine (G) residue at position 251 of the gene's protein product to a valine (V).

# Answer
> 2 votes
The most important thing to say is we are grateful this information has been released. It is important to note these are very large genomes of around 9000 amino acids in total.

The clades are assigned to single amino acid mutations. The tree is built using nucleotide data therefore synonymous (silent) mutations are resulting in phylogenetic resolution and the authors have superimposed amino acid mutations over this tree structure. They have represented the tree as a circular phylogram.

So the reason why there is loads of branches (of a phylogenetic tree) of different sizes is because there are different numbers of silent (mostly 3rd codon wobble) mutations between the isolates. RNA viruses mutate rapidly. There are only 4+ amino acid changes across the entire epidemic according to this analysis.

**Nomenclature** The clades are named after the amino acid mutation that arises within that group. Its a bit weird, but its okay. Hence S clade means a serine mutation, G clade a glycine mutation etc.... The thing is that S also equals Spike, which is a very important protein (below), so it is confusing because the S clade refers to ORF8. Which gene the amino acid mutation is occurring in is more important than the specific mutation. Secondly the "other" clade is dying out - which is not that informative, but others are holding their course.

**Mutation scheme and protein** So @terdon is right

* L-\>S at position 84 of ORF8 - Hmmmm I am not hugely excited by that because it is a very weird gene and is is not a radical mutation (the amino acids are related) which is absent in some closely related coronaviruses altogether. In other words some closesly related viruses can remove this gene altogether and remain funcational.
* D-\>G (aspartic acid is D) in the Spike protein at position 614. That is an interesting mutation because S protein contains the cellular receptor and that is an unusual mutation. S protein is critical for vaccination.
* G-\>V in the NS3 ... I dunno because glycine to valine means it will be a bit more hydrophobic. NS3 is in the structural protein region of the genome (including spike, capid, membrane and another structural protein) but I'm not sure precisely what is does

**Epidemiological interpretation**

Epidemiological interpretation is difficult because it doesn't really explain the possible differences in mortality rates between different countries, for example USA appears to have a low mortality rate in comparison to Italy but there is lots of mixing. It is also not easy to read the tree because it is a circular phylogram and I am struggling to read each "clade" (which has a precise phylogenetic diffinition) I have to take their word for it. Secondly is the colour scheme the 'red' swamps every other colour and gives a disproportionate red colour. This is an old epidemiologists trick in heat maps, there precise colour scheme albeit that is relative, appears to convey a "danger" message to the reader, even though mathematically its simply a consequence of the colour scheme.

Fundamentally it is not really clear whether the mutations confer phenotypic changes and these are huge genomes for an RNA virus. However, the G clade needs further investigation.

**Summary** The battle continues as to whether data collection skews, social skews, age-prevelence skews, climatic skews result in variation of mortality rates between countries? Versus is their a genetic difference in mortality. On this data I'd go for the non-genetic factors.

---

The answer to your second question is simple, you just translate the triplet codon into amino acids and identify which amino acid is present at the above 3 positions. You then simply denote the clade within the taxon ID of your nucleotide data. MEGA X has an alignment editor that will allow this translation.

---
Tags: rna-seq, phylogenetics, phylogeny, covid-19, sars-cov-2
---


# Question
Title: How to identify the reference genome assembly used to create a .hic file?
I have a bunch of .hic files of different species. But, I don't know which reference genome assembly (e.g. GRCh37 or GRCh38 in the case of human genome) the genomic co-ordinates of the .hic file correspond to. So I wonder if there is any way to figure out the version of the reference genome assembly from a .hic file.

**Therefore, my question is**: Does .hic file contain metadata such as the version of reference genome assembly? If yes, then how could I access it?

(I could obtain the information I need by digging into the methods sections of each paper corresponding to each .hic file. But it would be much better if I could get that right-away from the .hic file.).

# Answer
Turns out HiC file does store such metadata.

`straw` tool from Aiden Lab contains a script (read\_hic\_header.py) that prints out the metadata from a given HiC file. Usage: `read_hic_header.py <hic file or URL> [verbose]`

In order to obtain the metadata in the form of a dictionary, one may consider using a function of the same tool that is currently available in my fork.

```
from straw import straw_module
straw_module.read_metadata(infile,verbose=False)

Args
    infile: str, path to the HiC file
    verbose: bool, print the metadata if True else return the metadata
Returns
    metadata: dict, the metadata is returned if verbose is False.
                Keys of the metadata:
                HiC version,
                Master index,
                Genome ID (str),
                Attribute dictionary (dict),
                Chromosomes (dict),
                chromosomes2length (dict),
                Base pair-delimited resolutions (list),
                Fragment-delimited resolutions (list).

```
> 1 votes

---
Tags: reference-genome, chromosome-capture, hi-c
---


title: Differential Expression

text:
# Question
Title: Differential Expression
I have genes differentially expressed between two groups (case and control). I would like to annotate them by classifying them according to their biological functions.

I work on parasites of Plasmodium Falciparum type Dd2 which are divided into two groups (case and control). I carried out a Differential expression analysis with cuffdiff and CummeRbund then I annotated my genes on plasmodb.org.

Example of id genes: pfDd2\_14005040

Now I would like to do a pathway analysis to see which biological pathways these genes are involved in.

StupidWolf, I resumed my annotation with the 3D7 strain of plasmoduim Falciparum then I created my geneList as you indicated but when I make the dotplot nothing is displayed

```
setwd("C:/cygwin64/home/DIANGO/EXCELL")
library(org.Pf.plasmo.db)
library(clusterProfiler)

gene <- c ("PF3D7_0936800","PF3D7_1478900","PF3D7_1009700","PF3D7_0508500")
geneList = keys(org.Pf.plasmo.db)
enr <- enrichGO(gene,universe=geneList,OrgDb='org.Pf.plasmo.db',keyType ="SYMBOL",ont="BP")
head(enr@result,readtable=TRUE)
dotplot(enr)

```

The output gives a white rectangle

```
                   ID                           Description GeneRatio  BgRatio     pvalue
GO:0006323 GO:0006323                         DNA packaging       1/2  21/2615 0.01599974
GO:0071103 GO:0071103               DNA conformation change       1/2  30/2615 0.02281728
GO:0051276 GO:0051276               chromosome organization       1/2  58/2615 0.04387582
GO:0006996 GO:0006996                organelle organization       1/2 144/2615 0.10712138
GO:0016043 GO:0016043       cellular component organization       1/2 203/2615 0.14925925
GO:0006464 GO:0006464 cellular protein modification process       1/2 265/2615 0.19244223
            p.adjust     qvalue        geneID Count
GO:0006323 0.1026777 0.07205456 PF3D7_0508500     1
GO:0071103 0.1026777 0.07205456 PF3D7_0508500     1
GO:0051276 0.1316275 0.09237015 PF3D7_0508500     1
GO:0006996 0.2203602 0.15463872 PF3D7_0508500     1
GO:0016043 0.2203602 0.15463872 PF3D7_0508500     1
GO:0006464 0.2203602 0.15463872 PF3D7_1009700     1
> dotplot(enr)

```

Once again I need your help.

# Answer
> 1 votes
I guess plasmodium doesn't have that many online tools dedicated to it so you can use the annotation packages from bioconductor and try it. I used clusterProfiler for analysis and org.Pf.plasmo.db for annotation below.

From this paper, grabbed table 1, a table of supposed genes involved in something:

```
genes = c("PF3D7_1302100", "PF3D7_0501300", "PF3D7_0302500*", "PF3D7_1030200",
"PF3D7_0302200*", "PF3D7_1327300", "PF3D7_0410000", "PF3D7_0104300",
"PF3D7_0420300", "PF3D7_0810000*", "PF3D7_0522300", "PF3D7_0935700*",
"PF3D7_0424700*", "PF3D7_1477700*", "PF3D7_0202000*", "PF3D7_1371600",
"PF3D7_1036300", "PF3D7_0422900", "PF3D7_0501400", "PF3D7_0221700",
"PF3D7_0301000", "PF3D7_1454100", "PF3D7_0831800", "PF3D7_0424900",
"PF3D7_1447500", "PF3D7_1477400", "PF3D7_0424800", "PF3D7_0823700",
"PF3D7_1477300", "PF3D7_0814100", "PF3D7_1372100", "PF3D7_0201500",
"PF3D7_0825300", "PF3D7_0201700", "PF3D7_1470000", "PF3D7_1401000",
"PF3D7_0724900", "PF3D7_0601700", "PF3D7_1410000", "PF3D7_0109300",
"PF3D7_1347300", "PF3D7_0201800", "PF3D7_1242100", "PF3D7_0403300",
"PF3D7_0114500", "PF3D7_1301200", "PF3D7_0931200", "PF3D7_0832200.1",
"PF3D7_0305900", "PF3D7_0608310", "PF3D7_0221800", "PF3D7_0511100",
"PF3D7_1001300", "PF3D7_1225900", "PF3D7_1444600", "PF3D7_1349900",
"PF3D7_1008300", "PF3D7_0919600", "PF3D7_1210300", "PF3D7_0929700",
"PF3D7_0701600", "PF3D7_0201600", "PF3D7_0801600", "PF3D7_1463600",
"PF3D7_0201900", "PF3D7_0905000", "PF3D7_1310000", "PF3D7_0315200",
"PF3D7_1139200", "PF3D7_1023500", "PF3D7_1111300", "PF3D7_0930700",
"PF3D7_1350800", "PF3D7_0412300", "PF3D7_0831900", "PF3D7_0406600",
"PF3D7_1448100", "PF3D7_0114300", "PF3D7_0526100", "PF3D7_1243200",
"PF3D7_0921100", "PF3D7_0301900", "PF3D7_0622700", "PF3D7_0811800",
"PF3D7_0114200", "PF3D7_0221500", "PF3D7_1129850", "PF3D7_1142200",
"PF3D7_0114100", "PF3D7_0512900", "PF3D7_1346200", "PF3D7_1369400",
"PF3D7_0314300", "PF3D7_0324100", "PF3D7_0107700", "PF3D7_1122000",
"PF3D7_0618800", "PF3D7_0601200", "PF3D7_1360300", "PF3D7_1318500",
"PF3D7_1111700", "PF3D7_0901900", "PF3D7_1313900", "PF3D7_1101700",
"PF3D7_1313200", "PF3D7_1039700", "PF3D7_0222100", "PF3D7_1371500",
"PF3D7_1344700", "PF3D7_1408300", "PF3D7_0929800", "PF3D7_0115000",
"PF3D7_1462900", "PF3D7_1352600", "PF3D7_0626200", "PF3D7_0831750",
"PF3D7_0602300", "PF3D7_0521100", "PF3D7_1452100", "PF3D7_0718900",
"PF3D7_0101300", "PF3D7_0815700", "PF3D7_0113600", "PF3D7_0937200",
"PF3D7_1108800", "PF3D7_0221100", "PF3D7_0921500", "PF3D7_1038300",
"PF3D7_0201400", "PF3D7_0822200", "PF3D7_0937000", "PF3D7_0301100"
)

```

We can for example using clusterProfiler:

```
library(org.Pf.plasmo.db)
library(clusterProfiler)

```

Besides the list of interest, you need to get a "universe" of all genes tested, in your situation it will be genes that you tested for differential gene expression, below I used all the genes annotated:

```
geneList = keys(org.Pf.plasmo.db)

enr <- enrichGO(genes,universe=geneList, OrgDb='org.Pf.plasmo.db',
keyType ="SYMBOL",ont="BP")

```

Above is a GO term enrichment, biological processes, you can look at the result like this:

```
head(enr@result,readtable=TRUE)
                       ID
    GO:0042000 GO:0042000
    GO:0044417 GO:0044417
    GO:0044766 GO:0044766
    GO:0051808 GO:0051808
    GO:0051836 GO:0051836
    GO:1902579 GO:1902579
                                                                                Description
    GO:0042000                                             translocation of peptides or proteins into host
    GO:0044417                                                        translocation of molecules into host
    GO:0044766                                                                    multi-organism transport
    GO:0051808 translocation of peptides or proteins into other organism involved in symbiotic interaction
    GO:0051836            translocation of molecules into other organism involved in symbiotic interaction
    GO:1902579                                                                 multi-organism localization

```

And you can plot it:

```
dotplot(enr)

```

You can read more about the analysis you can do at the package website

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


# Question
Title: Marking optical or PCR duplicates with picard vs. samtools flagstat
I am trying to identify the best library prep method for noninvasive prenatal test samples, to be sequenced with Illumina Novaseq.

One metric that I am evaluating is the number/ percentage of PCR or optical duplicates. In order to achieve that I have done the following steps:

1. aligned the fastq reads against the reference genome with `bwa mem`
2. the `.sam` output from `bwa mem` was converted to `bam`, sorted and indexed
3. the duplicates were marked with the command `MarkDuplicates` from picard

Then if I call `samtools flagstat` on the sorted `bam` file which had the duplicates marked with picard I get:

```
26595942 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
466 + 0 supplementary
1636809 + 0 duplicates
24969064 + 0 mapped (93.88% : N/A)
26595476 + 0 paired in sequencing
13297738 + 0 read1
13297738 + 0 read2
21102678 + 0 properly paired (79.35% : N/A)
24049752 + 0 with itself and mate mapped
918846 + 0 singletons (3.45% : N/A)
462212 + 0 with mate mapped to a different chr
182138 + 0 with mate mapped to a different chr (mapQ>=5)

```

Which sounds a bit strange to me, as no duplicates were found.

However, the `MarkDuplicates` command also outputted a metrics file (`dups/dupsMetrics.txt`) that did identify some duplicates. This was achieved with the following command:

```
java -jar picard.jar MarkDuplicates I=alignment/sample.sorted.bam O=alignment/sample_markDup.bam M=dups/dupsMetrics.txt

```

I see that the value under the column `PERCENT_DUPLICATION` in `dups/dupsMetrics.txt` is 0.065555 (given the column name, I would expect that the value is indeed 0.06%, and not 6.55%)

What would be the explanation between the discrepancy in the output of `samtools flagstat` and `MarkDuplicates`. Is it reasonable to expect that the value displayed by `samtools flagstat` is **rounded** to zero, given the low percentage which was observed by `MarkDuplicates` (0.06%)?

# Answer
> 5 votes
There are duplicates, in this line:

1636809 + 0 duplicates, gives 1636809/26595942 = 0.06154356

According to samtools documentation for flagstat:

> Provides counts for each of 13 categories based primarily on bit flags in the FLAG field. Each category in the output is broken down into QC pass and QC fail. In the default output format, these are presented as "#PASS + #FAIL" followed by a description of the category.

So you can usually take whatever is in the first column

# Answer
> 2 votes
@StupidWolf's answer is correct -- that first number in the flagstat output is what you want to look at to see the number of reads marked as duplicates. I wanted to add that the number given in the Picard metrics file is in fact 6%, not 0.06%. This is due to an ambiguity that is widespread in Picard metrics files; in many places the program emits a fraction even when the column name is "percent of something". Soo Hee Lee wrote a very short GATK blog post about this way back in the day, see the archived version here.

---
Tags: sam, ngs, samtools, gatk, picard
---


# Question
Title: What samples can be used for a Mutect2 Panel of Normals
I'm working on calling somatic variants from solid tissue tumors. I have one normal and one tumor sample for each patient. I plan to use Mutect2 to call somatic variants after preprocessing the data with the GATK Best Practices pipeline. Mutect2 recommends the use of a Panel of Normals (PON), but I can't figure out what samples should be used to built the PON. The guide says

> Additionally, the samples should come from subjects that were young and healthy to minimize the chance of using as normal a sample from someone who has an undiagnosed tumor. **Normals are typically derived from blood samples.**

My normal samples are from solid tissue, not blood. Can I use the normal samples from each patient to build the panel of normals? Or would this cause problems because I'm using reusing the same normals for both the PON and as the normal sample for somatic variant calling? My first thought is that most somatic variants in individual normal samples should be filtered out by pooling all the normals together to build the PON but I wasn't sure if there are other reasons to avoid using non-blood based normal samples to build the PON

# Answer
> 2 votes
As suggested by @burger in their comment, you can indeed use solid tissue samples as normals for your PON. You are clearly aware, but I want to spell this out for others, that there are a few caveats that go with the practice of using solid tissue normals rather than blood normals. Mainly, if your normals are "tissue-adjacent", ie sampled from an area of tissue adjacent to the tumor, there is a higher risk of including contaminating tumor cells compared to blood normals. As a result you might have real somatic variants that end up filtered because they are seen in the normal. Though that typically acts more at the level of the direct Tumor-Normal comparison than at the level of the PON, since the PON draws from many normal samples and the chance of seeing the same contaminating somatic mutation in many such samples should be quite low. So ultimately that should be fine.

It's a very common practice to use whatever matched normals you have available to build your PON, and there doesn't seem to be any drawback to doing that as long as you have enough samples (at least a few dozen, ideally a few hundred or more). If you don't have enough samples in your own study, then you should seek to use other samples from young and healthy people that have been sequenced using the same technology (same library prep kits/protocols, same chemistry version etc) since the goal of the PON is to weed out recurrent technical artifacts.

---
Tags: ngs, variant-calling, quality-control, gatk, somatic
---


# Question
Title: Get start and end coordinates per chromosome
I will use GATK for SNP calling (HaplotypeCaller). I need to feed the interval file in the command, otherwise I get errors (even though I want to use the whole genome, not a subset - btw it's not the point of the question, so will not elaborate).

The interval file for GATK can be for example a `bed` file, with 3 columns: `chr_name chr_start chr_end`. I don't have this file, but have the genome and the reads.

To obtain the intervals per chr/scaf, I proceeded this way:

* converted the alignment to interval file, using `bam2bed` (BEDOPS)
* extract the columns $1,2,3 (using `cut`)
* get the min & max coordinates per chr/scaff with a small R script:

```
bd <- read.table('all.bed', h=FALSE)
bd_min <- aggregate( bd$V3 ~ bd$V1 , bd, function(x) min(x))
bd_max <- aggregate( bd$V2 ~ bd$V1 , bd, function(x) max(x))
bd_coord <- merge(bd_min, bd_max, by='bd$V1')
write.table(bd_coord, file='coords.bed', sep='\t', quote=FALSE, col.names=FALSE, row.names=FALSE)

```

These should thus represent the mapped intervals, conceptually same as not feeding the intervals at all (= considering all the sequences).

The ideas for this:

* in the bed file: `$1` is be the feature, `$2` the start position, `$3` the end position
* so, I need to get the min of `$2` for each unique item in `$1`(for the start coord) and the max of `$3` for each item in `$1` (for the end coord)

The output looks plausible:

```
$ cat coords.bed | head -3
Bla_chrm1   678 43860826
Bla_chrm10  181 20381540
Bla_chrm11  343 20367560

```

My question here:

1. is this a correct way to proceed?
2. is there a standardized way to perform this?

My main concern is that I will have to re-build this "coordinate" files for every GATK run, because I think the coordinates per chr/scaf will shift, even slightly, at each GATK run with different datasets.

My other minor concern, is that this approach is a bit slow (given the 160M lines in the current bed file), so an unix-tools solution will also be accepted. I tried to compute myself the max/min using awk (e.g. `awk '$3>max[$1]{max[$1]=$3; row[$1]=$0} END{for (i in row) print row[i]}'`), but I get different results than with the R approach...

# Answer
You can try this, since you have the bam file, all of them are aligned to the same reference genome (I hope)(and apologies for the awful awk.. been a while since I used it):

```
samtools view -H some.bam | awk '{if($1~/@SQ/){print $2"\t1\t"$3}}' | sed 's/[SL]N://g' > coords.bed

```

For the human genome I got this:

```
chr1    1   195471971
chr2    1   182113224
chr3    1   160039680
chr4    1   156508116
chr5    1   151834684
chr6    1   149736546
chr7    1   145441459
chr8    1   129401213
chr9    1   124595110

```

What I did is to get the chromosome length, and make every line start at 1.

If you use R, you can do:

```
library(Rsamtools)
CL = scanBamHeader(<bam file>)[[1]]$targets
write.table(data.frame(names(CL),start=1,end=CL),
"coords.bed",row.names=FALSE,quote=F,col.names=FALSE,sep="\t")

```
> 0 votes

# Answer
Not sure if I completely understood your question, but if you want simply an interval of your entire genome you could do instead:

```
samtools faidx MyGenome.fasta

```

Obviously that means that you need access to that, not sure if you have. The resulting fai file contains a format similar to that:

```
000000F 33203223        94      60      61
000001F 28828106        33756799        60      61
000002F 27810542        63065468        60      61

```

Where the first column is your chromosome and the second one your length. I guess you could then just do a `awk '{OFS=FS="\t"}{print $1,"0",$2-1}'` to get your file.
> 2 votes

# Answer
This probably comes too late, but for the record GATK supports multiple kinds of interval file format. BED is the most common, and there's the Picard format (.interval\_list) which is BED-like with a few tweaks. Both require coordinates. **However, GATK also supports a simpler format (.list or .intervals) that allows you to just give the contig names (eg chr1) without any coordinates. Sounds like that would satisfy your need.** Incidentally, GATK also supports using VCF files as interval file if you want to re-genotype specific regions based on variants located there (ideally paired with the interval padding arguments).

See the GATK doc article for more details: https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists

(The doc has a few formatting issues but I just notified the team so they'll fix it)
> 0 votes

---
Tags: file-formats, bed, gatk
---


title: GATK 4.1.4.0 Mutect stats output

text:
# Question
Title: GATK 4.1.4.0 Mutect stats output
I have adapted my pipeline to the new filtering strategies implemented in gatk 4.1.4.0 and while looking at its output I noticed that the stat file generated by mutect2 has a negative number:

```
callable -1.538687311E9

```

I did not find the documentation about this file content but in the gatk forum it is mentioned that is should represent the number of bases above a certain coverage (10 by default). This happened for a WGS sample, while for WES I obtain a more sensible value (well, at least positive and integer, even if maybe a tad high):

```
callable        8.6884303E7

```

Has this happened to anyone else? I compared the resulting calls obtained for the WGS with a previous version of mutect2 and they are comparable-ish (98% of the new "ones" are in common with the "old" ones), but I am still worried that something fishy is happening, expecially considering that the filtering step of the pipeline uses this stats file.

# Answer
This was answered by the OP -- it was a bug that has been fixed in a subsequent release, 4.1.5.0. See the release notes here: https://github.com/broadinstitute/gatk/releases/tag/4.1.5.0

> Made callable sites a Long to avoid integer overflow (#6303)

Issue ticket: https://github.com/broadinstitute/gatk/issues/6302
> 1 votes

---
Tags: wgs, gatk
---


# Question
Title: Set new Idents based on gene expression in Seurat and mix n match identities to compare using FindAllMarkers
I am relatively new to Bioinformatics and scRNA-seq data analysis. I am using Seurat V3 to analyze a scRNA-seq dataset in R. Currently, I have merged three scRNA-seq samples from the same donor into one Seurat object, *All\_Samples*. We'll call them:

1. Uninfected
2. Virus1
3. Virus2

```
All_Samples <- merge(x = Uninfected, y = c(Virus1, Virus2))

```

Uninfected cells did not receive virus (i.e. express no viral genes), the other two samples were infected with a virus introducing new genes into the cells. I would like to accomplish two things.

**First**: I want to create a new Ident for Virus1 and Virus2 samples based on expression of a viral gene of interest: "GeneA". All cells from Virus1 and Virus2 samples that have "GeneA" expression \> 0.5 would be labeled "Pos", those with "GeneA" \< 0.5 would be labeled "Neg" in a the new Column called GeneA.

I am able to subset the objects based on GeneA expression that applies to all samples in the object. For example:

```
All_Samples_GeneA_Pos <- subset(All_Samples, subset = GeneA > 0.5)
All_Samples_GeneA_Neg <- subset(All_Samples, subset = GeneA < 0.5)

```

I assume I would have to modify All\_Samples@meta.data based on this post and this post, but I admit I'm not entirely sure how to implement the suggested answers into my data.

**Second**: Once I have separated the data from Virus1 and Virus2 into "Pos" and "Neg" cells for "GeneA", I want to look for differentially expressed genes between all Uninfected cells and "Pos" cells from Virus1 or Virus2. I have been using something like this to compare samples:

```
 FindMarkers(object = All_Samples, group.by = 'virus', ident.1 = "1", ident.2 = "2")

```

But how would I do something equivalent to this:

```
 FindMarkers(object = All_Samples, ident.1 = "Uninfected", ident.2 = "Virus1_GeneA_Pos")

```

Or maybe

```
 FindMarkers(object = All_Samples, ident.1 = "Uninfected_GeneA_Neg", ident.2 = "Virus1_GeneA_Pos")

```

Let me know if I can clarify any points. I appreciate any help!

# Answer
For the first question, you can use `ifelse()` to create a new column in the `meta.data` slot:

```
All_Samples@meta.data$new_column <-
  ifelse(
    rownames(All_Samples@meta.data) %in% colnames(All_Samples_GeneA_Pos),
    "GeneA_Pos",
    "GeneB_Pos"
  )

```

`colnames(seurat_object)` provides a vector of cell names in a given `Seurat` object. Here whatever cell that is in the `All_Samples_GeneA_Pos` object would be `GeneA_Pos` and whatever is not `GeneB_Pos`. To better control the behavior, you can use a "nested" `ifelse()`; you can have another `ifelse()` instead of the `"GeneB_Pos"` bit above.

For the second question, the rationale is the same as above, you can define an extra column with the "identities of interest", then use these identities with `SetIdent()` and then in turn using two identities of interest in the differential expression call.
> 1 votes

---
Tags: r, scrnaseq, seurat, covid-19, sars-cov-2
---


# Question
Title: Attempting to generate a heat map like in this figure
I have a data frame

```
Response TO CHEMOTHERAPY    Gene    CCF
Responders                  POLQ    1
Responders                 CDKN2A   1
Responders                  TP53    1
Responders              CCDC102B    0.906729075
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders             TSHZ3    1
Non-responders             PIK3CA   1
Responders                 SLIT2    1
Non-responders              SLIT2   1

```

I have `two groups of patients`; `Responders to chemotherapy and non-responders`. For each group I have calculated `cancer cell faction (CCF)` for a fixed set of genes so that an individual gene has several reads in each group

If `CCF = 1` the gene clonal and if CCF \< 1 gene is sub clonal, for a given gene for each group, I therefore would have the number of clonal and sub clonal for this gene, for example 2 clonal 7 sub clonal which would be 2/7 something like the below heat map

By your kindly code I produced this

Which is not well informative

This is another potential example heat map

Which show the clonal frequency for each group

By the way thanks for your kindness

I don't know how to use `dput()` but this is my complete data

https://www.dropbox.com/s/ol7m4nlkuyuy85a/c1.txt?dl=0

# Answer
> 1 votes
Given the lack of effort from the OP, probably I shouldn't have wasted my time on this question with a lot of problems:

i) The question is not reproducible, the OP should have supplied their data with `dput()`.

ii) The data as pasted above does not quite make sense to me, why would you have three readings from the same gene and even in such a case how should these replicates be represented in a heatmap?

iii) The question is not specific enough and should not go like "how can I generate a heatmap" and should be more like "in an attempt to generate a heatmap like in this figure I have tried this piece of code but ..."

Having said these, and after downvoting the question, here is a solution leveraging the `ComplexHeatmap` package:

```
library(data.table)
library(ComplexHeatmap)

my_data <-
"
Sample  Gene    CCF
Sample1 POLQ    1
Sample1 CDKN2A  1
Sample1 TP53    1
Sample1 CCDC102B    0.906729075
Sample1 CCDC102B    1
Sample1 CCDC102B    1
Sample1 CCDC102B    1
Sample2 TSHZ3   1
Sample2 PIK3CA  1
Sample2 SLIT2   1
Sample2 SLIT2   1
"

my_data <- fread(my_data)

my_data <- dcast.data.table(my_data,
                            formula = Gene ~ Sample,
                            value.var = "CCF",
                            fun=mean)

# dcast generates NA values as the provided data does not have values for each gene
my_data[is.na(my_data)] <- 0

# Heatmap() expects a matrix input
my_matrix <- as.matrix(my_data[,c("Sample1", "Sample2")])
rownames(my_matrix) <- my_data$Gene

Heatmap(my_matrix,
        cell_fun = function(j, i, x, y, width, height, fill) {
          grid.text(my_matrix[i, j], x, y, gp = gpar(fontsize = 10))})

```

This is how the data look after `dcast`, in which multiple entries per gene (my point ii above) is summarized with `mean()`:

```
> head(my_matrix)
           Sample1 Sample2
CCDC102B 0.9766823       0
CDKN2A   1.0000000       0
PIK3CA   0.0000000       1
POLQ     1.0000000       0
SLIT2    0.0000000       1
TP53     1.0000000       0

```

Here is the resulting heatmap:

---
Tags: r, ggplot2
---


# Question
Title: How to filter the PDB for single-domain proteins with full length 3D-structures solved
I am trying to reproduce a machine learning model that has been developed here. As one of the datasets, they use single-domain proteins with full length 3D-structures solved. Since I don't have a background in biology, I am not sure how to filter the PDB for these kind of proteins. Any help would be appreciated.

# Answer
If you want to reproduce the model in that paper, you should email the corresponding author asking for the exact set of PDB IDs/chain IDs/residue ranges they used. There is too much flexibility in the criteria as given to get the exact set any other way.

If you want to assemble a similar set for your own work, the RCSB PDB advanced search is very powerful and would let you filter e.g. by max sequence length or for only monomers.

However the PDB doesn't deal with domains (regions that evolve and fold independently) directly. A PDB chain may contain part of a domain, one domain, or multiple domains.

If you want a dataset of protein domains you could use CATH, SCOP or ECOD.

As a side issue, please remember to remove overlap of homologous proteins between the training and test sets. People entering the field often neglect to do that. See for example this comment piece.
> 3 votes

---
Tags: proteins
---


# Question
Title: Relationship between the instantaneous rate of change of a continuous time Markov model and the rate of a global clock model?
If continuous-time Markov models contain a parameter q that denotes the instantaneous rate of change and a global/strict clock model is the overall rate of evolution, how are these two parameters related? That is, for a given substitution model (e.g., JC69, F81, GTR), is there a deterministic relationship between the rate parameter(s) of these substitution models and the rate of the clock model? For instance, should the individual instantaneous rates of change (e.g., A -\> C, T -\> C, etc.) sum up to the global clock rate? If this is not the case, how should one understand the relationship between the rate of a continuous-time Markov model and the overall evolutionary rate?

# Answer
> 1 votes
Essentially you are comparing a "local clock" with a global clock and whether the instance rate relates to the global it depends on the model and in MCMC clocks in most cases the answer is no. A global clock is a very specific and very limiting model (below). Everything below is a comparison between taxa not between alignment positions.

"Local clock" in the context of MCMC clocks is in inverted commas because they are using population models. I will use the nomenclature here for convenience to describe rate heterogeneity between branches, not within a branch (last example).

A "local clock" has mutation rate heterogeniety across the tree whilst a global clock is completely homogenious. You can't say how is a "local clock" related to global clock because there are numerous different models.

Global clocks never hold for biological data because there are always some taxa evolving at a faster rate than others. This is particularly true in infectious disease.

One historic model used in traditional local clocks (no inverted commas) was autocorrelation, where a local clock would relate to a global clock on a "clade by clade" basis. That is clades are correlated by their mutation rate, so some clades evolve faster and other clades slower. Statastically an autocorrelated clock always outperforms a global clock under a likelihood ratio test.

I remember autocorrelation clearly because I implemented a variant of this approach, so I needed to understand how they did it, rather than click a few buttons on a GUI, which is what happens in Beast.

Beast has at least 3 population models and off-hand I cannot remember the population scenarios. I do remember the population scenarios made sense and secondly they can give 10-fold differences in rate in the data set I looked at, You need to read on the background theory for each model and assess what is more biologically appropriate for a given context. One model (from dim memories) reminded me of a classic population growth equation, i.e. an exponential. Population growth and collapse in phylodynamics would be informative of the rate, but again you need to go to the underlying theory.

The molecular clock theory that the rate of evolution towards the tips of trees is much faster than at the deeper branches is the sort of model that could be implemented in a MCMC clock. In the latter case the relationship between a global clock and the rate heterogeneity clock (which is a better description in that model than "local clock") is quite different from global vs. "local clock".

The fundamental issue with clocks is exemplified in the 10-fold difference issue. The lesson that is important to me is very easy to get a "wrong answer". The lesson important to you is if two MCMC models can't agree the same rate how can both obtain the same answer as a global clock? I am personally a fan of maximum heterogeneity, but in this day and age they will ask you to perform a Beast model to support it.

---
Tags: fasta, phylogenetics, phylogeny, molecular-genetics, molecular-clock
---


# Question
Title: Autodetect max number of cores and pass as an argument in Nextflow
I am creating a pipeline in Nextflow. One step is creating a pangenome with Roary. Roary takes threads as an argument and if no number of threads is supplied as an argument it defaults to one.

Is there a way in Nextflow to pass the maximum number of available threads as an argument? At the moment the best I can come up with is defining an arbitrarily large number of threads as an argument which will cap out and run with the max number of threads. I would like to know, is there a 'cleaner' or more canonical way of doing this?

At the moment my process looks like this:

```

process roary {

    publishDir "${params.outdir}/roary", mode: 'copy'

    cpus 200

    input:
    file gff from gff.collect()

    output:
    file("*") into roary
    file("pan_genome_reference.fa") into pan_genome
    file("pan_genome_sequences/*") into alignment_files
    file("gene_presence_absence.Rtab") into gene_presence
    set file("*accessory*"),
        file("*.Rtab"),
        file("_*"),
        file("*.txt"),
        file("*csv") into roary_out

    script:
    """
    roary -p ${task.cpus} -e -n -v -z $gff
    """
}

```

Thanks in advance. Also, apologies as I can see why this might not strictly be bioinformatics, but it's where Nextflow has the most traction.

# Answer
> 5 votes
I have not read the full Nextflow manual in depth, but I think there is no Nextflow command to detect the maximum number of threads you can use.

However, you can try to automatize this task by using the `nproc` command.

`nproc` gives you the number of processing units available. An idea would be to include something such `nproc - 1` as the max number of threads for Roary.

# Answer
> 3 votes
Here's an example I wrote recently:

```
#!/usr/bin/env nextflow

params.requestedCPU = 80
params.requestedMem = 50

maxcpus = Runtime.runtime.availableProcessors()
requestedCpus = Math.round((params.requestedCPU * maxcpus) / 100)

maxmem = Runtime.runtime.totalMemory() / 10241024
requestedMem = (Math.round((params.requestedMem * maxmem) / 100))

process foo {
    echo true
    cpus = requestedCpus
    memory = requestedMem.GB

    script:
        """
        echo $requestedCpus
 echo $requestedMem
        """
}

```

You can pass the number of CPUs or Memory directly on the command line in percentage, e.g. `nextflow run main.nf --requestedCPU 100` to have all CPU usages. Hope this helps.

---
Tags: workflow-management, nextflow
---


# Question
Title: Snakemake on cluster produces empty output files
# Problem:

When I run my Snakefile on our cluster, all the processes complete successfully. I can see from the logs that my shell scripts are being run correctly and are claiming to complete successfully. About half-way through when analysis steps begins, the workflow crashes due to empty input files (which is normal for PLINK-1.9, what's not normal is that the input files are empty). I went back and checked, most of the data-containing files (`*.map` and `*.ped`, in some cases/processes Plinks binary file set of `*.bed`, `*.bim` and `*.fam`) are completely empty since the beginning of the workflow? (Logs still claim processes all completed and created output files successfully?)

# Nuances:

* When I was playing around, I accidentally messed up the dynamic naming variables (escaped the variable interpret syntax, causing literal usage) and caused creation of single file with actual variable names, not dataset names. When I did this however, the file was not empty?
* Some 'side' files produced by the software (E.g. PLINK-1.9 produces a `*.nosex` file indicating samples with missing sex information) are not empty but the actuall data-containing output is empty?
* Had a problem earlier using `cut -f 4 {params.prefix}.unlifted | sed '/^#/d' >> {params.prefix3}` *(Intention was to collate list of entries ID's to exclude in next workflow step)*, where the use of \>\> with Snakemake on cluster environment simply would not be allowed (Permission denied error despite it being created line before by the same script?). Had to fix by adding in `chmod a+rwx {params.prefix3}` just before. Might this be a symptom of the same problem?
* This is being executed on a PBS/Torque cluster environment where the scripts are run as user who submitted them so permissions are not *supposed* to be an issue. Cluster profile was built for Snakemake using this PBS/Torque Profile. *(I just stripped out dependency flags because they were breaking things)*

# Code:

Here is the first three processes. So far, it crashes at the third `CLEAN` rule as the software PLINK-1.9 does not play nice with empty inputs.

rule TRIM:
    input:
        expand("rawData/{{sample}}.{extension}", extension=['bed', 'bim', 'fam'])

    output:
        expand("Results/{{sample}}_{{location}}_TRIMMED.{fileExtensions}", fileExtensions=["map", "ped"])

    params:
        fromBP = lambda wildcards: config["locations"][wildcards.location]["GRCh37"]["from"],
        toBP = lambda wildcards: config["locations"][wildcards.location]["GRCh37"]["to"],
        chr = lambda wildcards: config["locations"][wildcards.location]["GRCh37"]["chromosome"]

    shell:
        """
        module load plink-1.9
        plink --bfile rawData/{wildcards.sample} --chr {params.chr} --set-missing-var-ids @_# --recode --keep-allele-order --from-bp {params.fromBP} --to-bp {params.toBP} --out Results/{wildcards.sample}_{wildcards.location}_TRIMMED
        """

#/********* LiftOverPlink (LIFTOVER) *********/

rule LIFTOVER:
    input:
        expand("Results/{{sample}}_{{location}}_TRIMMED.{fileExtensions}", fileExtensions=["map", "ped"])

    output:
        outMap="Results/{sample}_{location}_LIFTED.map",
        outPed="Results/{sample}_{location}_LIFTED.ped",
        exclusion="Results/{sample}_{location}_EXCLUDE.dat"

    params:
        prefix=lambda wildcards: "Results/{}_{}_PRE_FILTER".format(wildcards.sample, wildcards.location),
        prefix2=lambda wildcards: "Results/{}_{}_POST_FILTER".format(wildcards.sample, wildcards.location),
        prefix3=lambda wildcards: "Results/{}_{}_EXCLUDE.dat".format(wildcards.sample, wildcards.location),
        chainFile="Binaries/hg19ToHg38.over.chain",
        LiftOver="Binaries/liftOverPlink.py",
        rmBadLifts="Binaries/rmBadLifts.py",
        sexes="rawData/1000g.sexes"

    run:
        print("Determining Liftover requirements now...")
        if config['samples'][wildcards.sample]['refGenome'] != "GRCh38":
            print("Liftover required. Dataset {} is mapped to {}".format(wildcards.sample, config['samples'][wildcards.sample]['refGenome'])),
            shell("module load liftover"),
            if config['samples'][wildcards.sample]['refGenome'] == "GRCh37" or config['samples'][wildcards.sample]['refGenome'] == "Hg19":
                print("Lifting from GRCh37 to GRCh38."),
                shell("python {params.LiftOver} -m Results/{wildcards.sample}_{wildcards.location}_TRIMMED.map -p Results/{wildcards.sample}_{wildcards.location}_TRIMMED.ped -o {params.prefix} -c {params.chainFile}"),
                shell("ls Results"),
                shell("python {params.rmBadLifts} --map {params.prefix}.map --out {output.outMap} --log {params.prefix2}.log"),
                shell("cut -f 2 {params.prefix2}.log > {params.prefix3}"),
                shell("chmod a+rwx {params.prefix3}"),
                shell("cut -f 4 {params.prefix}.unlifted | sed '/^#/d' >> {params.prefix3}"),
                shell("less {params.prefix3} | {output.exclusion}"),
                shell("awk 'NR==FNR{{a[$1]=$4;next}}{{ if (a[$1]!= 'NULL') {{ if (a[$1] == 'female') {{ $5=2; print }} else if (a[$1] == 'male') {{ $5=1; print }} else {{$5=0; print}}}}}}' {params.sexes} {params.prefix}.ped > Results/{wildcards.sample}_{wildcards.location}_LIFTED.ped")
                # ToDo: Add conditionals for other human reference genome builds
        else:
            print("No liftover required. Dataset {} is already mapped to GRCh38.".format(wildcards.sample)),
            shell("cp Results/{wildcards.sample}_{wildcards.location}_TRIMMED.map Results/{wildcards.sample}_{wildcards.location}_LIFTED.map"),
            shell("cp Results/{wildcards.sample}_{wildcards.location}_TRIMMED.ped Results/{wildcards.sample}_{wildcards.location}_LIFTED.ped"),
            shell("touch Resuts/{wildcards.sample}_{wildcards.location}_EXCLUDE.dat")

rule CLEAN:
    input:
        mapFile="Results/{sample}_{location}_LIFTED.map",
        pedFile="Results/{sample}_{location}_LIFTED.ped",
        exclude="Results/{sample}_{location}_EXCLUDE.dat"

    output:
        expand("Results/{{sample}}_{{location}}_CLEANED.{extension}", extension=['bed', 'bim', 'fam'])

    shell:
        """
        module load plink-1.9
        plink --map {input.mapFile} --ped {input.pedFile} --out Results/{wildcards.sample}_{wildcards.location}_CLEANED --make-bed --keep-allele-order --exclude {input.exclude}
        """
```
</code></pre>

# Answer
Ok, so I feel like a right PawPaw. This is **NOT a snakemake problem**. The issue was multi-part but it was me, python and good ol' PBS/Torque Cluster environments.

### Empty input files:

I was using a wrapper python script by the name of liftOverPlink.py to convert my files from PLINK-1.9's binary file format into the format needed by the software liftOver by UCSC. In my infinite wisdom, I did make some custom alterations. Turns out in python, there is a big difference between `"w+"` and `"a+"` for the `open()` function *(`"w+"` = write mode and create if needed, `"a+"` = append mode and create if needed)*. `"w+"` opens or creates a file but places the cursor at the beginning of the file, effectively wiping the file if it already had content inside, wheras `"a+"` puts the cursor at the end of the files content (which is beginning of file if it had to create file). I used `"w+"` which was effectively erasing my input files content...

### Why alter the script in the first place?

As it turns out, our PBS/Torque cluster was not properly serving up the software at runtime (software in question being liftOver). From what I have been able to gather, normally liftOverPlink would create intermediate files, liftOver would run and create its files where necessary, and liftOverPlink would then step in and perform format conversion back into PLINK-1.9 format. Since liftOver was not firing, the files did not exist and liftOverPlink kicked up errors since it was never designed to handle that (functionality to detect missing files would have been redundant). I attempted to ugly-fix it by using

```
filename = open(filename, "w+")
filename.close()

```

to open and close the file, effectively forcefully creating it if needed. Unfortunately that let to the above problem. Providing an absolute path to the liftOver executable fixed runtime issues and my custom code was no-longer necessary.
> 1 votes

---
Tags: snakemake, workflow-management
---


# Question
Title: How to visualise the proportion of clonal genes to sub clonal genes
I have a data frame

```
Response_TO_CHEMOTHERAPY    Gene    CCF
Responders                  POLQ    1
Responders                 CDKN2A   1
Responders                  TP53    1
Responders              CCDC102B    0.906729075
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders             TSHZ3    1
Non-responders             PIK3CA   1
Responders                 SLIT2    1
Non-responders              SLIT2   1

```

I have `two groups of patients`; `Responders to chemotherapy and non-responders`. For each group I have calculated `cancer cell faction (CCF)` for a fixed set of genes so that an individual gene has several reads in each group

If `CCF = 1` the gene clonal and if CCF \< 1 gene is sub clonal, for a given gene for each group, I therefore would have the number of clonal and sub clonal for this gene, for example 2 clonal 7 sub clonal which would be 2/7 something like the below figure

This is another potential example which reveals what is in my mind

Which show the clonal frequency for each group

I don't know how to use `dput()` but this is my complete data

https://www.dropbox.com/s/ol7m4nlkuyuy85a/c1.txt?dl=0

# Answer
> 2 votes
I've used some dummy data for simplicity. It's fairly easy to do this with ggplot, but you have to format your matrix-like data for ggplot input.

eg :

```
# make data
set.seed(1)
ngenes <- 20
nsamps <- 10
my_dat <- matrix(rnorm(ngenes*nsamps), ncol = nsamps, nrow = ngenes)
colnames(my_dat) <- paste("samp", seq_len(nsamps))
rownames(my_dat) <- paste("gene", seq_len(ngenes))

library(ggplot2)
library(reshape2) # for melt

# format data for plotting
dat_toplot <- melt(my_dat)

#plot
ggplot(dat_toplot, aes(x = Var1, y = Var2, fill = value)) + geom_tile()

```

plots

If you need specifics on plot tweaking, google's your friend. There are plenty of tutorials on how to use ggplot for heatmaps. (eg. here, or here)

# EDIT

I loaded your data into the `dat` object.

```
library(tidyverse)
dat <- as_tibble(dat)

dat_counts <- dat %>% group_by(Sample, Gene) %>% summarise(count_subclonal = sum(CCF<1), count_clonal = sum(CCF>=1))
dat_counts$ratio <- dat_counts$count_subclonal/dat_counts$count_clonal

# set infinite values (/0) to NA
dat_counts$ratio[is.infinite(dat_counts$ratio)] <- NA

ggplot(dat_counts, aes(x = Sample, y = Gene, fill = ratio)) + geom_tile() +
  geom_text(aes(label = paste(count_subclonal, "/", count_clonal)), col = 'white') # add ratio text

```

Outputs :

---
Tags: r, ggplot2, cnv, heatmap, complexheatmap
---


# Question
Title: Getting sequences of one transcript per gene out of annotation and genome
I got a genome (`data/genome/genome.fasta`) and braker-based genome annotation (`data/genome/annotation.gff3`), now I would like to get sequences of one transcript per gene (let's say the longest one).

Cufflinks has a `gffread` program to get the transcriptome out of the genome and the annotation. Something like

```
gffread data/genome/annotation.gff3 -g data/genome/genome.fasta -w data/genome/transcripts.fasta -y data/genome/proteins.faa

```

Will do the job. However, it is pulling multiple transcripts per gene. I read the help page twice, found parameters that could potentially do what I would like to (collapsing alternative transcripts into single one), but even when I specify my candidate parameters `-M -K` I get still `.t1` and `.t2` sequences for some genes.

The task is surprisingly simple, and I got already a custom script solution, but it just feels wrong. There must be a better way!

What is a reasonable/elegant way to get nucleotide sequences of genes out of a genome and annotation? Am I missing something with cufflinks?

Related, not exactly identical:

# Answer
Here is the custom python script solution:

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

from collections import defaultdict
from Bio import SeqIO

# read file into dictionary of transcrpt->gene
# jg25690   jg25690.t1
# jg25690   jg25690.t2
# jg25691   jg25691.t1
transcripts2genes = dict()
gene2transcript_mapfilename = 'data/genome/transcripts2genes.map'
with open(gene2transcript_mapfilename, 'r') as g2t_file:
    for line in g2t_file:
        gene, trainscript = line.rstrip('\n').split('\t')
        transcripts2genes[trainscript] = gene

# construct a dictionary gene->list of transcripts
# this fasta file is output of default gffread -w
genes2transcript_seuqneces = defaultdict(list)
transcripts_fasta_filename = 'data/genome/transcripts.fasta'
transcripts_fasta_file = SeqIO.parse(transcripts_fasta_filename, 'fasta')
for seq_record in transcripts_fasta_file:
    gene = transcripts2genes[seq_record.name]
    genes2transcript_seuqneces[gene].append(seq_record)

genes_fasta_filename = 'data/genome/genes.fasta'
with open(genes_fasta_filename, 'w') as out_file:
    # iterate through genes
    for gene in genes2transcript_seuqneces.keys():
        # get all the associated transcripts
        transcripts = genes2transcript_seuqneces[gene]
        # find the longest, by default the one with index 0
        lengest = 0
        for index, tr in enumerate(transcripts):
            if len(transcripts[lengest].seq) > len(tr.seq):
                lengest = index
        # this would be nicer with Bio interface (with wrapping of sequence lines)
        # but I was lazy and blast does not care, only thing I do is converting sequence to uppercase (not even sure if necessary)
        fasta_string = '>' + gene + '\n' + str(transcripts[lengest].upper().seq) + '\n'
        out_file.write(fasta_string)

```
> 1 votes

---
Tags: genome, sequence-annotation, isoform
---


# Question
Title: Concerning the current outbreak, using bioinformatics
Dear Researches community

A question concerning the current outbreak, recently we attempting to contribute to the study of the Covid19, with a tight submission time ( 4 weeks).

Inquiring help because I'm not a virologist with an insufficient biological background. However, I assume able to do in silico analysis or ( dry lab). Moreover, **new sample** sequences from local **strain** will be provided. Assembly and annotation could be done at the first stage. Then perhaps a comparative genome investigation between the old SARS and the Cuvid19.

Furthermore, The virus estimation length that causes all this is only around 29,829 base pairs of RNA. The genome contains at least 10 predicted open reading frames (ORFs): ORF1a, ORF1b, S, 3, 4a, 4b, 5, E, M and N, with sixteen predicted nonstructural proteins being encoded by ORF1a/b. I am wondering if someone can assist by considering a particular section to focus on or a study suggests. Perhaps we can contribute at least 0.001 % to the communities.

I have been reading articles around, and seems are thousand of them every day.

so I have come up with starting list:

1. Genome assemblies, annotation
2. epitome studies Here you predict the protein sequences based on the genome that will sequence, then there are a number of tools that analyze the immunogenicity of different proteins That can help you identify parts of the COVID-19 proteins that are likely to be detected by the immune system and they, of course, can be targeted by vaccine developers.
3. comparison between old saras Vs newCoVID
4. The M protein is a good target. It is extremely abundant in the nucleocapsid and is the protein responsible for giving the virus it's the shape, Almost everyone is looking at the S protein, which is justified. But there is not as much focus on the M protein.

Your response is highly appreciated if you could attach the any reference related to the list above.

# Answer
My personal approach to something like this is to first try to reproduce existing work in order to learn to use relevant tools and understand the concepts involved. Then as you get more comfortable with the subject matter, you may spot an opportunity to add value.

If you'd like to get started with the phylogenomics side of things, you can use this public workspace developed by the viral genomics group at the Broad Institute, which includes viral genome sequences from NCBI SRA as well as some workflows for doing genome assembly and quality control: https://app.terra.bio/#workspaces/pathogen-genomic-surveillance/COVID-19
> 4 votes

# Answer
Please note that Cuvid19 is COVID-19, the WHO nomenclature. SARS2 is specifically called SARS-CoV-2 and is the virologists name for COVID-19 because that is the nomenclature of the International Committee on the Taxonomy of Viruses (ICTV) for the same virus following formal taxonomic investigation. It is important not to get the names confused.

Technically you are asking to compare COVID-19 with COVID-19. Whilst that is a standard investigation this would be technically rewritten "to investigate the molecular epidemiology of COVID-19". If you are working within the virological community then the virus is termed SARS2-CoV-2. Phylogeny is a key approach, rather than "comparative genomics".

I wish you luck with your study.
> 2 votes

---
Tags: covid-19, literature-search
---


# Question
Title: Use a specific python.exe instance with PyMOL/Windows?
I want to use RDKit with a python script in PyMOL. RDKit requires a conda environment, so I think I need to be able to specify that the instance of python that I want to use would be in that environment, instead of the default for PyMOL, so that I can run it through PyMOL. The script I made is working for me through PyCharm when I test it, by using the RDkit environment "rdkit-env" as my python interpreter, with all of the additional modules I need installed there as well.

I would like to run it through PyMOL so I can add some PyMOL specific functionality as well. I've tried using shebang to specify the RDkit environment that I made, but that does not seem to be working for me with Windows.

Could someone please help me figure out what other approach might work for this? I'm new to this area, so I would appreciate any insights. Thank you!

# Answer
> 1 votes
So it is really hard to make a given python instance run a library from a different instance. It is not just one `dist-packages` library folder, but multiple that you would need to link, and all that assuming they use the same version of python. So the best bet is to install a new instance of pymol the python module.

The module can be installed in a python that uses conda with

```
conda install -c schrodinger pymol

```

This is the incentive version of PyMOL 2 and does not require a licence. You installed rdkit, so I am assuming you used conda as `conda install -c rdkit rdkit`. But just in case you went through the pain of compiling rdkit for a given python instance without conda and want the same for PyMOL. There are guides of instructions required to compile open source. Such as https://pymolwiki.org/index.php/Windows\_Install. On a barebone OS like a Pi's Raspsian, I had to do a lot more than claimed say `sudo apt-get install python3-dev libglm-dev freeglut3-dev libglew-dev libpng12-dev libfreetype6-dev build-essential libxml++3.6-dev` and link a few things manually. So just like for rdkit, conda make life a breeze.

PS. RDKit and PyMOL? Add PyRosetta and you just my favourite modules!

---
Tags: python, pymol, conda, windows
---


# Question
Title: Extract autodocked protein-ligand connections programatically
I have 2 crystal structures, one for my protein and one for my ligand (I have several protein-ligand pairs). I am using AutoDock Vina to simulate docking, which returns another file of the ligand with updated 3D coordinates. Everything is in pdbqt format, but I can use OpenBabel to convert them if needed.

I am looking for a programmatic way (mainly python, but other languages work as well) to load both the protein and the docked ligand, then extract only the acting atoms which connect the 2 structures together (and the bonds between them).

If a direct answer is too complex, I would be happy if someone could suggest a tool/framework I could use to manipulate this sort of data easily so I can implement a solution on my own.

# Answer
With pymol as a python module (`conda install -c schrodinger pymol`, not the GUI one), it is very easy.

```
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load('protein.pdbqt', 'protein') # the second is the name within pymol instance.
     pymol.cmd.load('ligand.sdf', 'ligand')
     pymol.cmd.save('neighbours.pdb','ligand expand 3') # get atoms 3 Å around

```

Done! However, if you are doing stuff in python in the first place you might want the last line as a string: `pymol.cmd.get_pdbstr('ligand expand 3')`. Changing the selection `'ligand expand 3'` to `'byres ligand expand 3'` will save the full residues.

Alternatively you can get the atom coordinates themselves etc. with it. To see what attributes `Atom` has you have you can do this:

```
for atom in pymol.cmd.get_model('ligand expand 3').atom:
    print(dir(atom))
    break

```

Which will give a large amount of attributes, such as `atom.coord`, which is a 3-element list.
> 4 votes

---
Tags: python, protein-structure, docking
---


# Question
Title: How can I find the position of every mutation where the Allele Frequency is greater than X, in regards to a reference such as Hg19?
I have a Human genome sequenced in a BAM file (along with other files with the indels, snps, cnvs).

I want to find every mutation with regards to the reference Human genome. However the majority of "variants" are noise/read errors. I am therefore required to filter the variants via a minimum threshold based on allele freq (eg. 0,5).

Using a reference (eg Hg19), how can I find the position of every mutation where the allele frequency is greater than the threshold?

PS. I am just starting, so I don't know which are the best (command line?) tools for the job nor how to build a pipeline for this.

Help would be greatly appreciated

# Answer
Hello and welcome to the world of bioinformatics,

BAM files are compressed formats of SAM and you can read more about it here.

To answer your question, I presume that your BAM is generated from a Whole Genome library preparation. Calling mutations on a Whole Genome takes a longer time than a Whole Exome since you are looking at every position instead of the exons of your genome.

What you are looking for is a special format of files called Variant Caller Format (VCF) generated by Variant Callers tools listed below in comparison studies. There isn't a tool better than another, they all have their advantages and disadvantages. I suggest that you read about them as they also describe their pipeline approaches.

1. Comparison of three variant callers for human whole genome sequencing
2. Comparing the performance of selected variant callers using synthetic data and genome segmentation
3. Performance assessment of variant calling pipelines using human whole exome sequencing and simulated data
4. Systematic comparison of germline variant calling pipelines cross multiple next-generation sequencers
5. Variant Callers for Next-Generation Sequencing Data: A Comparison Study

Once the VCF is generated, you will also have a column named `AF` which corresponds to the allele frequency also computed when using these tools. There are tools that allow you to filter VCF files to extract the needed information with a wished criteria like SnpSift, vcftools, VariantFiltration. So, I suggest reading the above articles as a start.

As for building pipelines, it depends on what language you feel at ease at. However, if I could suggest, Nextflow is a good start. You can also use Snakemake if you like Python. Or maybe simply `bash` ? In any case, there are many pipelines already made to help you through your needs.

(If this helps, don't forget to either accept the answer, vote the answer or both).
> 2 votes

---
Tags: genome, ngs, variant-calling, allele-frequency, hgp
---


# Question
Title: Phylogenetics analysis between diverse species
Hope everyone is well.

So i was appointed with a task, and i need some guidance. Attention, the one that gave me the task may not have explained this very clearly. I will write what he exactly asked of me, and what he gave me as primary data.

So he gave me 4 bed files (human genome) with chr and coordinates(start/stop), and asked me to do: "sequencing/coordinates conservation between diverse species".

Is that doable with coordinates only or do i need sequences in fasta format (phast uses fasta format)? I just dont understand how i will find "conservation between diverse species" only with coordinates. It seems weird to me.

I have found these programs: TargetOrtho2, Phast and BUSCO(?). Are they any good? Should i continue with one of them that you may have experience with?

Thanks

# Answer
I second commenter Devon Ryan's suggestion, I'm just going to give a little more detail about how to do that. You can plug the BED files into the UCSC browser following these instructions: https://genome.ucsc.edu/goldenPath/help/customTrack.html

Now you can look at the overlap between your bed tracks (black boxes I believe) and conservation (red box).

Assuming that is all you want to do, that should do it. If you want to do further data analysis or processing, you can download the conservation track as a `.wib` file from UCSC. You can convert this file to BED as described here, and use BEDTools to analyze it together with your BED of regions of interest.
> 1 votes

---
Tags: phylogenetics, phylogeny
---


# Question
Title: Why do we scale the rate matrix of a substitution model to make the average rate of substitution equal 1?
The matrix of transition probabilities for a substitution model Q over time t is found as follows: $e^{Qt}$. Since Q and t only show up as a product, we cannot distinguish time and rate of change. So we multiple Q by a constant so that the average rate of substitution is 1. Time t is then measured by the expected number of substitutions per site.

But why does forcing the average rate of substitution to equal 1 have this effect?

# Answer
> 1 votes
@Michael's comment is correct but possibly a little difficult to understand, so I'm going to rephrase it here with an example (please let me know if I'm misrepresenting):

This matrix (pretending for a moment that we care only about purine/pyrimidine transitions and thus only two states) is directly interpretable in terms of the number of transitions per unit time: $$ Q = \begin{matrix} 0 & 1.5 \\ 0.5 & 0 \\ \end{matrix} $$

Whereas this matrix is harder to think about because it is talking about the number of transitions per (I think?) $\frac{1}{0.37}$ units of time: $$ Q = \begin{matrix} 0 & 0.37 \\ 0.37 & 0 \\ \end{matrix} $$

You will have to rescale any result anyways to talk about it in terms of "per unit time", so you might as well do it in the matrix itself so that at each point we are dealing with a more or less intuitively comprehensible numbers, given that the scaling doesn't actually change anything because time is dimensionless as you point out.

Of course, if you have a time-calibrated phylogeny then time is not dimensionless, so it might make sense to use a different scaling, such that the average rate in $Q$ is instead "per year" or "per million years" or whatever works within the allowances of floating-point precision and brains.

So it's not that rescaling has any particular effect on the math, but rather that rescaling makes it easier for human brains to think about, because we like to represent time in units.

---
Tags: phylogeny, phylogenetics, molecular-genetics
---


title: Normal score transformation

text:
# Question
Title: Normal score transformation
This question is related to a normal score transformation which is performed in the paper by RUAN, Quansong, et al. Local similarity analysis reveals unique associations among marine bacterioplankton species and environmental factors. Bioinformatics, 2006, 22.20: 2532-2538.

They claim that given a variable X with observations $x= (x\_1,\ldots,x\_n)$ that are not normally distributed, normal score transformation is performed by first creating a rank vector observations $x$, $R^x = (R\_1^x,\ldots, R\_n^x)$ and then using and then transforming to $x'=(x\_1',\ldots,x\_n)$, by $x\_i\ = \phi\_{-1} (\frac{R\_i^x}{n+1} )$ where $\phi\_{-1}$ is the inverse of cummulative normal distribution.

It's intuitively clear (but maybe a bit harder to prove?) that resulting variable $X'$ has observations normally distributed, but I have a question about time series dependence. If we had that observations $x$ were correlated, say, over 4 adjacent observations, , i.e. $x\_{n+1} = f(x\_{n},x\_{n-1},x\_{n-2})$, how does this change after the normal transformation? Can we still find a function g, so that $x'\_{n+1} = g(x'\_{n},x'\_{n-1},x'\_{n-2})$ would hold? Since we're transforming based on the rank of data, I expect this to be not possible?

# Answer
As you suspected, correlations will be heavily changed by such a transformation. The reason for this is also exactly as you described, namely that using ranks rather than values for the transformation eliminates all of the information needed to actually preserve these relationships. As an example, suppose we have a single observation with rank 1 across all observations that none the less shows a decrease in absolute value across them (i.e., there's a correlation with a negative slope over time). Such a transformation will result in a correlation of 0.

One might be able to find a different transformation that preserves this relationship, or might be able to simply pool all of the time points together before computing the ranks (this would then at least partially preserve some correlations, though depending on what's needed to properly normalize the data it may cause other problems), but in general there's no guarantee that this is possible.
> 1 votes

# Answer
It appears the original observations were not a timeseries and the OP is specifically asking about application to timeseries data.

In my experience these are two separate areas of analysis because a temporal skew may result in aggregate behaviour which violates the original transformation. For example hot weather versus cold weather.

What the OP appears to be doing is normalising each point in the time series, albeit I don't know what X represents. Personally I have doubts because the aggregate correlation is across all time points and one set of aggregrates, e.g. Ecoregion 1 is compared with another e.g. Ecoregion 2. This for example would be done via linear modelling and assessement of fit performed through many approaches, e.g. the residual is just one approach.

If data has large temporal skews your approach is unlikely to work, ultimately you don't know the variation in the time series which could far more heterogeneous than the agrregrate data. Seasonality is likely a major factor in your data and I would consider this in the first instance. There are lots of ways to transform a time series, standardisation is a classic approach. Again my personal experience is these are two separate areas of statistics and directly applying one to the other is something that should be assessed with care.

---

So just to try and explain in this scenario terms, when its cold if the variation in bacterioplankton is much greater than when its hot (smaller population sizes increases variance), i.e. at different times of year then your approach can't work because normalisation is about a fixed (constant) relationship between mean and standard deviation.
> 0 votes

---
Tags: statistics, normalization
---


# Question
Title: How to download RNAseq gene expression data from GTEx
I would like to download RNAseq gene expression data from GTEx from the liver only. What is the best way to do this?

I have tried a few things. I downloaded several files from the GTEx website (https://www.gtexportal.org/home/datasets, under "RNA-Seq Data"), but the problems I run into are:

* Not sure what is the best file: "Gene read counts" or "Gene TPMs"? What is the difference?
* I would like 1 value per sample per gene. There are many duplicate genes in these 2 files. How can I collapse them to 1 value per gene?
* How can I download the data for liver tissue only, instead of ALL data for ALL tissues?

Thank you for suggestions.

# Answer
> 1 votes
* Gene read counts are the number of reads or fragments that are mapped to each gene. So it is an integer and the values are not directly comparable between samples due to differences in sequencing depth. TPM stands transcripts per million. It is the normalized gene expression level. Basically, you first normalize read counts by gene length and then normalize the library size to 1e6 in all libraries. As for which one should you use, it depends on what analysis you want to do and what tool you are going to use.
* I have not noticed that there are duplicated genes in those files, can you show an example?
* I do not think you can download expression from one tissue, but you can always subset. The whole dataset is not very large.

---
Tags: rna-seq, public-databases, gene-expression, public-dataset
---


# Question
Title: how do i change a variable with iterator under a for loop in python
This question comes up in my bioinfo class.

I have a csv file containing gene names and the cluster they belong to. I want to run GESA on each individual cluster. So I would like to extract the list of gene names corresponding to each cluster. So I wrote code like this

```
import pandas as pd
gene_list = pd.read_csv('marker_genes.csv')
for i in range(12):
    "cluster_"+str(i) = gene_list.gene[gene_list.cluster==i]

```

I know this raises syntax error, so what is the correct way of doing it. In the end, I would want `cluster_0` to `cluster_12` each being a list of gene names.

# Answer
cluster\_0 to cluster\_12 are just strings. they cannot be the names of arrays or matrices. Strings can only be used in print statements and string processing statements.

The correct way to define an array or matrix is as follows: call the arrays as: cluster\[0\] to cluster\[12\]

then, assign:

```
  for i in range(12):
        cluster[i] = gene_list.gene[gene_list.cluster==i]

```

then whenever you need to print the names of these arrays, you use: "cluster\_"+str(i)
> 1 votes

---
Tags: python, pandas
---


# Question
Title: Loading all PPI information from Uniprot into a graph in NetworkX
I am looking for guidance to load all SwissProt PPI data for Homo sapiens into a graph (from NetworkX). In particular, the part that I'm struggling with is gathering the data from SwissProt for each human protein.

---

The data that I'm looking to extract is protein-protein interaction. This is not included in the files found here, which @rtviii (answer below) suggested to use. Rather, when say you click on an entry in SwissProt, for each entry, there is a section devoted to "Interaction" as shown in the image attached below. I would like to extract that information for each human protein, and load it into a graph data structure using NetworkX

---

EDIT 2:

How I ended up doing it:

```

    query = "reviewed:yes+AND+organism:9606"
    args = {'query': query, 'format': 'xls', 'columns':'id,entry name,interactor'}

    uniprot_url = "https://www.uniprot.org/uniprot/?{}".format(urllib.parse.urlencode(args))

    with open(<DESIRED FILE PATH>, "wb") as target:
          target.write(requests.get(url).content)

```

# Answer
The query link for all interactors would be thus : `https://www.uniprot.org/uniprot/?query=organism%3A9606+AND+interactor%3A*&sort=score`

---

There is always an option of downloading the static version of the service if you find it incovnient to deal with their API for your use case from here and have it all to yourself, pipe data into your dataframes locally. I find it useful to just broswe their repositories sometimes instead of reading the specifications.

In general, you can get a variety of datashapes from uniprot that you cannot get from their static database-image by specifying the appropriate parameters in the request. You can browse and explore parameters by expanding their GUI's **search fields** as described in this little guide. Then you can request your data like so, in python for example.

```
import requests
URL = "https://www.uniprot.org/uniprot/"
PARAMS = {'organism': 9606, 'format': 'tab', ....#Additional parameters#...}
response = requests.get(url=URL, params=PARAMS)

```
> 2 votes

---
Tags: python, uniprot
---


# Question
Title: What RNA-Seq expression value would be closest to Microarray equivalent?
I know this question may seem strange.

I'm using Spearman correlation between gene expression profiles for various reasons (I won't go into details here). As a result, I often compare RNA-Seq and Microarray samples. For preliminary analysis, I usually grab what version of the data is easily accessible (RPKM, FPKM,...), but I'd like to dig a bit more.

Intuitively, I'd think a value similar to RPKM or such would make more sense than raw counts, which is why I usually convert raw counts to RPKM (**I know they're considered obsolete for statistical analysis of the data**, what I'm doing relies on the correlation)

So I ask what RNA-Seq gene expression value would be closest to a microarray equivalent ?

*i.e.* What would (theoretically) maximize correlation between gene expression of the same sample profiled with RNA-seq and microarray. RPKM ? FPKM ? TPM ? etc.

# Answer
> 2 votes
I think it is very hard to say which are the closest because they are not really comparable. But since you are using Spearman correlation, I guess RPKM, FPKM, and TPM do not change the order of gene expression levels. You might also want to normalize RNA-seq and microarray data so that they are more comparable.

# Answer
> 2 votes
I did a comparison of cDNA count data against microarray data that was published a few years ago:

> For comparisons to published data (Fig. S2; Miller et al., 2012), a generalized linear model was fitted to the relationship between log-transformed microarray and VSTPk expression levels obtained from the ImmGen Project database, and was used to transform the microarray data into values comparable to our VSTPks.

I found that the Variance-Stabilizing Transformation that was carried out by DESeq2 was close to what I wanted, but there still seemed to be a length-based bias to the reads. I corrected this by dividing by the length of the longest gene isoform in kilobases (creating something that I called *VSTPk*).

After doing this, there were range differences between the microarray and cDNA data, so I did an additional linear transformation to get the data fitting as close as possible.

---
Tags: rna-seq, transcriptome, microarray, fpkm
---


# Question
Title: Problem with classification model of genomic data: every machine learing model predicts wrongly almost always the same subset of dataset
First of all, I'd like to apologize for any spelling or grammar mistakes.

I'm having a problem using R for a classification problem. My dataset contains ~300.000 genomic data, and the features are DNA-related features (number of dinucleotides, number of trinucleotides, the CG Content, and some more). In conclusion, I have a dataset of 300.000 rows and 84 columns (columns= features). The 84th feature is basically the classification variable (there are two classes: class 1 and class 2).

I split the dataset in train set and test set and fit the train set to create a model using different methods, like Random Forest, XGBoost, Adaboost, etc. and check the accuracy using the test set.

In every model, the accuracy of the prediction on the test set is always around 86-87%. What I noticed is that, in the test set, for EVERY MODEL, almost the same number of rows (data) is predicted right (the model predict its true class), and always almost the same data of rows (data) is predicted wrong (the model predict wrong class. For example its true class is 1 but the model predicts 2).

If I take the set of data that are predicted right for every model I've created (250.000 data from the 300.000), and build a new model out of them (splitting them again in a new train and test set), I have an accuracy of 98%.

That means that the other 50.000 data (from the 300.000), add an error at every model I'm trying to create.

But I can't just ignore them. I thought that it would be a way to create two differents models, but doesn't know how to figure this out.

Then I thought that there could be a way to find which variables from the two subsets have strong differences, but even if I figure this out, I don't know what's the way for example to scale the original dataset to prevent those differences from keeping the accuracy in 87% in every model (using the original dataset of the 300.000 rows).

In conclusion, I have a dataset that pretends like "two datasets". It consists of two subsets: the first one is always predicted right from any model (machine learning algorithm) I'm trying, and the second one is always predicted wrong. How could I build a strong model taking into account the differences?

I also took a look at the box plots of the features. The most important feature for every model is called "CG" (number of dinucleotide of cytosines and guanines in DNA). In the following image you can see the boxplot of CG feature separating the original dataset in two datasets (T = for the dataset that is predicted right in every model) and (F = for the dataset that is predicted wrong in every model). 1 stands for data of class 1 and 2 stands for data of class 2. So, T1 means data from the right predicted dataset of class 1, T2 means data from the right predicted dataset of classes 2, T means data from the right predicted dataset for both class 1 and 2, etc.

As I'm a new beginner in R and in machine learning, I hope you do understand my problem. Thanks sincerely.

# Answer
You could be seeing the effect of SNPs on haplotype-driven conditions (i.e. where a combination of multiple variants predicts the class, rather than a single variant).

Do you have any method to determine the age of particular variants (e.g. via an LD plot to find inconsistencies between variants and the surrounding block structure), so that you can see if that correlates with the likelihood of error?
> 0 votes

---
Tags: r, genomics, phylogenetics, machine-learning
---


# Question
Title: Questions regarding the detection of variants (SNPs and Indels) in mammal
Greetings,

I’m a newbie in bioinformatics working currently on my first project regarding detection of variants (SNPs and Indels) in mammal assembly genome. I used the following method to detect variants for non-model organism (Camelus dromedarius), and then predict variants effects:

1- Downloading two runs of raw sequencing reads in paired end mode using fastq-dump (version 2.9.3).

2- Raw reads were quality checked, then trimmed using Trimmomatic (version 0.36).

3- Clean reads were mapped to reference genome using BWA-MEM (version 0.7.15).

4- Merging SAM files to one big file, converting SAM to BAM, sorting reads by coordinate, and marking duplicate reads using Picard (version 2.20.1).

5- Create realignment targets using RealignTargetCreator, realign around those targets using IndelRealigner, calling variants for the Base Quality Scoring Recalibration BQSR process using HaplotypeCaller, apply BQSR using two steps of BaseRecalibrator then PrintReads, calling variants using HaplotypeCaller, previous processes were within GATK (version 3.8).

6- Calling variants again using BCFtools.

7- Overlapping the two VCF files from HaplotypeCaller and BCFtools mpileup using BEDtools.

8- Filtering variants using BCFtools.

9- Select variants SNPs and Indels in two separated files using GATK.

10- Predicting effects using SnpEff.

11- SnpEff database was built using reference genome + GTF file which was converted from reference genome GFF version 3 file using gffread.

Kindly, I would like to share my pipeline with those whom has an experience in related field, to review the above steps conducted on my task, whether it's accurate or still required additional steps, or rearrangement of the pipeline sequence?

In addition, I would like to share my results with you to make an assessment whether my results are suitable to be reported to my supervisor or not?

Your response is highly appreciated,

Kind regards,

# Answer
> 1 votes
This is close to what many people do. May I add:

* between 4 and 5 you might want to add read groups using Picard. This will allow GATK to work. For this purpose (in case picard gives you issues), some time ago I made this.
* filtering with BCFtools is somewhat uncomfortable (for me at least). I would suggest you to plot the distribution of every variable in the INFO and GENOTYPE fields of the VCF file and understand their distribution and meaning. It could vary a lot and there is no one-size-fits-all threshold. For that purpose, some time ago I made this and this.

# Answer
> 0 votes
It looks like you mostly followed GATK? It looks reasonable, though a lot will depend on how good your genome and gtf file are.

---
Tags: snp, variant-calling, snpeff
---


# Question
Title: Finding discriminating positions for a sequence alignment column
I don't have any theoretical background in biology, so forgive me if my question is a bit.. well..dumb. I'm trying to use a Monte Carlo approach to find the discriminating position of a given sequence alignment column. So essentially, given a sequence alignment, I want to find the position at which we can distinguish active and inactive conformations.

So, here is my code for the Monte Carlo method:

```
def MonteCarlo(S, max_iter, threshold):
    """
    Computes the discriminating position for a sequence column

    Input: S = sequence column, max_iter = max iterations, threshold = threshold of acceptance for boltzmann distribution (between 0 and 1)
    Output: 2 classes
    """
    n = len(S)
    S_copy = S.copy()
    num_class1 = int(0.30*n) #number in class1
    print(num_class1)
    for _ in range(num_class1): #intialisation
        rand_position = random.randint(0,len(S)-1)
        if S_copy[rand_position] == 'X':
            rand_position = random.randint(0,len(S)) #reselect if used position
        S_copy[rand_position] = 'X'

    #Monte Carlo with Insertions and Deletions

    s_entropy = entropy(S)
    n_iter = 0 #counter for iterations
    class_1 = dict() #class 1 determined by X's
    class_2 = dict() #class 2 is rest of S
    for i in range(len(S)):
        if S_copy[i] == 'X':
            class_1[i] =S[i]
        else:
            class_2[i] =S[i]

    res = score(S, class_1, class_2)

    while n_iter < max_iter:
        r = random.random()

        if r < 0.25: #25% chance to do an insertion
          (c1,c2,pos) = Insertion(S, class_1, class_2)

        elif r < 0.5: #25% chance to do a deletion
          (c1,c2,pos) = Deletion(S, class_1, class_2)

        else: # 50% chance to do a swap
          (c1,c2,pos) = Swap(S, class_1, class_2)

        res2 = score(S, c1, c2)

        #Boltzmann Distribution coefficient
        delta_S = abs(res2 - res)
        lambda_, N = 0.2, 19
        size = 6000
        ind1 = random.randint(0,size - 1)
        rv = boltzmann.rvs(lambda_, N, size=size)/N
        beta = rv[ind1] #pick float between 0 and 1 with boltzmann distribution
        if  (res2 > res) or (beta*delta_S > threshold): #if better score i.e. higher entropy gain
            res = res2
            class_1 = c1
            class_2 = c2

        n_iter += 1

    return class_1, class_2, pos

```

Here I'm using the Boltzmann distribution as a probability measure. All the other functions (swap, insertion, deletion, entropy) calculate pretty much what the name says (if you need more clarification on my side, I'd be happy to provide it!)

Let me explain what my code is trying to do: Given a sequence alignment (as a list) say \['M','M','M',...,'L','L','L'\] I first set a bunch of random positions and mark them out as 'X' (the number of Xs are determined by num\_class1). Next, in the while loop I do swaps,insertions, and deletions with the goal of obtaining a better score. This meansI'm adding, removing, or swapping positions between class1 and class2. Eventually, my code returns two classes and the discriminating position.

While my code works, I don't think it's giving the right result. When I use a base case example of a list with two letters 'M' and 'L with a clear discriminating position, I get random results.

Could someone explain theoretically (or via code) how one can find the discriminating position using a Monte Carlo method as above?

Thanks!

# Answer
Without asking further questions, my suspicion would be how indels (gaps) are counted in your code against the aligment. Its the classic place for alignment mis-position errors. In other words, based on my experience your calculation probably performs fine, but linking the calculated alignment position against the actual observed alignment position is classic area of logical errors and needs careful monitoring of variables using code breakpoints within your loops.
> 1 votes

---
Tags: python, sequence-alignment, proteins
---


title: Is it possible to get such a plot

text:
# Question
Title: Is it possible to get such a plot
I wish to reproduce a very informative plot in google

Given a group of samples and a list of genes for which we have data how each genes has mutated and whether the alleles are homozygous of heterozygous. These are expressed as boolean values

Genes enriched in 2 pathways like cell cycle and p53 as shown in this figure

Map of functional alterations for a group of patients. Genes (rows) encoding components p53–DNA repair pathway; are affected by selected functional events (percent of samples altered and types of alteration are represented by colored squares) across group of samples in column. Alterations of the pathway are observed stacked green bar plots at bottom

I have the boolean matrix of mutated genes per sample would be grateful for assistence.

An example of my copy number data

```
> dataset_gistic$annotations
     type                event
G2   "Homozygous Deletion"   "ACVR2A"
G13  "Homozygous Deletion"   "HIST1H3B"
G73  "Heterozygous Deletion" "CCNE1"
G74  "Heterozygous Deletion" "TSHZ3"
G78  "Low-level Amplification"    "CTNNB1"
G110 "Low-level Amplification"    "CCNE1"
G111 "Low-level Amplification"    "TSHZ3"
G119 "High-level Amplification"   "PIK3CA"
G121 "High-level Amplification"   "NIPBL"
G147 "High-level Amplification"   "CCNE1"
G148 "High-level Amplification"   "TSHZ3"
>

> dput(dataset_gistic$annotations)
structure(c("Homozygous Loss", "Homozygous Loss", "Homozygous Loss",
"Homozygous Loss", "Homozygous Loss", "Homozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Heterozygous Loss", "Heterozygous Loss",
"Heterozygous Loss", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "Low-level Gain", "Low-level Gain", "Low-level Gain",
"Low-level Gain", "High-level Gain", "High-level Gain", "High-level Gain",
"High-level Gain", "High-level Gain", "High-level Gain", "High-level Gain",
"High-level Gain", "High-level Gain", "High-level Gain", "High-level Gain",
"High-level Gain", "ACVR2A", "HIST1H3B", "CDKN2A", "PTEN", "NAV3",
"LIN7A", "ARID1A", "ACVR2A", "SCN3A", "CTNNB1", "PBRM1", "EPHA3",
"POLQ", "PIK3CA", "SLIT2", "NIPBL", "MSH3", "APC", "HIST1H3B",
"ARID1B", "CDKN2A", "NOTCH1", "PTEN", "CCND1", "KRAS", "LRRK2",
"ARID2", "NAV3", "LIN7A", "SIN3A", "TP53", "ERBB2", "GPATCH8",
"RNF43", "GATA6", "SMAD4", "CCDC102B", "STK11", "SMARCA4", "CCNE1",
"TSHZ3", "CTNNB1", "PIK3CA", "SLIT2", "NIPBL", "HIST1H3B", "ARID1B",
"EGFR", "KCNQ3", "NOTCH1", "PTEN", "CCND1", "KRAS", "LRRK2",
"ARID2", "SIN3A", "ERBB2", "GPATCH8", "RNF43", "GATA6", "CCDC102B",
"STK11", "SMARCA4", "CCNE1", "TSHZ3", "PIK3CA", "NIPBL", "EGFR",
"KCNQ3", "CCND1", "KRAS", "ERBB2", "GPATCH8", "RNF43", "GATA6",
"CCNE1", "TSHZ3"), .Dim = c(77L, 2L), .Dimnames = list(c("G2",
"G13", "G17", "G19", "G24", "G25", "G38", "G39", "G40", "G41",
"G42", "G43", "G44", "G45", "G46", "G47", "G48", "G49", "G50",
"G51", "G54", "G55", "G56", "G57", "G58", "G59", "G60", "G61",
"G62", "G63", "G64", "G65", "G66", "G67", "G68", "G69", "G70",
"G71", "G72", "G73", "G74", "G78", "G82", "G83", "G84", "G87",
"G88", "G89", "G90", "G92", "G93", "G94", "G95", "G96", "G97",
"G100", "G102", "G103", "G104", "G105", "G107", "G108", "G109",
"G110", "G111", "G119", "G121", "G126", "G127", "G131", "G132",
"G139", "G140", "G141", "G142", "G147", "G148"), c("type", "event"
)))
>

```

And this is sample of my mutational data the number of mutations for each gene in each sample

```
    om$numericMatrix
       LP6005500-DNA_D03 LP6008336-DNA_G01 LP6008460-DNA_D01 LP6005334-DNA_H01
TP53                   2                 2                 2                 2
CDKN2A                 3                 1                 2                 0
       LP6007600 LP6008334-DNA_A03 LP6008334-DNA_A04 LP6008334-DNA_B02
TP53           2                 2                 2                 2
CDKN2A         0                 0                 0                 0
       LP6008334-DNA_C02 LP6008334-DNA_D02 LP6008336-DNA_H01 LP6008337-DNA_A07
TP53                   2                 1                 2                 2
CDKN2A                 0                 0                 0                 0
       LP6008337-DNA_H06 LP6008460-DNA_A04
TP53                   2                 2
CDKN2A                 0                 0
>

> dput(om$numericMatrix)
structure(c(2, 3, 2, 1, 2, 2, 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, 2,
0, 1, 0, 2, 0, 2, 0, 2, 0, 2, 0), .Dim = c(2L, 14L), .Dimnames = list(
    c("TP53", "CDKN2A"), c("LP6005500-DNA_D03", "LP6008336-DNA_G01",
    "LP6008460-DNA_D01", "LP6005334-DNA_H01", "LP6007600", "LP6008334-DNA_A03",
    "LP6008334-DNA_A04", "LP6008334-DNA_B02", "LP6008334-DNA_C02",
    "LP6008334-DNA_D02", "LP6008336-DNA_H01", "LP6008337-DNA_A07",
    "LP6008337-DNA_H06", "LP6008460-DNA_A04")))
>

```

# Answer
While it was easy to get your data after your edit, the data was not quite useful, the second data frame/matrix has 2 columns only!

I had to create a fake data set, that looks like this:

```
my_df <- data.frame(
  gene_class = rep(letters[1:4],3),
  some_value = rep(1:3*10,4)
)
rownames(my_df) <- paste("Gene-",rep(LETTERS[1:4],each = 3),1:4,sep="")

> my_df
        gene_class some_value
Gene-A1          a         10
Gene-A2          b         20
Gene-A3          c         30
Gene-B4          d         10
Gene-B1          a         20
Gene-B2          b         30
Gene-C3          c         10
Gene-C4          d         20
Gene-C1          a         30
Gene-D2          b         10
Gene-D3          c         20
Gene-D4          d         30

```

Using the `ComplexHeatmap` package, I got something similar to the figure you have posted.

```
# data frame for annotation
my_annot <- my_df[,"gene_class",drop=FALSE]

# the Heatmap() function requires a matrix
my_matrix <- as.matrix(my_df[,"some_value",drop=FALSE])

dev.off() # clear the plotting area in case not empty
ha = rowAnnotation(gene_class = my_annot$gene_class)

ht <-  Heatmap(my_matrix,
               left_annotation = ha, # gene_class
               row_names_side = "left", # position of row (gene) names
               heatmap_width = unit(0.3, "cm")*nrow(my_matrix), # cell width
               cluster_rows = FALSE, # no clustering
               cluster_columns = FALSE, # no clustering
               cell_fun = function(j, i, x, y, w, h, col) { # add text ("some value column") to each grid
                 grid.text(my_matrix[i, j], x, y)
               })
draw(ht, heatmap_legend_side = "bottom")

```

The `ComplexHeatmap` package offers enormous flexibility for customization. I strongly recommend devoting some time to it. The way I use it, I try to put all of my data/information in a single data frame/table, which I subsequently subset to create data frames and matrices for annotations and the actual heatmap(s) respectively.
> 6 votes

---
Tags: r, heatmap, figure-reproduction
---