# Question
Title: No variant found using GATK 4.0 HaplotypeCaller
I am doing variant calling on RNA-seq datasets from wheat which is hexaploid,the binary alignment (BAM) files were created using STAR version 2.6.0c and variant calling was done using GATK 4.0 HaplotypeCaller.The whole pipeline is as follows:

Alignment through **STAR 2.6.0c** :

```
STAR --runMode alignReads --runThreadN 20  --outSAMtype BAM  SortedByCoordinate --sjdbOverhang 100  --genomeDir wheat --readFilesIn G1_cleaned_R1.fastq G1_cleaned_R2.fastq --sjdbGTFtagExonParentTranscript wheat.gff3

```

Then **AddOrReplaceReadGroups** from **picard** was used with default settings:

```
 AddOrReplaceReadGroups \
   I=input.bam \
   O=Read_groups_added.bam \
   RGID=4 \
   RGLB=lib1 \
   RGPL=illumina \
   RGPU=unit1 \
   RGSM=20

```

Finally **GATK 4.0 haplotypeCaller** was applied:

```
gatk HaplotypeCaller --native-pair-hmm-threads 28  \
    -I Read_groups_added.bam -R wheat_ref.fa -O final.vcf

```

The command ran successfully as shown below:

```
20:33:01.018 INFO  ProgressMeter - TGACv1_scaffold_470693_5DS:301            102.2              45108670         441320.0
20:33:11.039 INFO  ProgressMeter - TGACv1_scaffold_724195_U:301            102.4              45111270         440625.5
20:33:21.072 INFO  ProgressMeter - TGACv1_scaffold_724579_U:301            102.5              45113880         439932.4
20:33:31.108 INFO  ProgressMeter - TGACv1_scaffold_248869_3B:301            102.7              45116440         439240.9
20:33:41.143 INFO  ProgressMeter - TGACv1_scaffold_248918_3B:301            102.9              45118990         438551.6
20:33:51.150 INFO  ProgressMeter - TGACv1_scaffold_306137_4AL:301            103.0              45121550         437866.7
20:34:01.176 INFO  ProgressMeter - TGACv1_scaffold_193559_2DS:301            103.2              45124100         437182.5
20:34:03.134 INFO  HaplotypeCaller - 100112154 read(s) filtered by: ((((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter) AND GoodCigarReadFilter) AND WellformedReadFilter)
100112154 read(s) filtered by: (((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter) AND GoodCigarReadFilter)
100112154 read(s) filtered by: ((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter)
100112154 read(s) filtered by: (((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter)
100112154 read(s) filtered by: ((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter)
100112154 read(s) filtered by: (((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter)
100112154 read(s) filtered by: ((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter)
100112154 read(s) filtered by: (MappingQualityReadFilter AND MappingQualityAvailableReadFilter)
21066679 read(s) filtered by: MappingQualityReadFilter
79045475 read(s) filtered by: MappingQualityAvailableReadFilter

20:34:03.134 INFO  ProgressMeter - TGACv1_scaffold_725922_U:301            103.2              45124592         437049.1
20:34:03.134 INFO  ProgressMeter - Traversal complete. Processed 45124592 total regions in 103.2 minutes.
20:34:04.749 INFO  VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0
20:34:04.749 INFO  PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.0
20:34:04.749 INFO  SmithWatermanAligner - Total compute time in java Smith-Waterman : 0.00 sec
20:34:04.749 INFO  HaplotypeCaller - Shutting down engine
[July 3, 2018 8:34:04 PM PKT]   org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 103.99 minutes.
Runtime.totalMemory()=4789370880

```

but still the generated vcf file i.e final.vcf does not contain any variants as shown below:

```
[dellemc@localhost bilal]$ tail final.vcf
##contig=<ID=TGACv1_scaffold_725916_U,length=500>
##contig=<ID=TGACv1_scaffold_725917_U,length=500>
##contig=<ID=TGACv1_scaffold_725918_U,length=500>
##contig=<ID=TGACv1_scaffold_725919_U,length=500>
##contig=<ID=TGACv1_scaffold_725920_U,length=500>
##contig=<ID=TGACv1_scaffold_725921_U,length=500>
##contig=<ID=TGACv1_scaffold_725922_U,length=500>
##contig=<ID=TGACv1_scaffold_731344_7BL,length=500>
##source=HaplotypeCaller
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  20

```

Is there any mistake in the haplotypeCaller command or any argument is missing ?

# Answer
> 3 votes
There are steps of the 'best practices' workflow missing, which I think is leading to the majority of reads being filtered by the haplotype caller, notice these lines in the output:

```
21066679 read(s) filtered by: MappingQualityReadFilter
79045475 read(s) filtered by: MappingQualityAvailableReadFilter

```

That's \>100m reads being filtered by the HC.

Specifically, you're missing `SplitNCigarReads`, which separates reads that bridge a splice junction (padded by Ns in the STAR alignment) and reassigns the mapping quality. At the moment this element of the workflow is not validated or available in GATK4 (RNA-Seq variant calling is largely abandonware in newer versions of GATK, AFAIK), but in the best practices WDL utilises GATK3.5 - see https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/master/rna-germline-variant-calling.wdl, line 440-477.

There are other elements of the recommended workflow missing too, but I think this is the critical bit. I won't harp on here about the unsuitability of RNA-Seq data for variant calling, as I'm all too aware there often isn't much option.

# Answer
> 4 votes
As well as the issues @sjcockell mentions, you have another potential problem you've got a hexaploid organism, by default the HaplotypeCaller (HC) is hardwired to a ploidy of 2, you need to set this to 6 to match your organisms ploidy using `-ploidy 6` as an argument. Since the ploidy is a key term in the underlying Bayesian stats the HC employs (see this link for an explanation) it's likely that some of the lower allelic frequencies that hexaploidy would generate will mean quite a few variants won't be called as they'd be statistically unlikely for a diploid organism.

---
Tags: vcf, variant-calling, gatk
---


# Question
Title: Evaluate clusters of individuals by using their sequence data
For a dataset of several hundred individuals, I applied a hierarchical clustering to generate clusters based on a functional trait that sets them apart.

My task is now to evaluate if these clusters can be supported by the nucleotide sequence data of the corresponding gene, i.e. is there more genetic similarity within each cluster than between the clusters.

For a first approach, I created a multiple sequence alignment for every cluster and calculated the the % identical sites value. There are more identical sites within each cluster than in an alignment of all sequences.

Do you know any good tools (preferably python) that can perform a more sophisticated evaluation of existing clusters based on sequence data? Are there other scoring methods with which I can evaluate the sequence similarity within the clusters (again preferably python)?

# Answer
You can calculate allele frequencies for each cluster you have to further verify if they belong to similar population, however if size of your dataset is rather small this may not work for you.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3813878/, here is an article about differentiating populations based on sequencing data.
> 2 votes

---
Tags: python, sequence-alignment, clustering
---


title: Using sickle for quality trimming

text:
# Question
Title: Using sickle for quality trimming
I am using sickle package from bioconda to trim my reads from the Illumina Miseq.

When I tried following code, it failed to run

```
sickle pe -t illumina -f file1   -r file2  -o output1  -p output2  -s combined_output -q 27 -l 25

```

When I replaced **illumina** with **sanger**, it works. But my reads are from the illumina. Could anyone clarify if using this option is acceptable ? Also, what is the optimum value of -q and -l

---

Sickle examines quality and length quality thresholds via a sliding window. This assesses whether to trim the 3'-end of reads, for example removing a sudden drop in quality. It can also assess whether the quality of the 5' end is sufficiently high for trimming.

# Answer
> 2 votes
It'd be rare to run into Illumina encoded quality scores any more, I doubt there are machines running that still produce those. Everything produced at the moment uses Sanger encoded (phred + 33) scores.

There is no optimum value for `-q` and `-l`, though reasonable values would typically be \<=5 for `-q` and ~5 for `-l`.

---
Tags: ngs, covid-19, trimming, virus
---


# Question
Title: rJava issue when trying to use ATACseq pipeline
I am trying to use the following `ATACseq` pipeline.

Unfortunately when I am importing it with `library(esATAC)`, it is giving me an error:

```
Error: package or namespace load failed for ‘rJava’:
 .onLoad failed in loadNamespace() for 'rJava', details:
  call: dyn.load(file, DLLpath = DLLpath, ...)
  error: unable to load shared object '/Library/Frameworks/R.framework/Versions/3.5/Resources/library/rJava/libs/rJava.so':
  dlopen(/Library/Frameworks/R.framework/Versions/3.5/Resources/library/rJava/libs/rJava.so, 6): Library not loaded: /Library/Java/JavaVirtualMachines/jdk-9.jdk/Contents/Home/lib/server/libjvm.dylib
  Referenced from: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rJava/libs/rJava.so
  Reason: image not found

```

There are threads that discuss similar problems  and :

But I checked the Java version and it is the right one: `64bit`:

```
d-172-25-60-75:~ nikitavlasenko$ java -version
java version "1.8.0_172"
Java(TM) SE Runtime Environment (build 1.8.0_172-b11)
Java HotSpot(TM) 64-Bit Server VM (build 25.172-b11, mixed mode)

```

I tried to do the following in `Rstudio`:

```
Sys.setenv(JAVA_HOME = "/usr/bin/java")
library(rJava)
library(esATAC)

```

But `library(rJava)` fails. I am using macOS High Sierra 10.13.4. I am not sure what can I try further. Any suggestions would be greatly appreciated.

# Answer
I uninstalled `Java` the way it is written here:

> 1. Click on the Finder icon located in your dock
> 2. Click on Applications tab on the sidebar
> 3. In the Search box enter JavaAppletPlugin.plugin
> 4. This will find the JavaAppletPlugin.plugin file
> 5. Right click on JavaAppletPlugin.plugin and select Move to Trash

Then installed the recommended latest version for mac on the official `Java` website and it fixed the issue.
> 0 votes

---
Tags: atac-seq
---


# Question
Title: Label protein strucuture residues in Chimera
I have generated proteins structure of capsid protein of two strains of viruses. I have modeled them showing Coulombic electrostatic potential described below where positive charge is blue and negative charge is red. I have opened them in two separate windows for comparison. I have following two issues:

1. How can I keep two molecules on the same plane when they are opened in two windows. Is there any function that lets us know which side or which angle the molecule is facing so that I can keep both molecules at same position for comparison.
2. I was able to label residue on the surface. But whenever I move the molecule some of the residue start moving and some remain static. I want to show that there is difference is electrostatic potential because of difference in the residue at some particular site. I have not been able to do that because those residue keep moving haphazardly?

It would be great if someone could help me with this.

# Answer
You cant do this in Chimera, but in pymol you can write a script to do this.

In terms of the view, you can use get\_view on the command line. This returns a section that you can copy and paste in - which will give you the same view every time. I recommend labelling the residue you want, then finding the position you want to view it from. Then using get\_view to get that view, which you can use every time (just copy it into a text file) for a reference

This does presume your different molecules are aligned, but you can do that in pymol as well to make sure that when you use the output of get\_view it is consistent.
> 3 votes

---
Tags: proteins, protein-structure, covid-19, 3d-structure, virus
---


# Question
Title: gene-level versus transcript-level analysis
Traditionally, RNA-seq data was quantified on gene level. Newer methods quantify on transcript/isoform level. For example, Kallisto only outputs transcript-level abundances. From the DESeq2 vignette:

> A newer and recommended pipeline is to use fast transcript abundance quantifiers upstream of DESeq2, and then to create gene-level count matrices for use with DESeq2 by importing the quantification data using the tximport package.

Is it just for consistency and/or simplicity that the values are converted to gene level? Would it better to proceed with transcript-level table or does that violate some DESeq2 (or similar tools) assumptions?

# Answer
> 6 votes
There are a variety of reasons people use gene-level quantitations.

1. Transcript-level differences are difficult to biologically interpret. Let's be honest, few groups are likely to put in the work required to determine what these might mean. Most genes have at least some level of characterization, so it's much more biologically tractable to think in terms of gene-level changes. Granted, "genes" are largely artificial constructs and in an ideal world one would always think in terms of transcripts, since they're what's really changing. But we don't live in an ideal world.
2. Statistical power increases with counts, so if you're splitting your counts among a handful of isoforms you could be missing something. In practice this isn't a big problem, since there are typically obvious major isoforms, but that won't always be the case.
3. "Mo transcripts mo problems", to channel the bioinformatics notorious B.I.G. It's common practice for people to use a subset of transcripts when using tools like Kallisto and Salmon, since the results usually end up being better. But which subset of transcripts should you use? Those in refseq? Gencode basic? A subset of one of those? Do we really understand all the consequences of choosing one over the others? (No, we don't) We all have our own preferences, but this isn't a solved problem.

# Answer
> 4 votes
To add to the list that Devon Ryan outlined (or perhaps to elaborate on point 2?):

Although Salmon/Kallisto/RSEM are the more accurate in their transcript quantification than the methods they superseded, transcript level quantification is still not as accurate as gene level quantification (see the tximport paper, which would also be the tool i'd recommend for getting gene level estimates from the transcript-level output of these tools).

That said, the tximport authors do claim that considering the transcript-level quantification in your gene-level analysis improves the accuracy of the gene-level analysis.

---
Tags: rna-seq, transcriptome, deseq2
---


# Question
Title: finding RNA-protein physical interaction
We all know that for physical protein-protein interaction, we need to find the distance between residues from PDB file of that interaction (finding distance between carbon alpha, carbon beta or centroid of two residues in PDB data of two proteins). I need to investigate the physical interaction between a protein and an RNA and I have PDB ID for this interaction but I don't exactly know to find the distance between which parts of residues (calpha, cbeta, centroid and etc.) with which parts of nucleotides (bases, phosphate group, or ribose backbone) from a PDB ID. Is there anyone who knows what can I do for this purpose?

# Answer
> 1 votes
If you were to treat the nucleotide like a ligand you could follow these steps on pymol, and find the measurements you desire.

Calculating distances to centroids is a little bit more complex, as you have to make a pseudo atom to measure against.

---
Tags: pdb, interactions
---


# Question
Title: Plot to show the expression of genes between tumor and normal
I have RNA-seq raw counts data for 50 samples. 20 Normal and 30 tumor. After differential analysis I got 30 DEGs. I want to make a violin plot showing the expression of each gene. I transformed counts to logCPM.

counts:

```
Genes   Tumor1  Tumor2  Normal1 Normal2
RP11-351J23.1   0   5   6   0
MIR17HG 989 896 8   0
RP11-563N12.2   0   0   0   0
LINC04393   0   2   16  0
RP11-336A10.4   0   0   0   0
DRH6-AB1    53  13  39  9
RP11-115J16.1   0   0   50  6
LINC70518   2   65  0   0

logCPM <- cpm(counts, prior.count=2, log = TRUE)

                Tumor1   Tumor2  Normal1  Normal2
RP11-351J23.1 11.84477 13.09301 15.71337 11.84477
MIR17HG      19.84847 19.79600 16.10350 11.84477
RP11-563N12.2 11.84477 11.84477 11.84477 11.84477
LINC04393     11.84477 12.47723 17.06532 11.84477
RP11-336A10.4 11.84477 11.84477 11.84477 11.84477
DRH6-AB1      15.72257 14.03896 18.32772 19.19283
RP11-115J16.1 11.84477 11.84477 18.68262 18.61229
LINC70518     12.44599 16.08368 11.84477 11.84477

```

I want to make violin plot for each gene showing expression for Tumors and Normals with this.

I first transformed logCPM. then tried like below:

```
ggplot2.violinplot(data=logCPM, xName='RP11-351J23.1',
    groupName='Tumor', position=position_dodge(0.8),
    backgroundColor="white", groupColors=c('#999999','#E69F00'),
    legendPosition="top")

```

# Answer
You'll want to use long-form data for everything:

```
library(dplyr)
library(ggplot2)
logCPM$gene = row.names(logCPM)
d = logCPM %>% gather(Sample, logCPM, -gene)
d$group = c(rep("Tumor", 2), rep("Normal", 2))
geneOfInterest = d %>% filter(gene == 'RP11-351J23.1')
ggplot(geneOfInterest, aes(x=group, y=logCPM)) + geom_violin()

```

That's a simple example, you can tweak it to meet your needs.
> 2 votes

---
Tags: r, rna-seq, rna, visualization, ggplot2
---


# Question
Title: Supplying adapter sequences for ATACseq pipeline
I have several `atac_seq.gz` files and I need to run the pipeline with them.

I have single end `ATACseq`, so the function `atacPipe` requires me to have some kind of `adapter sequence` for it to run. What is it? Where could I get one?

> Update

I found the adapter sequences, supplied them, but it is writing me that

> Error in atacPipe(case = list("mm10", fastqInput1 = paste0(base\_folder, : adapter1 should not be NULL for single end sequencing data

```
library(esATAC)
base_folder <- "/Users/nikitavlasenko/Downloads/sams_dev/"
conclusion <- atacPipe(case=list("mm10", fastqInput1 = paste0(base_folder, "Micro_atac12.fastq.gz"),
                   adapter1 =
"CAAGCAGAAGACGGCATACGAGATGTAGAGAGGTCTCGTGGGCTCGGAGATGT",
                   createReport = TRUE))

```

I supplied it and it is not null, so I am not sure what I am doing wrong.

# Answer
I was calling the function incorrectly. Here is how it worked:

```
conclusion <- atacPipe("mm10", fastqInput1 = paste0(base_folder, "Micro_atac12.fastq.gz"),
                   adapter1 = "CAAGCAGAAGACGGCATACGAGATGTAGAGAGGTCTCGTGGGCTCGGAGATGT",
                   createReport = TRUE)

```
> 0 votes

---
Tags: atac-seq
---


# Question
Title: Get Gene Expression Matrix from GEOquery
I am trying to get gene expression matrix for a list of genes I have for my list: Gene ID, Gene Symbol How can I get for each gene corresponding expression as array?? I suppose that I will have GENE\_ID GENE\_NAME SAMPLE\_1 SAMPLE\_2 ..... I am new in that, I don't know more details I am trying to study relation between diseases and genes(human) using gene expression as features. Thanks in advance.

# Answer
The answer to this really depends on the type of data you're retrieving from GEO. Microarray data sets should have a normalised matrix of expression values uploaded as part of the entry. `getGEO()` from the `GEOquery` package will return a list of ExpressionSets - in many cases this list will have a length of 1. An ExpressionSet has an `exprs()` accessor for retrieving the expression matrix. An example (GSE76510, an Illumina HT12 v4 BeadArray dataset):

```
> library(GEOquery)
> gse76510 <- getGEO('GSE76510')
> names(gse76510)
[1] "GSE76510_series_matrix.txt.gz"
> class(gse76510$GSE76510_series_matrix.txt.gz)
[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"
> exprs(gse76510$GSE76510_series_matrix.txt.gz) %>% head
             GSM2027280 GSM2027281 GSM2027282 GSM2027283 GSM2027284 GSM2027285 GSM2027286
ILMN_1343291  14.936287  14.983927  15.008199  14.830549  15.008199  14.789385  14.891532
ILMN_1343295  13.404168  13.329903  13.289727  13.222221  13.386299  13.240755  13.125056
ILMN_1651199   6.779805   7.189868   7.046952   6.962957   6.697924   6.850950   6.792647
ILMN_1651209   6.983803   6.769702   7.016746   6.954751   6.890472   6.982947   6.915671
ILMN_1651210   6.837042   6.999098   6.895225   6.944250   6.955914   6.917959   6.837368
ILMN_1651221   7.034363   7.025626   6.606457   6.977622   7.051120   6.971884   7.066445
             GSM2027287 GSM2027288 GSM2027289 GSM2027290 GSM2027291
ILMN_1343291  14.891532  14.983927  14.848452  14.913474  15.008199
ILMN_1343295  13.240755  13.010757  13.525999  13.009661  12.998406
ILMN_1651199   6.786784   6.820451   6.851160   7.013586   6.864316
ILMN_1651209   6.983701   6.995199   7.270760   6.901377   6.950876
ILMN_1651210   6.804955   7.033074   6.894117   6.770691   7.000786
ILMN_1651221   7.003035   6.995199   7.050040   6.923861   6.929933

```

You will be able to filter out a list of required genes using the probe annotation (available in the appropriate annotation.db - in this case, for example, `illuminaHumanv4.db`).

If the dataset you're interested in is RNA-Seq, you're probably out of luck. There's currently no requirement on data submitters to provide expression estimates with an RNA-Seq dataset, so the `ExpressionSet` that `GEOquery` returns will most likely have an empty `exprs` slot.

Again, an example (GSE117133):

```
> gse117133 <- getGEO('GSE117133')
> exprs(gse117133$GSE117133_series_matrix.txt.gz)
 GSM3271892 GSM3271893 GSM3271894 GSM3271895 GSM3271896 GSM3271897

```
> 3 votes

---
Tags: gene, gene-expression, geoquery
---


# Question
Title: Simplest way to work out structural variant type?
In VCF 4.2, a structural variant (SV) can be described with the `BND` keyword in `SVTYPE`. For example, the following example is an insertion (from https://samtools.github.io/hts-specs/VCFv4.2.pdf):

```
#CHROM    POS ID REF ALT                      QUAL FILTER INFO
2 321682  bnd V  T   ]13 : 123456]AGTNNNNNCAT 6    PASS   SVTYPE=BND;MATEID=bnd U
13 123456 bnd U  C   CAGTNNNNNCA[2 : 321682[  6    PASS   SVTYPE=BND;MATEID=bnd V

```

The keyword can be used to describe any SV event.

I would like to work out how to find out if an event is an insertion, deletion, duplication, inversion, etc. The algorithm/program will tell me the example above is an insertion.

I'm not able to find a tool that can do the conversion for me.

**PS:** The program I'm using is GRIDSS. This program would only generate `BND` in the SVTYPE field. It's not strictly a bug (BND can be used for any SV arrangement), but it's annoying in analysis.

# Answer
> 3 votes
GRIDSS's author had posted his own solution on Github. The code is:

```
simpleEventType <- function(gr) {
    return(ifelse(seqnames(gr) != seqnames(partner(gr)), "ITX", # inter-chromosomosal
           ifelse(gr$insLen >= abs(gr$svLen) * 0.7, "INS",
           ifelse(strand(gr) == strand(partner(gr)), "INV",
           ifelse(xor(start(gr) < start(partner(gr)), strand(gr) == "-"), "DEL", "DUP")))))
}

```

---
Tags: vcf, genome, variation, structural-variation
---


# Question
Title: Plotting Ramachandran Plot from more than one PDB file
I want to plot a Ramachandran Plot for around 5000 PDB files in python that is saved in a folder in my system. I am searching for a way to make it without accessing each file by hard-coding each file name. Please tell if there is a way around this.

# Answer
You could probably use PyRAMA for that. The only problem is that by PyRAMA plots to display, while you probably would like to save to a file.

One way to get around this is to:

1. Patch PyRAMA to allow `ramachandran` function to accept output image filename.

```
def ramachandran(file_name_list, output_image=None):
    ...
    plt.tight_layout()
    if output_image is not None:
        plt.savefig(output_file, dpi=300)
    else:
        plt.show()

```

If this works as expected I would really **encourage** you to submit a **pull-request** (or at least open an issue) on github.
2. Now you can write your main script:

```
from glob import iglob

from rama_lib.ramachandran_calc import ramachandran

for pdb_filename in iglob('/my/dir/with/pdb/', '*.pdb'):
    img_filename = pdb_filename.rsplit('.', 1)[0] + '.png'
    ramachandran((pdb_filename,), img_filename)

```

I did not test the code, but I hope you get an idea.

There are also other ways to make a ramachandran plot in python. Please check: , , .
> 2 votes

---
Tags: python, file-formats, pdb
---


# Question
Title: How to merge transcript sequence with same name in a FASTA file.
Suppose we have a fasta file like

```
>Seq1
GTTGAGAGGTGTATGGACACGAAAAACGAAACTGTATCCCGTGTTTAGCAAAGAAATCAT
>Seq1
AAAAACGAAACTGTATCCCGTGTTT
>Seq2
CGTGTTTAGCAAAGAAAT

```

I want to produce

```
>Seq1
GTTGAGAGGTGTATGGACACGAAAAACGAAACTGTATCCCGTGTTTAGCAAAGAAATCATAAAAACGAAACTGTATCCCGTGTTT
>Seq2
CGTGTTTAGCAAAGAAAT

```

# Answer
> 3 votes
Certainly this python script is not the most elegant solution but it will give you the desired result:

```
import sys

seqs = {}
with open(sys.argv[1], "r") as fh:
    curr = ""
    for line in fh:
        if line.startswith(">"):
            curr = line
            if curr not in seqs:
                seqs[curr] = ""
        else:
            seqs[curr] += line.strip()

for ident, seq in seqs.items():
    print("{}{}".format(ident, seq))

```

If you saved the script in seq\_merge.py just run it with:

```
python seq_merge.py your_file.fasta > result.fasta

```

# Answer
> 4 votes
This is a scenario where awk excels and yields a much simpler solution than the other examples here (though the logic is essentially equivalent to terdon’s Perl solution):

```
awk '
    />/  { id = $0 }
    !/>/ { seq[id] = seq[id] $0 }
    END  { for (id in seq) print id "\n" seq[id] }
' input.fasta \
> output.fasta

```

That said, for anything beyond a quick’n’dirty application I strongly recommend using a *proper* FASTA parser, not an ad-hoc solution.

# Answer
> 1 votes
assuming two lines per fasta record (or you can linearize) and using datamash:

```
$ cat jeter.fa | paste - - | sort -t $'\t' -k1,1 | datamash -g 1 collapse 2 | tr -d ',' | tr "\t" "\n"

>Seq1
AAAAACGAAACTGTATCCCGTGTTTGTTGAGAGGTGTATGGACACGAAAAACGAAACTGTATCCCGTGTTTAGCAAAGAAATCAT
>Seq2
CGTGTTTAGCAAAGAAAT

```

# Answer
> 1 votes
If you don't care about the order of the sequences and your file is small enough to fit in RAM, you can do:

```
$ perl -nle 'if(/^>/){
                $n=$_
             }
             else{
                $seq{$n}.=$_
            }
            END{
                print "$_\n$seq{$_}" for keys(%seq)
            }' file.fa
>Seq1
GTTGAGAGGTGTATGGACACGAAAAACGAAACTGTATCCCGTGTTTAGCAAAGAAATCATAAAAACGAAACTGTATCCCGTGTTT
>Seq2
CGTGTTTAGCAAAGAAAT

```

Or, if you are into the whole brevity thing<sup> </sup>:

```
perl -nle '/^>/?($n=$_):($s{$n}.=$_);}{print"$_\n$s{$_}"for keys%s' file.fa

```

---
Tags: fasta
---


# Question
Title: Comparing two files with overlapping regions and get their associated information
I have two file that have overlapping regions. One file has chromosomes, position, CNV data and samples.Other file has chromosomes, positions and ensemble ids of genes.

Here is some lines of two files.

```
contig  start   end     type  bi-allelic  S1  S2  S3
chr1    12048   177033    DUP FALSE   -1  22  24
chr1    12048   89237     DUP FALSE   12  10  -1
chr1    2985712 3355300   DEL   TRUE  2   2   2

```

Second file is like this

```
chr1    901877  911245  gene_id "ENSG00000187583.6"
chr1    22013   41670 gene_id "ENSG00000160075.10"
chr1    2985732 3355185 gene_id "ENSG00000142611.12"

```

I want output like this.

```
chr1    22013   41670 gene_id "ENSG00000160075.10"  DUP FALSE   12  10  -1
chr1    2985732 3355185 gene_id "ENSG00000142611.12" DEL   TRUE  2   2   2

```

Is there any way to do this?

# Answer
You can use `bedtools intersect` and cut the results:

```
$ cat file1.bed
chr1    12048   177033  DUP FALSE   -1  22  24
chr1    12048   89237   DUP FALSE   12  10  -1
chr1    2985712 3355300 DEL TRUE    2   2   2

$ cat file2.gtf
chr1    22013   41670   gene_id "ENSG00000160075.10"
chr1    22013   41670   gene_id "ENSG00000160075.10"
chr1    2985732 3355185 gene_id "ENSG00000142611.12"

$ bedtools intersect -a file2.gtf -b file1.bed -wb | cut -f 1-4,8-11
chr1    22013   41670   gene_id "ENSG00000160075.10"    DUP FALSE   12  10
chr1    22013   41670   gene_id "ENSG00000160075.10"    DUP FALSE   -1  22
chr1    2985732 3355185 gene_id "ENSG00000142611.12"    DEL TRUE    2   2

```
> 3 votes

---
Tags: bedtools, ensembl, cnv
---


# Question
Title: How to interpret PCA output statistically and biologically?
How can I interpret the PCA results statistically for biological data?

I have used FactoMineR and factoextra libraries for PCA

Scripts used:

library(FactoMineR)

```
res.PCA = PCA(df, scale.unit=TRUE, ncp=4, graph=F )
par(mfrow=c(1,2))
plot.PCA(res.PCA, axes=c(1, 2), choix="ind")
plot.PCA(res.PCA, axes=c(1, 2), choix="var")
dimdesc(res.PCA, axes=c(1,2))

library("factoextra")
fviz_pca_var(res.PCA, arrowsize = 1, labelsize = 3, repel = TRUE, col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

```

Here is my PCA output of differnt KO genotypes); how can I explain/interpret the PCA output statistically and if possible biologically: esp, how this **colors /contrib explain**

# Answer
> 5 votes
All that plots like this are telling you is that there are some genes that contribute more to the variability seen between your various samples than others. In an ideal world these genes will also be differentially expressed between your groups, but since we don't live in an ideal world that might not be the case. In general, the longer the line and the closer it is to the X axis the more it contributes to the overall variation (since it's contributing more the PC1, which is 53.5% of the variation in your data).

Since you've asked in the comments for a brief "results and description" section style summary of such a PCA I'll include something brief below. Note that one should not over-interpret PCA plots. The axes are (mathematically) orthogonal and don't necessarily correspond to anything biologically relevant.

> We performed a PCA on the variance-stabilized counts to check for batch effects and overall clustering of the data. As can be seen, the "3T" and "5T" groups cluster together along the first principal component, while the "0T" and "1T" samples cluster on the opposite side. This suggests that these groups are more similar to each other, as further evidenced by ...

Assuming you had used the `vst()` function in DESeq2 (I note you've tagged your post with RNAseq, so perhaps you're using DESeq2) then something like that above would be a reasonable description.

---
Tags: r, pca, proteomics
---


# Question
Title: Variants in Exome that affect miRNAs related to Schizophrenia
I have 3 exomes results (VCF files), I would like to study variants (SNP, INDEL) that can affect the miRNAs, which are related to a complex disease (Schizophrenia). What is the best way to assess the effects of variants on miRNAs and how can I link between the damaged miRNAs and Schizophrenia. Can I find a help or advice please?

# Answer
> 1 votes
I would run them through the Ensembl VEP, which will tell you which genes they hit and the effect they have on them. Then I would filter by consequence "mature\_miRNA\_variant", which will find you all the variants in your file that fall in the mature miRNAs. Or you could filter instead by the gene biotype "miRNA" which will tell you if the variants fall anywhere in the miRNAs (which will include the full stem loop structure). I would then export the list of miRNAs and try to find out about their targets using miRbase.

---
Tags: sequence-annotation
---


# Question
Title: How to scale the size of heat map and row names font size?
I have an expression data matrix (120X15; 15 samples and 120 genes), my heatmap looks blurred and raw names (gene names) looks very small and can not read. How can I improve my scripts?

Here is the example data

```
df<-structure(list(X_T0 = c(1.32606106633333, 0.222768019,  1.27419565966667, 0.293489379),
X_T1 = c(1.499956424, 0.293117217, 1.11955345266667, 0.324453864666667),
X_T2 = c(1.118528738, 0.548228048, 0.936319105, 0.846378238666667),
X_T3 = c(1.025082136, 0.841192647333333, 0.957743318666667, 1.10810667833333),
Y_T0 = c(0.424537206666667, 3.144197864, 1.01435004766667, 2.195831453),
Y_T1 = c(0.723243112666667, 0.576764958333333, 0.604810037, 0.896297411666667),
Y_T2 = c(0.335509156333333, 1.32037215366667, 1.209961935, 0.290895345333333),
Y_T3 = c(0.328275209, 1.15039119233333, 1.409758678, 0.0922298696666667),
Z_T0 = c(0.788329993666667, 1.03539976366667, 1.10523299, 0.885098106666667),
Z_T1 = c(1.125292329, 1.00032109266667, 0.877359233, 1.385935293),
Z_T2 = c(2.357924224, 0.740699933666667, 0.585095342, 1.98251017733333),
Z_T3 = c(0.678921448, 0.687992671666667, 1.005347728, 1.05829030233333)), row.names = c("AA", "PQ", "UBQ", "XY"), class = "data.frame")

```

Script being used

```
col <- colorRampPalette(c("red","yellow","darkgreen"))(30)
png("S7.png")
heatmap.2(as.matrix(df), Rowv = T, Colv = FALSE, dendrogram = "row", #scale = "row", col = col, density.info = "none", trace = "none", margins = c(7, 15) )
dev.off()

```

For sample data, this script generates perfect heatmap, but when I have more raw data(genes), this problem comes I've tried to change these margins = c(7, 15) value, doesn't work well and not sure how to proceed. Also, I want to reduce the height of the color key

I'd appreciate your help with this.

# Answer
> 5 votes
`heatmap.2` is very configurable, and has options to adjust the things you want to fix:

`cexRow`: changes the size of the row label font.

`keysize`: numeric value indicating the size of the key.

The size of the key is also affected by the layout of the plot. `heatmap.2` splits your plotting device into 4 panes (see the picture below), and you can control the size of the key partly by controlling the size and layout of these 4 panes.

The layout and scale of these panes can be controlled using the `lmat`, `lwid` and `lhei` arguments. For example, in your case, because you're not displaying the column dendrogram in your plot, you could set the height of the first row of the layout to be smaller, which would reduce the size of the colour key: `lhei=c(2, 10)` (this is just an example, you'd need to experiment to find values which worked well for your specific heatmap).

# Answer
> 3 votes
Make the image itself bigger (e.g., `png("S7.png", width=1000, height=1000)`).

Having said that, rethink the utility of having the labels there. Are you really going to go through and look at them if you have hundreds of lines? Probably not.

---
Tags: r, rna-seq, proteins, clustering
---


# Question
Title: Which is the amino acids classification when analyzing the solvent accessibility of a protein?
I came across an article which has defined 3 categories: Buried \[A, L, F, C, G, I, V, W\], Exposed \[ R, K, Q, E, N, D\] and Intermediate \[M, S, P, T, H, Y\] when analyzing the solvent accessibility of a protein. But I could not find any reference validating the claim.

Another article has come up with 4 divisions:

B=completely buried (0-4% exposed)
b=partly buried (4-25% exposed)
e=partly exposed (25-50% exposed)
E=completely exposed (50+% exposed)

Which is authentic when dealing with solvent accessibility, and which should be used? Also for the second article \[B,b,e,E\], which amino acid corresponds to which class?

# Answer
There is some confusion in the question.

Solvent accessibility of a residue, i.e. the area accessible to the solvent (water), depends on the 3D structure of the protein. It's not a property of an amino-acid (chemical compound such as alanine, leucine, ...).

There may be a correlation -- amino-acids with hydrophobic side chains are more likely to be buried, but in general amino-acids (chemical compounds) cannot be classified as buried or exposed.

I don't have access to the first article, but I tried running the prediction server described in the second linked article, and I got a result

```
GSGMQIFVKTFRDRLRNLLPQTITLEVEPSDTIENVKAKIQDKEGIPPDQQVLIFSRKRL
CCCCEEEEEECEEEEEECCCCEEEEEECCCCCHHHHHHHHHHHHCCCHHHEEEEECCEEC
EEEBeBbBebbEbEbEebEEeEbeBEBeeEeeBeeBbEeBeEeEEbEeEeBeBbbEEeeb

...

Porter (Secondary Structure):
H = Helix   (DSSP classes H, G and I)
E = Strand  (DSSP classes E and B)
C = Coil    (DSSP classes S, T and .)

PaleAle (Relative Solvent Accessibility):
B = very buried      (<=4% accessible)
b = somewhat buried  (>4% and <=25% accessible)
e = somewhat exposed (>25% and <=50% accessible)
E = very exposed     (>50% accessible)

```

In one place in the sequence you can see `QQ`, predicted as `eB`. So clearly glutamine is NOT assigned to a particular RSA class.

Dividing RSA (relative solvent accessibility, 0-100%) into categories is arbitrary, I don't think there is any widely used convention. The easiest way would be to just use numeric values.
> 2 votes

---
Tags: proteins, protein-structure, 3d-structure, computation, secondary-structure
---


title: To find p-value of compounds (pairwise)

text:
# Question
Title: To find p-value of compounds (pairwise)
I have next `SNPs` data.frame structure

```
> str(SNPs)
'data.frame':   1703 obs. of  4 variables:
 $ group: Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ rs1  : Factor w/ 3 levels "D/D","I/D","I/I": 1 1 2 3 3 2 1 1 1 1 ...
 $ rs2  : Factor w/ 3 levels "a/a","a/b","b/b": 3 3 2 3 3 2 2 3 3 2 ...
 $ rs3  : Factor w/ 3 levels "G/G","G/T","T/T": 2 1 2 1 1 3 1 1 2 1 ...
 ...other rs

> head(SNPs)
  group rs1 rs2 rs3 ...other rs
1     A D/D b/b G/T
2     A D/D b/b G/G
3     A I/D a/b G/T
4     A I/I b/b G/G
5     A I/I b/b G/G
6     A I/D a/b T/T

```

For example, I noticed that `rs5 4b/4a` and `rs6 G/G` in `group A` very often occur together (see below). In `group B` they very often occur together too. So I want to know - is it statistical regularity or not.

I can create table with all pairs in both groups

```
> SNPs$rs5_rs6 <- paste(SNPs$rs5, SNPs$rs6)
> tmp <- table(SNPs$rs5_rs6, SNPs$group)
> tmp

              A   B
  4a/4a G/G   1  20
  4b/4a G/G  31  83
  4b/4a G/T  14  51
  4b/4a T/T   1   0
  4b/4b G/G  37 106
  4b/4b G/T  35 119
  4b/4b T/T  11  31

```

So, now I need compare (find p-value) group A and group B: `4a/4a G/G` in group A vs `4a/4a G/G` in group B, `4b/4a G/G` in group A vs `4b/4a G/G` in group B, `4b/4a G/T` in group A vs `4b/4a G/T` in group B and etc. Only this pairs/

How can I do that? I need p-values for all

```
> tmp

              A   B
  4a/4a G/G   1  20 - p-value?
  4b/4a G/G  31  83 - p-value?
  4b/4a G/T  14  51 - p-value?
  4b/4a T/T   1   0 - p-value?
  4b/4b G/G  37 106 - p-value?
  4b/4b G/T  35 119 - p-value?
  4b/4b T/T  11  31 - p-value?

```

# Answer
You can get standardised residuals from `chisq.test` and convert these to p-values, but it's not recommended. The frequencies of each genotype depend on the counts of all other genotypes, so the p-values aren't independent and they aren't corrected for multiple testing, but if you must:

```
pnorm(-abs(chisq.test(tmp[,c('A','B')])$stdres)))

```
> 2 votes

---
Tags: r, gene, snp, statistics
---


title: Annotation with Prokka or RAST

text:
# Question
Title: Annotation with Prokka or RAST
I was experimenting Prokka and RAST annotation tools. So, I took a well-annotated swinepox virus genome from genebank *(NCBI Reference Sequence: NC\_003389.1)*.

I ran those sequences on Prokka and RAST Seed server at the same time. I can see that only a few (may be around 1%) of the genes were annotated. Most of them were predicted as hypothetical protein. And the results were comparable between Prokka and RAST.

I would assume that these tools look for similar sequences in NCBI and find the best-match protein. But looks like that is not the case.

I am looking for prediction tool(s) that will:

* find a well annotated swinepox virus genome in Genbank to predict most of the proteins.
* improve on the prediction of Prokka and Rast to assign functions to otherwise **hypothetical protein**

Would prediction tools such as Genemark or other genome annotation tools be beneficial.

Alternatively have I misunderstood the concept of annotation?

Assistance would be very welcome.

*I have attached the image for comparison. Left one the Swinepox genome in .gb format and right one is the same genome annotated with Prokka.*

# Answer
> 1 votes
Both of the annotation tools you used are designed for prokaryotic genome annotation. I would not expect them to work very well for viruses.

# Answer
> 1 votes
I was able to get good annotation using prokka by using following commands in prokka.

```
prokka --proteins reference.gb --outdir annotation  --prefix myprotein contigs.fa

```

All the reading frames were annotated. Actually the ones that are hypothetical in reference genomes are also annotated.

# Answer
> 1 votes
If you are looking for another tool, you could try eggnog mapper: it is an annotation tool based on orthology assignments. The EggNOG database happens to have a section for viral Orthologous Groups so you might obtain better results with it.

Note that emapper annotates coding sequences, not whole genomes. You need an extra step to predict CDS, for instance using Prodigal.

Basic usage with the optimized database for viral models:

```
prodigal -i my.genome.fna -o my.genes -a my.proteins.faa

```

then

```
python emapper.py -i my.proteins.faa --output polb_viruses  -d viruses

```

---
Tags: sequence-annotation, ngs, covid-19, virus
---


# Question
Title: How is PCR duplication rate computed in scATAC-seq?
Reading Cusanovich et al. (2015) I encountered the sentence:

> We mixed pairs of cell lines (HEK293T or HL-60 versus GM12878), performed combinatorial cellular indexing, and sequenced the resulting libraries to saturation (65% duplicate rate).

How is PCR duplication rate computed?

1. Cusanovich, D. A. et al. Multiplex single-cell profiling of chromatin accessibility by combinatorial cellular indexing. Science 348, 910–914 (2015).

# Answer
> 3 votes
The answer is in the supplementary methods page 3.

> Finally, all fragments in the same library with duplicate start and end coordinates were removed using Picard

---
Tags: single-cell, atac-seq, pcr, duplications
---


# Question
Title: How to calculate Moreau-Broto autocorrelation, Moran autocorrelation, Geary autocorrelation from protein (amino acid) sequences?
I want to code Moreau-Broto autocorrelation, Moran autocorrelation, Geary autocorrelation in Python 3.6 from the scratch for amino acid sequences without using any package. Unfortunately, I have not been able to come across a proper literature which shows how these autocorrelation indices have to be calculated. Can anyone give a clear explanation as to how to calculate it or at least give a link which has got a clear description of how to calculate these?

# Answer
Rajarshi Guha describes the formulas to calculate Moreau Broto autocorrelation, Moran coefficient and Geary coefficient and has a bibliography.

But why reinvent the wheel? It seems that the Autocorrelation module of the PyBioMed package does already what you want: it is a pure Python implementation without any external dependencies.
> 2 votes

---
Tags: python, gene, proteins, literature-search
---


# Question
Title: Why aren't clustal alignments stable under deletion of some of the sequences?
I'm new to pairwise and multiple sequence alignment in general, but I thought I understood how clustal works -- k-tuple distance is a cheap metric that's 'probably approximately' monotonic in the real global pairwise-alignment score, so we can just calculate it between all C(n, 2) sequences and do a simple UPGMA. Clustalw algorithm description here

But I'm observing odd results. Suppose I have a file with 60 sequences and run clustal on it, and observe that seqA and seqB are paired terminal child nodes. If I then delete a handful (say 10) of the sequences, but keep seqA and seqB, then in the *new* clustal alignment, seqA and seqB may not remain paired.

This is strange to me, because if Clustal was actually doing the k-tuple distance + UPGMA thing I thought it was doing, then all terminal pairings should always be stable under deletions of other sequences.

So... what am I missing?

# Answer
> 3 votes
As the graph of sequences depends on the sequences in the input, if you change the input sequences (by deleting some sequences) you change the graph.
The alignment between the sequences changes, and thus the relative distance between them changes too.

The deletion of some sequences won't cause any change on other sequences only if the full space of the C(n 2) pairwise comparisons would be done. However for a large number (and long) of sequences, this would require too much computational effort and time, making such methods as the k-tuple distance competitive.

---
Tags: sequence-alignment
---


# Question
Title: How to get all chemical compounds that interact with a given protein target?
I want to get all Pubchem Chemical compounds (CID) that interact with a given UniProt target. I have a list of UniProt targets and I need to study the relationship between these targets and chemical compounds.

Thanks

# Answer
The easiest solution to the problem of mapping chemical compounds to protein IDs is to use the ChEMBL compound ID to UniProt mapping.

Doing this with PubChem is not so straight forward in part because its downloadable data format is a hideous obfuscated and badly documented mess split over multiple compressed files which in turn contain multiple compressed files...... For a publicly funded project they should really do better.
> 0 votes

---
Tags: uniprot, interactions
---


# Question
Title: STAR Indexing Diference for Small and Large Genome File Output
I have a quick doubt on the output of the Genome Indexing, I have used the STAR program along with genome .fasta file and GFF file.

Genome size is 3GB, here is the file output

```
chrLength.txt
chrNameLength.txt
chrName.txt
chrStart.txt
genomeParameters.txt

```

I have another small Genome 60MB in size, I did the genome indexing, here is the file output

```
chrLength.txt
chrNameLength.txt
chrName.txt
chrStart.txt
exonGeTrInfo.tab
exonInfo.tab
geneInfo.tab
Genome
genomeParameters.txt
SA
SAindex
sjdbInfo.txt
sjdbList.fromGTF.out.tab
sjdbList.out.tab
transcriptInfo.tab

```

My point here is that, why I got the extra information for my small genome size, but I didn't get the same for the big size genome. I do apply the same procedure for the both.

here is the below information. Only difference I made for the large Genome size is (`--sjdbOverhang 99 --genomeChrBinNbits 15`) to reduce the memory, but the rest of things are same for small genome.

```
#!/bin/bash
NUM_THREADS=12
mkdir DB
STAR --runMode genomeGenerate --genomeDir DB \
    --runThreadN $NUM_THREADS \
    --genomeFastaFiles XL9_2.fa \
    --sjdbGTFfile XENLA_Frog.gtf+gff3 \
    --sjdbOverhang 99 \
    --genomeChrBinNbits 15

```

Could anyone give an idea, why there is different, I am new to this field, so I am wondering about the difference in this.

Thanks in advance.

Cheer San

# Answer
> 1 votes
You didn't get extra information with the small genome, the index for the larger one just wasn't made. My guess is that you ran out of RAM. For mammalian sized genomes ensure you have at least 30GB of RAM, likely more like 60GB. Also, make sure to look for error messages.

---
Tags: genome, reference-genome, star, genome-index
---


# Question
Title: Deciding which samples go in which batch
I have 370 samples to sequence, we probably will end up using only 96 samples per run (due to the barcode with the primers we'll use). This means running 4 batches. To minimize the batch effect I need to decide which samples are sequenced in which run.

I have several variables from these samples, some of them are categoric (sex, disease status, treatment stage, ...) and some of them are numeric (age, clinical score...).

To minimize the batch effect I wouldn't group all samples from the same category in one run, and I would distribute as uniformly as possible the samples according to the numeric variables. However, I could find many papers about adjusting, taking into account, but I couldn't find a paper about software to help designing the runs.

Is there any software to help designing each batch?

If there are several I would like to point out that I would define best as the one with easier user interface

# Answer
However it is specific for Illumina and for a limited number of plates: "plates that hold 2, 4, 8 IlluminaBeadChip chips and have 24, 48, 96 wells, respectively."
It doesn't take in consideration the numeric values:

```
library(OSAT)
inPath <- system.file('extdata', package='OSAT')
pheno <- read.table(file.path(inPath, 'samples.txt'),
                    header=TRUE, sep="\t", colClasses="factor")
pheno$AgeGrp <- runif(576, 0, 100)
gs <- setup.sample(pheno, optimal=c("SampleType", "Race", "AgeGrp"))
gc <- setup.container(IlluminaBeadChip96Plate, 6, batch='plates')
gSetup <- create.optimized.setup(sample=gs, container=gc, nSim=1000)

```

It handles both numeric and categorical data, but it can only optimize for two nominal criteria.

It can handle both numeric and categorical data for an arbitrary number of variables. It is not restricted to an even number of samples per batch.
> 3 votes

# Answer
I've had to go through this pain, and found the best way to do it was to simulate the model fit to start assessing the balance of covariates. I've got two notebooks available on Github, which should help you deal with this issue, although I've not thought about how to accurately assess continuous variables.

While it's not wrapped in a package, you can use these notebooks as a template to fit to your problem. I may wrap this in a package at some point, but I haven't had the chance to do so yet.

This Notebook will allow you to set up all your variables, max size of a batch (may vary due to multiplexing, library set up, plate size, etc). Bin them into appropriate batches, then create a design matrix, model fit, and test.

Apologies in advance that it's a bit of a mess, but if anything isn't clear, open an issue on the repo.
> 2 votes

---
Tags: ngs, software-recommendation, design
---


title: Batch Effects in RNA Seq Sample

text:
# Question
Title: Batch Effects in RNA Seq Sample
I am using the R (using EdgeR) for the RNA Seq analysis, I had few batch effect samples like Control vs treatment. Could anyone tell me the best way to remove the batch effects.

I have looked into the cluster dendrogram and MDS post (using EdgeR), the Control and Treatments are not grouping together, and their is differences in the different batch within Control and Treatment.

I tried this syntax to remove the batch effects for the rest of the analysis to get the DE Genelist. I feel it didn't able to remove the batch effects from the rest of the analysis.

```
# define the experiment
Batch <- factor(c(1,2,3,1,2,3))
Treatment <- factor(c(0,0,0,1,1,1))
design_noI <- model.matrix(~ 0 + Batch + Treatment)

# define our DGEList
San <- DGEList(counts=x,genes=gene_info)

A <- aveLogCPM(San)
y2 <- San[A>1,]
y2 <- calcNormFactors(y2)
logCPM <- cpm(y2, log=TRUE, prior.count=5)
plotMDS(logCPM)
logCPMR <- removeBatchEffect(y2, Clutch = Clutch)
plotMDS(logCPMR)

```

I can see the difference in MDS plot on removing the batch effects but I couldn't able to carry the same data flow into my rest of analysis.

```
# filter the lowly expresed genes and normalize the data
keep <- rowSums(cpm(San/logCPMR)>1) >= 3
keep2 <- removeBatchEffect(keep, Batch)
y_group <- San/logCPMR[keep2, , keep2.lib.size=FALSE]
y_group <- calcNormFactors(San/logCPMR)

```

I have mention the `San/logCPMR`, as `San` take as data without remove the batch effects, whereas I think `logCPMR` as batch effected removed data frame.

Any suggestion on this more useful to me.

Please let me know, any other information needed to deal with this batch effects.

# Answer
The help page for `removeBatchEffect()` explicitly states that it should not be used in conjunction with differential expression testing, so don't do that. You already have the batch effect included in your model, so you're already correcting for it when testing for differential expression. If your MDS plots suggest that samples are variably affected by the batch effect then look into the `combat()` function from the `sva` package.
> 3 votes

---
Tags: r, rna-seq, edger, batch-effects
---


title: Genome scaffolding

text:
# Question
Title: Genome scaffolding
I have assembled a virus genome using Ray resulting in approximately 5000 contigs. Now I want to build a scaffold of those contigs to get the full genome (I am aware of the fact that there might be some gaps in the genome).

I found one tool called Medusa to do this. Unfortunately, this tool cannot take a file above 50 MB and my contig file is about 156 MB.

Are there any other tools that can perform reference based scaffolding?

It would be great if anyone knows about genome scaffolding packages in Bioconda.

# Answer
The number of contigs and total assembly size you have suggest that there was probably more in the sequencing run than a single virus strain. Does the total assembly length correspond to expected genome size of the reference strain? Have you tried just to blast some your contigs? Could you have a contamination or something?

Maybe you can just sort out the contigs by the reference using something like artemis or quast compares your assembly to reference and showing nice stats (just to be sure that your assembly makes sense so far). Or you could use directly MUMmer (which the alignment tool that is running on the backend of of Medusa and Quast).

This is not directly answering your question "How to scaffold", however I think you should first figure out if your contigs make sense before you scaffold them. Also sorting your contigs by artemis using reference is practically the same thing as scaffolding since you won't know the sequence in between in either of the cases.
> 2 votes

# Answer
What you mean is assembly based scaffolding, as opposed to using reads with long distance information such as mate pair/long jumping distance to scaffold.

You could reduce your contig file to just the longest contigs, since they have the most information. Probably 156MB of data is overkill for your virus.

There are other programs out there.

https://github.com/ksahlin/BESST

https://github.com/institut-de-genomique/MaGuS

Your mileage will vary.
> 1 votes

# Answer
I found that Viral NGS Pipeline from broad institute can actually do this. There is a python script called **order\_and\_orient** which serves this purpose. Here is the link viral ngs assembly
> 1 votes

---
Tags: ngs, covid-19, sars-cov-2, reference-genome, scaffold
---


# Question
Title: Are mitochondrial genes to exclude in scRNA-seq such as ribosomal genes?
In this answer, it is stated that ribosomal genes should be excluded prior to normalization in scRNA-seq as contaminants.

Do mitochondrial genes have to be excluded as well? I plotted the top 50 expressed genes for a specific dataset and they tend to appear often (for example MT-ATP6). My assumption is that, given that they work for mitochondrial function and may be highly expressed, they can dilute the signal of gene differential across cell types but expressed at lower levels. Is this biologically sound?

Additionally, in this course, mitochondrial genes are used to filter cells when they contribute above a certain ratio of total RNA of single cells. However, I could not find an explanation for this filtering. What is the underlying process this procedure is used to account for?

# Answer
> 6 votes
According to Ilicic et al. (2016), on upregulation of mtRNA in broken cells:

> There is an extensive literature on the relationship between mtDNA, mitochondrially localized proteins, and cell death \[34, 35\]. However, upregulation of RNA levels of mtDNA in broken cells suggests losses in cytoplasmic content. In a situation where cell membrane is broken, cytoplasmic RNA will be lost, but RNAs enclosed in the mitochondria will be retained \[...\]

According to this, I think that mtRNA levels relative to endogenous RNA can be used as a control for low quality (broken) cells, in a similar way to ERCC spike-ins.

# Answer
> 2 votes
Mitochondrial genes do not always have to be removed, though you **may** need to do so for a variety of technical reasons.

Mitochondrial reads are innately different from rRNA due to rRNA typically being excluded during normal library preparation (so any of it getting carried over is essentially a contaminant). Some scRNAseq programs (e.g., RaceID) can use sampling at various stages, which can break if you have a lot of your signal going to mitochondrial transcripts. Then it's not that they should be excluded because they're from mitochondria, but rather that the results are heavily affected by the presence of VERY highly expressed genes and your results are suspect if you don't remove these RNAs. When your methods aren't affected by that (e.g., you use more robust scaling methods) or you don't have super high expression of mtRNA then that's not explicitly needed.

Further, sometimes you very much **want** to keep mtRNAs, since you might be interested in metabolic changes. Since I work in an institute partially dedicated to immunology, we often look at things like metabolic changes involved in immune activation and removing mtRNAs would hinder that. But if you're *a priori* not interested in metabolism then excluding genes you won't care about isn't unreasonable.

---
Tags: rna-seq, scrnaseq, quality-control
---


# Question
Title: Using a Bash Script to search TaxIDs against NCBI's Taxonomy yields "400 Bad Request" error?
I've been searching TaxIDs against NCBI's Taxonomy DB to get taxonomic lineages for species. I have successfully done this for 1,000's of TaxIDs that were returned to me by Blast+ blasts in a CSV. (The first column in the CSV was the TaxID for the top hit of a sequence, and the other columns were for GI, Scientific Name, e-value, etc.). I simply wrote a script that would take the first field, which is always the TaxID, and search it against Taxonomy to get lineages for the species those TaxIDs represented. Everything worked fine!

However, some of the records returned from my original blast (that gave me the TaxIDs in the first place) would return two or more TaxIDs per record. (I suspect because a sequence matched equally well with two or more different species?) That's fine, but I had to remove these few exceptions before I searched the rest of the regular results against Taxonomy. So I've got most of my data/lineages, but for about 1:10,0000 records that gave multiple TaxIDs, I would like to know the lineages of the multiple species that the blast suggested.

I now took all of those "more-than-one-taxid" hits, separated them out into columns in an excel sheet and saved it as a CSV, and am trying to search these against Taxonomy with essentially the same script I used for the rest of the data, but am getting an error for *every* taxID I search (see below).

**Complete Error** (for each taxid searched against Taxonomy):

```
400 Bad Request
&rettype=full&retmode=xml&edirect_os=linux&edirect=7.90&tool=edirect&email=x@xxx.localdomain'my&id=464126
Result of do_post http request is
$VAR1 = bless( {
                 '_protocol' => 'HTTP/1.1',
                 '_content' => '<?xml version="1.0" encoding="UTF-8" ?>
<!DOCTYPE eEfetchResult PUBLIC "-//NLM//DTD efetch 20131226//EN" "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20131226/efetch.dtd">
<eFetchResult>
 <ERROR>WebEnv parameter is required</ERROR>
</eFetchResult>
',
                 '_rc' => 400,
                 '_headers' => bless( {
                                        'connection' => 'close',
                                        'set-cookie' => 'ncbi_sid=9FF573148F458940_CEC6SID; domain=.nih.gov; path=/; expires=Thu, 14 Mar 2019 22:36:03 GMT',
                                        'cache-control' => 'private',
                                        'access-control-allow-origin' => '*',
                                        'date' => 'Wed, 14 Mar 2018 22:36:03 GMT',
                                        'vary' => 'Accept-Encoding',
                                        'client-ssl-cert-issuer' => '/C=US/O=DigiCert Inc/OU=www.digicert.com/CN=DigiCert SHA2 High Assurance Server CA',
                                        'client-ssl-cipher' => 'ECDHE-RSA-AES256-GCM-SHA384',
                                        'client-peer' => '130.14.29.110:443',
                                        'strict-transport-security' => 'max-age=31536000; includeSubDomains; preload',
                                        'client-date' => 'Wed, 14 Mar 2018 22:36:03 GMT',
                                        'content-type' => 'text/xml; charset=UTF-8',
                                        'x-xss-protection' => '1; mode=block',
                                        'client-transfer-encoding' => [
                                                                        'chunked'
                                                                      ],
                                        'server' => 'Finatra',
                                        'client-ssl-socket-class' => 'IO::Socket::SSL',
                                        'client-response-num' => 1,
                                        'ncbi-phid' => 'DAF6EB0134DF5AEC039F11D248910BE7.1',
                                        'ncbi-sid' => '9FF573148F458940_CEC6SID',
                                        '::std_case' => {
                                                          'set-cookie' => 'Set-Cookie',
                                                          'access-control-allow-origin' => 'Access-Control-Allow-Origin',
                                                          'client-ssl-cert-issuer' => 'Client-SSL-Cert-Issuer',
                                                          'client-ssl-cipher' => 'Client-SSL-Cipher',
                                                          'client-peer' => 'Client-Peer',
                                                          'strict-transport-security' => 'Strict-Transport-Security',
                                                          'client-date' => 'Client-Date',
                                                          'x-xss-protection' => 'X-XSS-Protection',
                                                          'client-transfer-encoding' => 'Client-Transfer-Encoding',
                                                          'client-ssl-socket-class' => 'Client-SSL-Socket-Class',
                                                          'client-response-num' => 'Client-Response-Num',
                                                          'ncbi-phid' => 'NCBI-PHID',
                                                          'ncbi-sid' => 'NCBI-SID',
                                                          'x-ua-compatible' => 'X-UA-Compatible',
                                                          'client-ssl-cert-subject' => 'Client-SSL-Cert-Subject',
                                                          'content-security-policy' => 'Content-Security-Policy'
                                                        },
                                        'x-ua-compatible' => 'IE=Edge',
                                        'client-ssl-cert-subject' => '/C=US/ST=Maryland/L=Bethesda/O=National Library of Medicine/CN=*.ncbi.nlm.nih.gov',
                                        'content-security-policy' => 'upgrade-insecure-requests'
                                      }, 'HTTP::Headers' ),
                 '_msg' => 'Bad Request',
                 '_request' => bless( {
&rettype=full&retmode=xml&edirect_os=linux&edirect=7.90&tool=edirect&email=x@xxx.localdomain',
                                        '_uri' => bless( do{\(my $o = 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi')}, 'URI::https' ),
                                        '_headers' => bless( {
                                                               'user-agent' => 'libwww-perl/6.05',
                                                               'content-type' => 'application/x-www-form-urlencoded',
                                                               '::std_case' => {
                                                                                'if-ssl-cert-subject' => 'If-SSL-Cert-Subject'
                                                                               }
                                                             }, 'HTTP::Headers' ),
                                        '_method' => 'POST',
                                        '_uri_canonical' => $VAR1->{'_request'}{'_uri'}
                                      }, 'HTTP::Request' )
               }, 'HTTP::Response' );

```

**What I've done:**

Here is a sample of the data I used that worked correctly and yielded lineages when searching the taxIDs (taxIDs in 1st column):

```
549298,Contig1,Francisella halioticida,5.44e-126,gi|1213460577|gb|CP022132.1|,CP022132
240395,Contig2,Glandirana emeljanovi,4.74e-131,gi|1025812536|gb|KU641020.1|,KU641020
1069680,Contig3,Pneumocystis murina B123,2.34e-85,gi|1124120202|ref|XM_019757798.1|,XM_019757798
3871,Contig4,Lupinus angustifolius,1.89e-93,gi|1239725835|gb|CP023132.1|,CP023132

```

Here is a sample of the **data that I am trying to get to work** (when I search these manually against Taxonomy from their web page they seem to work and give lineages; it's just that I have just enough of them to make a non-programmatic solution really time consuming and tedious):

```
258234
1316279
1316315
710996
711000
225338
985895
710996
1382731
77133
344338

```

This is the business end of **the script that I used that worked for the working data**:

```
#!/bin/bash
while IFS=',' read -r f1 f2 f3 f4 f5 f6
do
        #Firstly echo the Contig# (Query Seq. ID) followed by lineage information.
        echo -n "$f2," >> $2

        #Now get the lineage, pipe to xtract and parse XML for readability in output file
        efetch -db taxonomy \
        -id $f1 \
        -format xml \
        | xtract -pattern Taxon \
                -tab ',' -sep ',' \
                -def "NONE" \
                -element TaxId,ScientificName \
                -division LineageEx \
                -group Taxon \
                -if Rank -equals superkingdom \
                -or Rank -equals kingdom \
                ...
                -or Rank -equals species \
                -tab ',' -sep ',' \
                -def "NOT_AVAILABLE" \
                -element Rank,ScientificName \
                >> $2
done < "$1"

```

This is the slightly altered script that I am trying to use to make the not-working set of data yield results (and that produces the error):

```
#!/bin/bash
while IFS='' read -r line || [[ -n "$line" ]]; do
        efetch -db taxonomy \
        -id $line \
        -format xml \
        | xtract -pattern Taxon \
                -tab ',' -sep ',' \
                -def "NONE" \
                -element TaxId,ScientificName \
                -division LineageEx \
                -group Taxon \
                -if Rank -equals superkingdom \
                -or Rank -equals kingdom \
                ...
                -or Rank -equals species \
                -or Rank -equals subspecies \
                -tab ',' -sep ',' \
                -def "NOT_AVAILABLE" \
                -element Rank,ScientificName \
                >> $2
done < "$1"

```

As you can see, **the only difference** is that I am reading in entire lines for the TaxIDs list because it is simply a list with line breaks between the IDs as opposed to a true CSV with other values as in the example that works. I've saved the taxIDs list as both a CSV and TXT and it doesn't make a difference. Also: I can echo what $line is in the second script and it is exactly I'd expect--just the numbers of the taxID, **so why would I be getting a 400 error for bad syntax on my query if it looks exactly like what I am sending in the first, working script?**

Thanks!

# Answer
If I understand correctly, your question boils down to: parsing a bunch of numbers in a row from a .csv file gives a different result from parsing a bunch of numbers as *rows* of a text file. I that right? I suggest that you look at the source of your text file. For example, was the text file created on a Windows host (running bash, or in the command terminal?)?

Try using this modified version of the unix cat command: `cat -et text_file.txt`

The very end of each line will have a '\\$' sign and that will show you if there are any extra non-visible characters in your rows. For example, a DOS/Windows formatted file will usually have a ^M at the end of each row before the '$'. If you don't strip off those line endings then they get added to your queries, even if you cannot see them. You can strip off those linefeed (or are they carriage returns?) using a Perl one-liner like: `perl -pe 's/\r//' badfile.txt > goodfile.txt`
> 3 votes

---
Tags: taxonomy, shell
---


# Question
Title: How to efficiently extract coordinates from a `pybedtools.Interval`?
I'm using pybedtools to process data, and at some point, I loop over a `BedTool` object and use *unpacking* to extract coordinates from the `Interval`:

```
for interval in my_bedtool:
    (chrom, start, end, _, _, strand) = interval

```

Suppose I want to compute the interval length, as follows:

```
    length = end - start

```

This results in a `TypeError: unsupported operand type(s) for -: 'str' and 'str'`.

It happens that the unpacked `end` and `start` variables are strings. Converting them to ints is probably suboptimal given that `Interval`s can be created by providing ints for the coordinates:

```
from pybedtool import Interval
interval = Interval("X", 0, 100, strand="+")

```

The `start` and `end` attributes of an `Interval` are ints.

I could do:

```
    (chrom, start, end, strand) = (
        interval.chrom, interval.start, interval.end, interval.strand)

```

But I think this looks a bit ugly.

Is there a better approach?

# Answer
The reason why unpacking an `Interval` results in strings instead of ints might be related to how indexing by position works.

Even without considering the string/int issue, it seems slightly more efficient to access the interval's attributes rather than unpacking them directly:

```
In [311]: %timeit [(chrom, start, end, strand) for (chrom, start, end, _, _, strand) in my_bedtool]
306 ms ± 26.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [312]: %timeit [(interval.chrom, interval.start, interval.end, interval.strand) for interval in my_bedtool]
238 ms ± 24.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

```

To resolve the ugliness issue, attribute access can be made more compact through the use of `operator.attrgetter`:

```
In [313]: from operator import attrgetter

In [314]: get_interval_attrs = attrgetter("chrom", "start", "stop", "strand")

In [315]: %timeit [get_interval_attrs(interval) for interval in my_bedtool]
239 ms ± 26.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

```
> 2 votes

---
Tags: python, pybedtools
---


title: TRAL does not find "Phobos result file"

text:
# Question
Title: TRAL does not find "Phobos result file"
I want to use TRAL to annotate tandem repeats in the reference genome of *Caenorhabditis elegans*. For this, I need to install some external software, such as Phobos. I've downloaded Phobos and I am using it by typing the path to it's executable like this:

`input/software/phobos-v3.3.12-mac/bin/phobos-mac-intel-64bit-start-from-command-line genome.fa output.txt`

This works quite well.

In order for TRAL to find Phobos, I changed TRAL's configuration file as explained here.

I've replaced `PHOBOS = phobos` with `PHOBOS = /Users/user/Documents/.../project/input/software/phobos-v3.3.12-mac/bin/phobos-mac-intel-64bit-start-from-command-line`

Then, I followed the tutorial Run and parse de novo repeat detection software and typed the following commands in Python 3 in the Terminal:

```
import os
from tral.sequence import sequence
from tral.paths import PACKAGE_DIRECTORY

proteome_HIV = os.path.join(PACKAGE_DIRECTORY, "examples", "data", "HIV-1_388796.faa")
sequences_HIV = sequence.Sequence.create(file = proteome_HIV, input_format = 'fasta')

```

So far so good. But the following command creates an error:

```
tandem_repeats = sequences_HIV[0].detect(denovo = True, detection = {"detectors": ["Phobos"]})

Did not find Phobos result file in /var/folders/z2/rt6qxntd3tg5bjmqlqnm0mdm0000gn/T/tmpvc20y3i4/001/PHOBOS/phobos.o

```

Note that if I use `PHOBOS` instead of `Phobos` in the command, I get the following error:

```
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/sequence.py", line 162, in detect
    **kwargs['detection'])[0]
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/repeat_detection_run.py", line 962, in run_detector
    Detectors(detectors, sequence_type)
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/repeat_detection_run.py", line 872, in Detectors
    FINDER_LIST.keys()))
Exception: Unknown TR detector supplied (Supplied: ['PHOBOS']. Known TR detectors: dict_keys(['HHrepID', 'Phobos', 'TRED', 'T-REKS', 'TRF', 'TRUST', 'XSTREAM']))

```

I am working on Mac OS X. I installed TRAL with `pip install tral`.

# Answer
TRAL 0.3.5 has bugs running and parsing Phobos. You can follow progress towards a fix on this issue: https://github.com/acg-team/tral/issues/16

Until the issue is fixed, you should remove Phobos from the repeat detector list in your tral config file (~/.tral/config.ini) and use a different repeat detector.

```
[sequence]
    [[repeat_detection]]
        DNA = TRED, T-REKS, TRF, XSTREAM

```
> 3 votes

---
Tags: sequence-annotation, software-installation, repeat-elements
---


# Question
Title: Change sample names while using cummeRbund
I want to analyze RNA seq data using R package cummeRbund. Using this package I have a problem changing sample names. When I run this command in cummeRbund

```
samples(cuff)
sample_index sample_name sample_name parameter value
1            1          q1        <NA>      <NA>  <NA>
2            2          q2        <NA>      <NA>  <NA>
3            3          q3        <NA>      <NA>  <NA>
4            4          q4        <NA>      <NA>  <NA>
5            5          q5        <NA>      <NA>  <NA>

```

It shows q1,q2,q3,q4 and q5 as sample names for the five samples but I want to change the sample names and put following names for the samples Batis, Blue\_silver, Chakwal\_50, Local\_white and Syn\_22 for q1,q2,q3,q4 ad q5 respectively. To accomplish this I ran following command:

```
> samples(cuff) <- c("Batis", "Blue_silver","Chakwal_50","Local_white","Syn_22")
Error in samples(cuff) <- c("Batis", "Blue_silver", "Chakwal_50", "Local_white",  :
could not find function "samples<-"

```

but this command gave error as shown above.

Also the BAM files have same names for all the samples as mentioned above.

How could I change the sample names to desired names in R using cummeRbund?

Or I need to re-run cuffdiff and change the names somewhere else?

# Answer
Don't know which file was used by cummeRbund for sample names but replacing q1,q2,q3 and q4 with Batis, Blue\_silver, Chakwal\_50 and Syn\_22 for all the files in the directory containing cuffdiff output worked for me.

I used sed language for this as shown below:

```
sed -i 's/q4/Syn_22/g' *

```

Or more generally you can use

```
sed -i 's/sample_name_shown/sample_name_to_change/g' *

```

After that I deleted the sqlite object created by cummeRbund and again created this object using **readCufflinks()** in R. This time sample names are changed as desired:

```
samples(cuff)
sample_index sample_name sample_name parameter value
1            1       Batis        <NA>      <NA>  <NA>
2            2 Blue_silver        <NA>      <NA>  <NA>
3            3  Chakwal_50        <NA>      <NA>  <NA>
4            4      Syn_22        <NA>      <NA>  <NA>

```
> 1 votes

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


# Question
Title: A basic question about phylogenetic inference with Bayes-MCMC
I am new to the world world of Bayesian-MCMC inference and have been working my way through various texts. I am having some difficulty understanding how MCMC works. I understand that it is a way of sampling from the posterior distribution. What is puzzling me, however, is how the probabilities of the posterior distribution have been calculated. That is, if the Markov Chain is sampling from the posterior distribution, that posterior distribution must have probability values on some axis. How were those calculated?

# Answer
In short, the posterior is a probability value between 0 and 1 which is calculated from Bayes rule, which In phylogenetic usage is given as the product of the prior probability and the likelihood divided by a normalising constant. The posterior probability allows you to state the probability of some hypothesis given the data at hand (as opposed to the likelihood parlance, which is probability of data given hypothesis).

In a phylogenetic setting, the sampler works (roughly) like this: propose some value of a parameter (topology, branch lengths, etc.) and use the specified prior and calculated likelihood for that parameter value given the data to calculate the unnormalised posterior. Next, propose a new value, and repeat the same step. Compare the two values, and if the posterior probability of the second is greater, accept this parameter value and repeat the cycle again. In this way the Markov Chain moves through posterior space until it reaches a reasonably stable optimal approximation of the posterior distribution. Here, we want to get as many samples as possible to maximise the accuracy of our posterior. The maximum a posteriori (MAP) value of the parameter and its credible intervals, etc, will come from the distribution of these samples.
> 2 votes

---
Tags: statistics, phylogenetics, phylogeny
---


title: The algorithm behind makeblastdb?

text:
# Question
Title: The algorithm behind makeblastdb?
> What is the algorithm/concept behind‚`makeblastdb`?

I know `makeblastdb` produces a `.phr`, `.pin`, and a `.psq` file, that represents an index or a data structure that is better optimised for searching a given query.

> What is exactly stored in these files?

Seems like `.phr` and `.pin` are storing some metadata (based on the size) and `.psq` the actual index. It can't be the table of all k-mers since the command doesn't take a value of k? (or maybe it is fixed to 3)

I could not find any direct discussion of the algorithm since google is flooded with people asking how to build such a database... Maybe someone can give me a link or hint?

# Answer
> 2 votes
From looking at this publication "Database indexing for production MegaBLAST searches" which although deals with MegaBLAST and it's indexing method, it would appear that this was rolled into the NCBI C++ toolkit as the command `makembindex` so I'm inclined to think the current implementation must have been derived from this strategy. Looking at the relevant part of the current NCBI C++ toolkit I find no reference to `makembindex` but only `makeblastdb` which I assume has superseded the older command: blastdb Directory Reference

Ulimtately the true details of the index generation lie in the C++ source code for `makeblastdb` can be found here. However following on from that I think the further details of the index generation may be in one of the following files which are includes for the main `.cpp` file:

`seqdb.hpp`

`build_db.hpp`

`seqdbcommon.hpp`

`writedb.hpp`

It would appear from `seqdb.hpp` that:

The index files have extension `.pin` or `.nin` and sequence files `.psq` or `.nsq`,

However it's not at all clear how the current implementation differs to the 2007 reference or where the `.phr` file is generated from! Some degree of reverse engineering the C++ source code will be needed!

---
Tags: blast, makeblastdb
---


# Question
Title: How to uninstall hmmpress from my conda environment?
I have installed `hmmpress` on my Mac OS X. It is here:

```
which hmmpress
/Users/user/anaconda/bin/hmmpress

```

It is available in my root environment:

```
hmmpress -h
# hmmpress :: prepare an HMM database for faster hmmscan searches
# HMMER 3.1b2 (February 2015); http://hmmer.org/
# Copyright (C) 2015 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Usage: hmmpress [-options] <hmmfile>

Options:
  -h : show brief help on version and usage
  -f : force: overwrite any previous pressed files

```

I don't remember exactly how I installed it but I thought that I used conda because it is in the `anaconda` folder. However, `conda search hmmpress` returns no packages. I thought it could be part of the `hmmer package`. Indeed, the hmmer package exists in the bioconda channel. However, if I try `conda uninstall hmmer`, it says:

```
conda uninstall hmmer
Solving environment: failed

PackagesNotFoundError: The following packages are missing from the target environment:
  - hmmer

```

Does anyone know how I can remove `hmmpress`? I thought of `rm -rf` the `/Users/user/anaconda/bin/hmmpress` but it does not seem like a clean solution. There might be dependencies elsewhere that I could also get rid of.

# Answer
> 1 votes
This time `conda uninstall hmmer` in the root environment did the job. The sequence of commands that led me there is shown below. Thanks Devon Ryan.

```
conda list | grep hmmer
hmmer                     3.1b2                         3    bioconda

conda uninstall hmmer
Solving environment: done

## Package Plan ##

  environment location: /Users/user/anaconda

  removed specs:
    - hmmer

The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    conda-4.5.9                |           py36_0         626 KB  conda-forge

The following packages will be REMOVED:

    barrnap: 0.9-0        bioconda
    busco:   3.0.2-py36_6 bioconda
    hmmer:   3.1b2-3      bioconda
    prokka:  1.12-3       bioconda

The following packages will be UPDATED:

    conda:   4.5.8-py36_1 conda-forge --> 4.5.9-py36_0      conda-forge
    openssl: 1.0.2o-0     conda-forge --> 1.0.2o-h470a237_1 conda-forge

Proceed ([y]/n)? y

Downloading and Extracting Packages
conda-4.5.9          |  626 KB | ################################################################################################################# | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

```

Checking that `hmmpress` is not available anymore:

```
hmmpress
zsh: command not found: hmmpress

```

---
Tags: software-installation, hmmer, bioconda
---


# Question
Title: Should PCA be standardized for gene expression?
This is a theory/good practice question more than a technical one. If samples are being plotted on a PCA projection of gene expression data, I'm wondering whether it is standard (and if so, why) to center and scale the PCs.

The reason I ask is that in this case the variables (genes) are being measured in the same scale, so I don't know if it would be necessary to center/scale data.

# Answer
> 4 votes
When the gene expression is scaled and centred you reduce the difference between genes.

Imagine you have one gene A that is highly expressed usually and has a standard deviation of 500 units compared to a gene B that is not much expressed and only have a standard deviation of 5.

In the scaled and centred genes both contribute the same because A usually is expressed 10000 and B is usually expressed 100 units. So, for both the standard deviation is 5% of their expression. Meaning that a variation in one is as important as in the other.

If not scaled (and centred) the first gene A will contribute more to the variation than gene B, because the expression variation in absolute numbers is bigger.

Both are used in publications, but I think the scaled and centred is more used, because the first dimension reflects "better" the differences between samples. Of course you can do both and select the one that is better to show what you want.

---
Tags: statistics, pca
---


# Question
Title: Omics data: How to interpret heatmap and dendrogram output?
How to interpret heat map and dendrogram output for biological data (omics) in words (when writing results and discussion)?

What should I consider (statistics behind?) and what is the best approach?

Here is one of my HM for proteomics data.

Script being used

```
col <- colorRampPalette(c("red","yellow","darkgreen"))(30)
png("HM1.png")
heatmap.2(as.matrix(df), Rowv = T, Colv = FALSE, dendrogram = "row", #scale = "row", col = col, density.info = "none", trace = "none", margins = c(7, 15) )
dev.off()

```

Heat Map output.

I would be greatful if you would give an example.

# Answer
> 3 votes
The dendrogram summarize the information of a group of values and sort them according to the similarity they have. It can be applied to both, samples and features.

The dendrogram allows to visualize features that are more similar together, usually revealing patterns that wouldn't have been seen otherwise.

In an article usually it is used something like: "we can see that X and Y are clustered together revealing that ... (they have something in common), while Z and A are clustered far"

# Answer
> 3 votes
You may be interested in reading up on heatmaps. For a history perspective (pre the biological introduction by Eisen et al. ) read *The History of the Cluster Heat Map* by Wilkinson and Friendly

---
Tags: r, clustering, proteomics
---


# Question
Title: A tool or webserver for building PSSM matrix
I have some protein sequences and I want to build a position-specific scoring matrix (PSSM) for them and then upload this PSSM to NCBI PSI-BLAST. I used CHAPS program for this pupose but uploading the output PSSM gave me an error in NCBI PSI-BLAST. Do you know any tool or webserver for getting PSSMs for a group of protein sequences which then can work in NCBI PSI-BLAST? This is a part of my PSSM file:

```
PssmWithParameters ::= {
pssm {
    isProtein TRUE ,
    numRows 28 ,
    numColumns 131 ,
    byRow FALSE ,
    query
    seq {
        id {
            other {
                accession "WP_000208753" } } ,
        inst {
            repr raw ,
            mol aa ,
            length 131 ,
            seq-data
            ncbieaa "MTTKRKPYVRPMTSTWWKKLPFYRFYMLREGTAVPAVWFSIELIFGLFALKNGPEAW
AGFIDFLQNPVIVIINLITLAAALLHTKTWFELAPKAANIIVKDEKMGPEPIIKSLWAVTVVATIVILFVALYW" } } ,
    intermediateData {
        freqRatios {
            { 0, 10, 0 } ,
            { 564418841, 10, -10 } ,
            { 0, 10, 0 } ,
            { 11768571, 10, -9 } ,
            { 185838265, 10, -10 } ,
            { 31496547, 10, -9 } ,
            { 3872857, 10, -8 } ,
            { 291750464, 10, -10 } ,
            { 128450763, 10, -10 } ,
            { 759856221, 10, -10 } ,
            { 359173937, 10, -10 } ,
            { 179865517, 10, -9 } ,
            { 14537895, 10, -8 } ,
            { 212921456, 10, -10 } ,
            { 220554141, 10, -10 } ,
            { 368516324, 10, -10 } ,
            { 319343525, 10, -10 } ,
            { 426173953, 10, -10 } ,
            { 463659523, 10, -10 } ,
            { 817322186, 10, -10 } ,
            { 811693964, 10, -11 } ,
            { 0, 10, 0 } ,
            { 227810064, 10, -10 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 768325726, 10, -10 } ,
            { 0, 10, 0 } ,
            { 1425562, 10, -8 } ,
            { 372685285, 10, -10 } ,
            { 466727421, 10, -10 } ,
            { 185741084, 10, -10 } ,
            { 427328646, 10, -10 } ,
            { 122594965, 10, -10 } ,

```

# Answer
If you take a look at my answer in this BioStars post, you can generate a PSSM using `AlignIO` in Biopython:
> 2 votes

# Answer
You can use biopython's code for generation of PSSM. Here is the link for the source code which can help you.
> 1 votes

# Answer
I'm not 100% sure what you're trying to do with regards to your search in the first place, but I recommend using Sean Eddy's program Jackhammer it's similar to PSI-BLAST but uses Hammer HMM underpinnings.

As input it can take a single sequence, a multiple sequence alignment, or an existing HMM and then allow you to search in an iterative model building way. You can also download an HMM too at any stage in the iterative search, as well as the results / alignment.
> 1 votes

---
Tags: blast, phylogenetics, phylogeny
---


# Question
Title: Generating list of known fusions for Arriba from COSMIC
The RNA-Seq fusion detection program Arriba features the ability to white-list known gene fusions, specifically the authors recommend using COSMIC as a source of these here. Specifically they state:

> A comprehensive list of known fusions can be obtained from CancerGeneCensus in the section titled "Complete Fusion Export". *Depending on the gene annotation that is used to run arriba, some gene names need to be adjusted*.

Emphasis is mine to highlight the issue here!

The state the format of the input file they require to be:

> The file has two columns separated by a tab. Each line lists a pair of genes. The order of the genes is irrelevant. arriba searches for both genes as the 5' end and the 3' end of a fusion. Lines starting with a hash (#) are treated as comments. Optionally, the file can be gzip-compressed.

However the format of the download from COSMIC is somewhat different:

```
Sample ID       Sample name     Primary site    Site subtype 1  Site subtype 2  Site subtype 3  Primary histology       Histology subtype 1     Histology subtype 2     Histology subtype 3     Fusion ID       Translocation Name      Fusion Type     PUBMED_PMID     ID_STUDY
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16784981
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1271    CCDC6{ENST00000263102}:r.1_535_RET{ENST00000355710}:r.2369_5659 Observed mRNA   16784981
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1272    CCDC6{ENST00000263102}:r.1_535_RET{ENST00000355710}:r.2369_5659 Inferred Breakpoint     16784981
XXXXXXX XXXXXXX thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16784981
YYYYYYY YYYYYYY thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16728576
YYYYYYY YYYYYYY thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1216    PAX8{NM_003466.2}:r.?_PPARG{NM_015869.2}:r.?    Inferred Breakpoint     16728576

```

*Note: numeric samples IDs have been deliberately replaced with letters here.*

Using this file, as is, will cause arriba to spurt out many warnings that the gene names in the first two columns are not recognised:

`WARNING: unknown gene in known fusions list: XXXXXXX`

From the looks of the rest of my arriba output the currency of gene IDs its using would appear to be gene symbols taken from the `gencode.gft` I'm supplying it.

I'm thus assuming the format arriba requires for a known gene list to be:

```
SYMB1   SYMB2
SYMB3   SYMB4

```

Where the record separator is a tab character `\t`

My question is how can I reprocess the COSMIC "Complete Fusion Export" `.tsv` file to the format required for arriba to act as a white-list?

# Answer
I managed to solve this with a one-liner:

```
grep -v 'Sample ID' CosmicFusionExport.tsv | cut -f 12 | grep -v -e '^$' | perl -lne '/^(\S+?)\{\S+\d+}\S+_(\S+?)\{/; print "$1\t$2"' | sort -u > CosmicFusionList.tsv

```

The various parts of the piped command are as follows:

1. Remove header:`grep -v 'Sample ID' CosmicFusionExport.tsv |`
2. Extract 12th field (Translocation Name):`cut -f 12 |`
3. Remove blank entries:`grep -v -e '^$' |`
4. Extract gene symbols from 12th column:`perl -lne '/^(\S+?)\{\S+\d+}\S+_(\S+?)\{/; print "$1\t$2"' |`
5. Remove duplicate lines:`sort -u > CosmicFusionList.tsv`
6. Optionally the final two stages could be:`sort -u | gzip -c > CosmicFusionList.tsv.gz` to make use of file compression (for what it's worth on a 3.9 kB file)

Since the authors of arriba have neglected to document this necessary data munging step, I've posted it here in Q and A format in the hope that it will be useful to others, and other solutions can be proposed.
> 1 votes

---
Tags: rna-seq, fusions
---


# Question
Title: How to force pilon to correct a reference with low coverage mapped reads?
I have mapped some reads to a reference and want to change SNPs in the reference to SNPs found in the reads, even at low coverage. For example, here is a screenshot of one reference in IGV that I would like to correct.

When I use pilon to try to correct the reference, however, no bases are ever changed, even if using the option `--mindepth 1`.

Here is the call that I am using:

```
java -Xmx100G -jar {params.pilon}\
  --genome {input.ref}\
  --frags {input.bam}\
  --output {params.outname}\
  --changes --fix all,amb,breaks,bases\
  --threads {threads}\
  --mindepth 1

```

How can I use even low-coverage information to correct my reference?

# Answer
> 1 votes
Pilon will not consider fragments that are not proper pairs if the `--frags` option is used to pass the alignment `.bam` file. The reads mapped in the above IGV screenshot are unpaired.

Instead, if you want Pilon to consider *all* reads that map as legitimate and usable for correction, pass the `.bam` file with the `--unpaired` option instead.

To above corrected example would be:

```
java -Xmx100G -jar {params.pilon}\
  --genome {input.ref}\
  --unpaired {input.bam}\
  --output {params.outname}\
  --changes --fix all,amb,breaks,bases\
  --threads {threads}\
  --mindepth 1

```

---
Tags: sam, reference-genome, read-correction
---