# Question
Title: How to perform DE analysis for each sample
I am new to R and biocondunctor. I have the normalized expression values for 20 samples for a disease and for 10 controls. I wanted to get the differential expressed values for each sample with all 10 controls.

How could I do that with limma.

Any help will be appreciated.

Edit: What I want to achieve is single-sample t-test between the control population and the single sample values to get DE values

# Answer
> 1 votes
Treat your disease samples as individual groups and then follow the normal routine in limma to use contrasts to compare two group (the control group versus individual disease samples). Note that you need to think long and hard about what these results then mean. Normally we compare groups because then the results should generalize to other samples. That will not be the case when you look at individual disease samples, so I suggest you do some clustering of the samples to see if you perhaps have any coherent clustering of the disease samples into subsets that could then be used as more useful groups in DE testing.

---
Tags: differential-expression, microarray, limma
---


# Question
Title: How can I find out if my gene of interest is duplicated?
To identify homologous genes, I performed tBLASTn analyses against my genome of interest resulting in sequences of a closely related species. Using this method, I found two (and in some cases up to five copies) of my different genes of interest in the target genome.

*How can I determine if these gene copies correspond to true duplication events?*

I aligned the newly found sequences to the reference genes and I found good alignments in more than one region (20% - 50% identity), but I am unsure if this is enough to consider them as duplicates.

Is there a **threshold** to help me determine if they are indeed duplications? How can I distinguish between a true duplication event and a pseudogene or a paralogous gene? The target genome belongs to a polyploid organism, although the genome assembly type is listed as haploid.

Any help will be greatly appreciated.

---

Thank you all very much for your responses, I am now beginning to understand how to find out if my genes are duplicated.

I performed some of the analyses that you suggested and I would very much appreciate your input on the results.

First, I followed @terdon's advice and I took the best hit for one of the query (source) genes in the target genome and I used this sequence as a query against the query (source) genome, i.e., I performed a reciprocal BLAST analysis (in my case it was a reciprocal BLASTp analysis). The first hit was my original query gene, which had an E-value of 0.0. For other sequences the reciprocal BLASTp hits were more variable, ranging from 3e-42 to 2e-154, but the first hit was always the query sequence. I also looked at the protein domains and they are conserved, i.e., both the query and the target sequences have the same domains. As for the exonic structure, it is fairly similar, for e.g., the query has 12 exons while the target sequence has 10. Exons are similar in length. So, in summary, when the potential duplicates from my target genome are blasted against the original source genome, the best hit is the original query sequence.

As for @zhanf's suggestion, I did a basic phylogenetic analysis using Clustal Omega (neighbour-joining tree without distance corrections) using four sequences: the original query sequence and the three gene copies (seq\_5-1, seq\_5-2 and seq\_5-3) that I found in the target genome. This is the result:

I also compared each of the sequences against the original query:

* seq\_5-1 vs query: 34% identity (query length 780 aa, seq\_5-1 length 400 aa)
* seq\_5-2 vs query: 42% identity (query length 780 aa, seq\_5-2 length 580 aa)
* seq\_5-3 vs query: 40% identity (query length 780 aa, seq\_5-3 length 490 aa)

So in summary, these sequences would appear to be paralogues, right?

Again, thank you very much for your help, I greatly appreciate it!

Leah

# Answer
> 2 votes
You might want to have a look at What is the difference between orthologs, paralogs and homologs?. Paralogs are duplicates by definition, and whether a gene is a pseudogene or still active has no bearing on whether it comes from a duplication.

In any case, to figure out what's going in in your target genome, you should use the target's genes. So, first get the best hit for each of your query genes in the target genome. Then, to investigate its evolutionary history, use that sequence (the one from the target genome) as query against the target genome. That will let you assess how similar the target's hits are to one another and that's how you can attempt to figure out their evolutionary history.

Note that 20% identity is very low for closely related species. While there is no golden threshold, I would expect recent duplication events (and in closely related species, a duplication event that is in one genome and not the other will be recent) to produce much more similar sequences. In any case, you can't do this based on sequence identity alone.

Are the protein domains of the gene products conserved? Is the exonic structure of the gene conserved? You should also try taking the putative duplicate of geneA from genome2 and blast it against genome1 (your original source genome). Is the best hit the geneA from genome1, the query you used to find the putative duplicate in genome2? Or is there a better hit for it against another gene?

These are all questions to investigate and take into account when reconstructing the evolutionary history. You can't just look at identity and use some threshold to decide this sort of thing.

If you blast the potential dupicate from your target genome against your original source genome, is the best hit the duplicate?

# Answer
> -1 votes
I don't have enough reputation to comment but if a sequence is duplicated then there has obviously been a duplication event. Doesn't matter how scientists call it now.It can be older or more recent but it's yet a duplication.

---
Tags: sequence-alignment, gene, phylogenetics, phylogeny, duplications
---


title: public multi-modal single-cell data

text:
# Question
Title: public multi-modal single-cell data
There is a scRNAseq Bioconductor package with a few different example scRNA-seq datasets. Are there any R packages that offer multiple modalities of single-cell data? For example, hashtags or ADTs or ATAC-seq. I am trying to find something that can be used for testing or for a quick tutorial. Preferably something that is on CRAN or Bioconductor.

# Answer
These tutorials on Seurat multimodal data and the wrapper Seurat data are easy ways to start.

The wrapper has some cite-seq data preinstalled making it easy to work with benchmarked data sets

If you are using Seurat u can just as well have a look at Signac to start working with some sc-atac seq
> 5 votes

---
Tags: scrnaseq, single-cell
---


# Question
Title: Dmpfold error "No such file or directory: 'best_qdope.pdb'" when running "run_dmpfold.sh"?
I am running DMPfold, "seq2maps.sh" ran successfully, but whrn running "run\_dmpfold.sh", it gives me this error:

```
   Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
seed =  30603
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/run_dmpfold.sh: line 177: c600m2-tail-before-lysin-sequence_[1-9]*.pdb: No such file or directory
No atoms read from PDB file
0 models read from PDB file
Outliers:

Traceback (most recent call last):
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-softmax/pytorch_dmp_iterdistpred_cb.py", line 82, in <module>
    main()
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-softmax/pytorch_dmp_iterdistpred_cb.py", line 55, in main
    with open(sys.argv[3], 'r') as refpdbfile:
FileNotFoundError: [Errno 2] No such file or directory: 'best_qdope.pdb'
Traceback (most recent call last):
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-hb/pytorch_iterdmp_hbpred.py", line 75, in <module>
    main()
  File "/media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/nn/dmp-hb/pytorch_iterdmp_hbpred.py", line 53, in main
    with open(sys.argv[3], 'r') as refpdbfile:
FileNotFoundError: [Errno 2] No such file or directory: 'best_qdope.pdb'
No valid current contacts.

```

Can you please let me know how this can be solved?

the comand line is:

```
run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./dmpfold-lysin 10 500 &> dmpfold-log.log

```

Thanks!

# Answer
The `target` variable, which is set to the first part of the sequence file name, should be `sequence` according to your given command line. However the error message thinks it is `c600m2-tail-before-lysin-sequence`. Have you edited the command line for simplicity?

I'm a little worried about the floating point errors, which are almost certainly from CNS as there are 3 CNS calls and the error shows up 3 times. I would test your CNS installation to see if everything is okay - maybe try some examples from their website. An error there would explain why no output is generated.

You can ignore any errors below line 5 of the above program output for now, as they are follow-on errors. In particular the question title isn't the 'real' error here, the real error is that no PDB files are generated on line 160.

By the way, I wouldn't use `10 500` as parameters for iteration number and number of models. It will take ages to run. I'd just use the defaults (omit the arguments), which we find gives good results.

Finally, I'll point out that the DMPfold repository has an issue tracker where you can raise any future issues. Though we may as well finish debugging this one here since we've started!
> 1 votes

# Answer
thanks! I figured out the solution, it was the cns. I downloaded the binaries instead of compiling from source and it works!
> 0 votes

---
Tags: protein-structure
---


# Question
Title: DMPfold: is it suitable to get 5000 models to compare with rosetta
I ran this command for 10 iteraions and 500 models.

```
run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./dmpfold-lysin 10 500 &> dmpfold-log.log

```

I am looking to generate 5000 dmpfold models to compare rosetta results with dmpfold; is it reasonable to do this? and, 5000 models can be generated by using "10 5000" parameter or "10 500" when running dmpfold as stated in the command line? Thanks!

# Answer
DMPfold tends to produce low-variability models - it produces a small number of distinct conformations, and either the right one is in there or it isn't. This arises from the way the models are generated from distance constraints. Rosetta, by contrast, produces many different conformations which is why 5000 or more models tend to be generated. This arises from Rosetta's method of assembling small fragments, which allows for many more combinations.

You *could* generate 5000 models with DMPfold. Running with `10 5000` would give you 5000 models as output from the last iteration. Running with `10 500` would give you 500 models from the last iteration but you could also use the 500 models from each of the previous iterations (saved as `ensemble.1.pdb` etc.) to get 5000 models overall.

However I wouldn't recommend it, it sounds like a waste of disk space and compute time. Just run with default parameters `3 50` and use the final models saved as `final_*.pdb`.

For some numbers on how DMPfold compares to larger numbers of Rosetta models, and how much better DMPfold gets when you use more models, see Table 1 of the paper.
> 3 votes

---
Tags: protein-structure
---


# Question
Title: How to make a robust amino acid phylogeny with adequate robustness support
This is the second part of this thread: Duplicated genes and genome fragmentation

As @Michael G. suggested, I did a pseudo-annotation of sequences 5-2 and 5-3 (i.e., I joined both sequences into one) and I did the phylogenetic analysis again. Here is the tree that I obtained:

Both of my sequences are shown in a red rectangle. The length of the branch seems to have decreased a bit. Also, I looked at the branches near sequence 5-1. One of them is the original query sequence from species ***A*** and the other one is a hypothetical protein from the same species (probably related to the query sequence, but I couldn't find any domains because it is only a few amino acids long). I still need to find the root of the phylogeny but I am not sure how to do this...

# Answer
To finish you need bootstrap support using a model-based tree building algorhithm, I use RAxML, specifically raxmlHPC (easily downloadable and compiles on Linux and OSX). The codes are quite complicated and I've given them below.

A robust maximum likelihood tree is,

```
raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup-p 10000 -s alignment_file.phy

```

To make a bootstrap tree,

```
raxmlHPC -b 459 -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -s /path_to_alignment.phy -N 500

```

Make a consensus tree of the bootstraps,

```
raxmlHPC -f b -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -t /path_to_master_tree(best_tree) -z path_to_bootfile

```

You require bootstrap support \>80% and please repeat this with and without 5-2/5-3 (it still looks long)

You need a fairly robust desktop, dedicated for around 24 hours to complete one calculation. You have enough results to draw a conclusion.

Viewing the tree, FigTree is easy.

**Root**

Check the work of Cavalier-Smith, he heavily focused on basal protists to find a species which evolved, using external evidence, before any other member of your group. One of the ingroups may have sufficient evidence to say it evolved before the rest and that is your root, which you select as an outgroup.
> 1 votes

---
Tags: phylogenetics, phylogeny, duplications
---


# Question
Title: Plotting distance tree from blastn output
I'm trying to plot a simple distance tree of my blastn output with nj (like the tree view on NCBI). From what I understood, what I think I should do is

1. extract all the hsps from each alignment
2. re-align them using MUSCLE
3. generate a distance matrix from the msa file
4. plot the tree using nj

However, there is a couple problems I am not sure about:

1) When there are multiple hsps in one alignment then which one do I include in the tree?

2) If two sequences map to two different parts of the query sequence and do not overlap, then what should I do?

# Answer
> 1 votes
Specific questions

1. When there are multiple hsps in one alignment then which one do I include in the tree?

You can included them all because NJ is not computer intensive, so running loads is fine. You can remove identical sequences if you want either programmatically (these days I use pandas Python) or something like RAxML. Its honestly not critical for an NJ tree (and in fact the more the better) providing you can read the labels. A tree viewing program such as FigTree allows you to collapse clades to make the tree more readable to users.

2. If two sequences map to two different parts of the query sequence and do not overlap, then what should I do?

Hmmmm.... good question. Do 3 alignments one with one e.g. 5' fragment, the other with the 3' fragment and attempt a joint alignment. Under other non-pairwise tree buidling algorithms this would be disasterous, NJ gets around the problem because it is pairwise. HOWEVER, ask yourself whether your feel the resulting tree is meaningful, I would argue not because there are no homologous sites between the 5' and 3' partial sequences and phylogeny models point mutations. NJ will give you a tree however and will even bootstrap it.

Generic questions, you can actually do this as one of the new features of NCBI's Blast

* Go to Blast here, https://blast.ncbi.nlm.nih.gov/Blast.cgi
* Enter your sequence into the box (it doesn't accept PDB codes alone)
* Enter the protein database - when I first did this calculation I used SwissProt, thinking there would be alot of sequences - I then used "nr"
* Under the algorithm parameters enter "50" (default is too many)
* Hit "Blast"
* Once the search is complete at the top of the page are the hyperlinks: "Other reports:
* Search Summary \[Taxonomy reports\] \[Distance tree of results\]"
* Click on "Distance tree of results"
* The tree page will load automatically, automatically aligning your sequences and producing, in this case a parsimony based tree, but there is also the option of a nj tree (recommended) ....
* Click "Tool", "Download", "PDF" ...

You can generate a distance and matrix and/or perform a neighbor-joining phylogeny in MEGA10 and using Timura 3-parameter model is a good model. The bootstrap option is important. Its very easy to use, the format is a bit weird

```
# MEGA
Title: Stuff
# Sequence1
ACTAGACGT
# Sequence2
ACCTTAGGA

etc ...

```

---
Tags: ngs, blast, blastn
---


title: Read alignment using Bowtie2

text:
# Question
Title: Read alignment using Bowtie2
So this is related to CRISPR-CAS9. I am working with off-target predictions for my thesis and was looking at all scientific papers related to CRISPR. I found one and decided to use their datasets. The datasets were not available directly and the method used to obtain them was specified which I had mentioned below:

They have around 30 sgRNAs which are each 23bp long and are aligning the sgRNAs to human reference genome with a tolerance of upto 6 mismatches using Bowtie2.

In order to access the exact dataset, I thought I need to replicate the method they used. To do so, I collected the sgRNAs from their paper and tried to use Bowtie2 to align them against human genome with upto 6 mismatches.

Since my laptop does not have the capacity to perform this operation (it has 8GB RAM), I am searching for interfaces that provide access to Bowtie2. I found Galaxy to be useful. Does anyone know of any other tool that provides an option for Bowtie2?

Thanks!

# Answer
> 0 votes
One possible workaround would be to blast your sequences against the genome. Then, when you have found where they align, you can extract those regions only into a new fasta file (make sure to take a few kb around each target region) and then use that fasta file as the "genome" you pass to bowtie2.

That would give bowtie a far smaller file to work with and should reduce the memory needed significantly.

---
Tags: sequence-alignment, bowtie2, crispr
---


# Question
Title: DMPfold: Understand the run of DMPfold "run_dmpfold.sh"?
I ran dmpfold using this command:

```
run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./dmpfold-lysin 10 500 &> dmpfold-log.log

```

It produced files like this:

sequence\_sub\_embed\_1.pdb ..... sequence\_sub\_embed\_500.pdb

And, now I see some more files like this being produced:

sequence\_1.pdb .... sequence\_2.pdb .....etc

Can you please let me know what is dmpfold doing here? And can I stop dmpfold in between when desired number of models are generated? And, while it is running, can I move the previously made pdbs to new location becuase of disk space issues or all pdbs have to be in the "dmpfold-lysin" as stated in the command-line above? Thanks.

# Answer
From the docs:

> The final model is `final_1.pdb` and other structures may or may not be generated as `final_2.pdb` to `final_5.pdb` if they are significantly different. Many other files are generated totalling around 100 MB - these should be deleted to save disk space if you are running DMPfold on many sequences.

These other files include intermediate CNS files (`sequence_sub_embed_1.pdb` etc.), output files before clustering (`sequence_1.pdb` etc.) and saved ensembles of previous iterations (`ensemble.1.pdb` etc.).

Let the run finish or you won't get the final output files. If you want to use all the generated structures, use the `ensemble.n.pdb` file from the last iteration (i.e. the highest number). You can also use the other ensemble files if you want to include models from previous iterations.

I would strongly recommend using default "3 50" parameters until you get things working. Don't move files during the run. You really don't get much benefit of going beyond 50 models per iteration, trust me. We tried it relatively extensively!
> 1 votes

---
Tags: protein-structure
---


title: Fast and reliable alternatives to blast

text:
# Question
Title: Fast and reliable alternatives to blast
After some unexpected results (and previously reported) I heard that there are other tools for finding similar sequences besides blast that are faster and more accurate.

I only found hmmer, but I don't know if it is faster and more accurate.

Is there any review comparing such tools? A quick google scholar search didn't report anything useful.

I want to search with short DNA sequences (~400), mainly interested for taxonomic imputation so I'll be looking for similar sequences and not distant homologs. I won't try to model genes and the sequence will be unspliced as it is from direct sequencing of a a region of the 16S gene.

# Answer
> 4 votes
Besides generic nucleotide aligners, there are also more specialized tools for the alignment of 16S amplicon sequences, e.g. SINA (article, software) which is part of SILVA

# Answer
> 4 votes
As I understand, the software tool Lambda is a viable, yet lesser known alternative to BLAST in the context of taxonomic classification of NGS data.

# Answer
> 0 votes
Why not shift the parameter inputs of Blast to exclude the species/subspecies which are being queried? Happy to show you how to do it, but that will pick up homologues. Alternative parse the blast output for e.g. 20000 outputs using E or identity within the threshold required (that's the way I'd do it).

# Answer
> 0 votes
All tools I’m aware of will be more-or-less as fast as blast. With kmer-based methods, the investment will be the up-front kmer database build.

Not sure if you’ve tried this, but NCBI has a 16S microbial blast database (lots of species, 16S rRNA only) just for this purpose. No point to search millions of nts for the 16S sequence every time if you know it’s SOME 16S sequence, so limit the search to only16S sequences.

I’d recommend you build a local blast database from ftp://ftp.ncbi.nih.gov/blast/db/16SMicrobial.tar.gz

BTW, another answer says blast won’t give you the alignment. Blast actually will give you alignments, you just have to ask it for that output format.

# Answer
> -2 votes
You can give a try to this tool: dnaservic.es It's actually available experimentaly. It is a real local aligner, not just a search tool like BLAST. It really returns an alignment. It is usually slower than BLAST, but bot much slower in most cases. In contrast it usually returns alignments with better scores and not only partial alignments.

---
Tags: sequence-alignment, blast, literature-search
---


# Question
Title: How to I get the infix expression for the body in a SBML function defintion using python libsbml?
How do I translate the body of a function definition in SBML into infix using python libsbml?

For example after loading the model using:

```
reader = libsbml.SBMLReader()
doc = reader.readSBMLFromString (r.getSBML())
model = doc.model

```

# Answer
> 1 votes
To do the translation you need to access the function definition. First, make sure there are function definitions in your model by using:

```
if model.getNumFunctionDefinitions() > 0:
   # continue

```

To get the infix of the body for the ith function definition use the following code:

```
f = model.getFunctionDeinfition(i)
body = f.getBody()
print (libsbml.formulaToL3String(body))

```

As a side note you can also get the jth argument for the function definition using:

```
f = model.getFunctionDeinfition(i)
arg = f.getArgument(j)
print (libsbml.formulaToL3String(arg))

```

---
Tags: sbml
---


# Question
Title: can we show the orig.ident in separate windows by using Seurat
I'm curious about can we show the orig.ident in sprat windows by using Seurat? Here are my code and the fig I got. Is it possible that I can separate each orig.ident into a single panel? `DimPlot(whole,group.by = "orig.ident", reduction = "umap")`

# Answer
The `DimPlot()` function of the new version of Seurat, `Seurat v3` has a `split.by` argument, which splits the plot based on the levels of the variable provided. You basically need to provide a column from the `@metadata` slot:

```
DimPlot(whole, split.by = "orig.ident", reduction = "umap")

```

If you are using `Seurat v2`, you can have the `DimPlot()` function return the `ggplot` object created and then use `facet_wratp()` or `facet_grid()` to generate facets based on your variable of interest.
> 1 votes

---
Tags: single-cell
---


# Question
Title: How to best detect the "peaks" in RNA-seq data that are not assigned to any gene?
I encountered that many reads from single-cell RNA seq data were lost in the analysis because not assigned to any gene (genome: galgal6). I am trying to find an approach than could give me all the "peaks" in this data, that is to say regions where there is a high density of reads. Once I get there, I would extract "peaks" that are not assigned to any gene, and get the distance to the closest genes.

I am thinking of using first `bedtools genomecov` to get the "peaks", and then `bedtools intersect` to get all reads assigned to a gene. Finally, I extract reads that are not in the output of `bedtools intersect` to get a new file with all my reads from the peak, that are not assigned to any gene.

Does it make sense this way ? Can you see any other way of doing ?

Here are some tiny data:

```
> cat reads.bed
chr9    505479  505498
chr9    508014  508037
chr9    514603  514633
chr9    529519  529540
chr9    529519  529540
chr9    529519  529540

> cat tiny.galGal6.chrom.sizes
chr9    24153086

> bedtools genomecov -bg -split -i reads.bed -g tiny.galGal6.chrom.sizes
chr9    505479  505498  1
chr9    508014  508037  1
chr9    514603  514633  1
chr9    529519  529540  3

```

The last line would correspond to a peak. However, another difficulty is how should I define the thresholds ?

Any help would be more than welcome.

# Answer
> 5 votes
Basically what you're discovering is that there are unannotated expressed features, so your task isn't really finding peaks, but rather finding novel expressed transcripts. For that, you can use stringTie or similar programs. Ensure that the BAM file you give to stringTie have all of the cells together.

# Answer
> 2 votes
Many pipelines have been set up for this such as software developed for the FANTOM Project at RIKEN. The software for this is available from our website. These were set up for analysis of CAGE sequencing data for exactly this purpose of calling peaks of previously unannotated gene expression. You can also find data for many species (including chicken) from FANTOM5 where this has already been done.

---
Tags: rna-seq, scrnaseq, ngs, bedtools, coverage
---


# Question
Title: Oxford Nanopore mapping Quality and sequencing error
I'm analyzing some Oxford Nanopore sequencing reads, with several articles stating the sequencing error is between 5% and 15%. Since its 1D and not 1D^2 sequencing i would assume the sequencing error is closer to 15%.

I'm using Minimap2 aligner to align the reads to the reference genome, and then im choosing those reads with a mapQ of 60, to analyze.

Where i know with the phred score mapQ of 60 means 99.9999% of the bases are called correctly.

So my question is, if the mapQ is 60, then does it really matter if my sequencing error is 15%? Or how is the sequencing error correlated with mapQ? Becuase i guess you could still obtain a good mapping quality with a few mismatches?

Or am i completely off target.

Thanks for your time and help.

# Answer
> Where i know with the phred score mapQ of 60 means 99.9999% of the bases are called correctly.

That's not at all what the MAPQ means, it means that the sequence is very likely to have originated from that place in the genome, whether the underlying sequence is correct is another question all together.

For long nanopore reads, sequence quality is not well correlated with MAPQ. MAPQ will tend to be either very low (0 or 1) or very high (60) since you tend to have kb+ reads and those will largely be uniquely assignable in the genome. The sequencing quality of nanopore reads is only an issue when you care about something at the base level, such as SNP calling or finding exact splice junctions. Even then, you mostly need sufficient coverage to get around this.
> 3 votes

---
Tags: nanopore, read-mapping, quality-control, errors
---


# Question
Title: In-sample and across samples normalized expression
I want to get the expression data that is ***in-sample normalized*** like FPKM and also ***across samples normalized*** as obtained using DESeq2 or else.

What I am currently doing is that I first normalize the data across samples (using DESeq) and from the resultant expression I calculate the FPKM. Does it make sense or am I missing something here?

This is how I am doing (not the exact code but the idea)

```
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), ~ condition)
dds <- estimateSizeFactors(dds)
norm_data <- counts(dds, normalized=TRUE)

foreach sample
{
    foreach transcript
    {
        FPKM = (Normalized read count * 10^9) / (transcript length * total mapped normalized read count)
    }

}

```

# Answer
Yes, this is a standard way of obtaining RPKM/FPKM/CPM values for plotting. Not that you do not need to use a `for` loop for any of the computations in R. You have a matrix of normalized values and things like transcript length are constant across samples (at least unless you're using something like salmon...although then you'd have TPMs to begin with).
> 1 votes

# Answer
You might want to check the EDASeq() package, which attempts to address both gene-specific and sample-specific effects (I believe this is what you are trying to achieve). From their vignette:

> Following (Risso et al. 2011), we consider two main types of effects on gene-level counts: (1) within-lane gene-specific (and possibly lane-specific) effects, e.g., related to gene length or GC-content, and (2) effects related to between-lane distributional differences, e.g., sequencing depth. Accordingly, withinLaneNormalization and betweenLaneNormalization adjust for the first and second type of effects, respectively. **We recommend to normalize for within-lane effects prior to between-lane normalization**.
> 0 votes

---
Tags: rna-seq, normalization, reads, fpkm
---


title: Interpreting this PCA plot for RNA-seq

text:
# Question
Title: Interpreting this PCA plot for RNA-seq
I have RNA-seq from two sequencing batches; Lab technician says that he has run the RNA expression quantification two times in bathes 1 and 2 for example `tumor 1 in batch 1 and tumor 1 in batch 2 , normal 2 in batch1 and normal 2 in batch 2`. This is my design for DESeq2

```
> head(mycols)
        condition batch
N_1_305         N     1
N_1_310         N     1
N_1_337         N     1
N_1_353         N     1
T_1_305         T     1
T_1_310         T     1
> tail(mycols)
        condition batch
T_2_337         T     2
T_2_338         T     2
T_2_344         T     2
T_2_346         T     2
T_2_349         T     2
T_2_353         T     2
>

```

I got this PCA plot

And this is biplot of samples

In PCA plot I am seeing for instance , `T_1_337 (batch1) has been placed too close to T_2_337 (batch2)`

Then I used svd for detecting hidden batch

Does this mean that there is no big batch effect between experimental runs and I can concatenate the fastqs from both batches for each sample (technical replication) or collapse technical replicates afterwards?

Please help me to interpret these Thank you

---

> EDITED

---

Sorry @Devon, I have 4 lanes for each samples (paired end) in each of experimental runs; For concatenating fastq files can I do like this ?

```
cat fastq1_lane1_batch1 fastq1_lane1_batch2 fastq1_lane2_batch1   fastq1_lane2_batch2  fastq1_lane3_batch1 fastq1_lane3_batch2 fastq1_lane4_batch1 fastq1_lane4_batch2  > fastq1

cat fastq2_lane1_batch1 fastq2_lane1_batch2 fastq2_lane2_batch1   fastq2_lane2_batch2  fastq2_lane3_batch1 fastq2_lane3_batch2 fastq2_lane4_batch1 fastq2_lane4_batch2 > fastq2

```

# Answer
> 3 votes
Yes, you can safely concatenate the technical replicates. Odds are good that these are even the same libraries just sequenced twice, so even labeling them as replicates is a bit of a stretch. As an aside, it would be surprising if you actually had a batch effect in a situation like this. You will commonly see sequencing facilities just sequence a given sample a second time if it didn't get the number of reads requested by the client to begin with. They then just have the client concatenate the files.

# Answer
> 1 votes
Prepping the RNA on different days, or making Illumina libraries on different days, or having different technicians handle different samples; that can lead to batch effects. Running samples on two different days does not cause a significant batch effect, as you can plainly see in your PCA. You should just combine the fastqs.

---
Tags: r, rna-seq, deseq2, pca, batch-effects
---


# Question
Title: Download protein structure from gene name
I have a list of gene name, I want to download and save each protein 3D structure corresponding to each gene name.

# Answer
> 1 votes
There may also be a misunderstanding, from your previous post. To reiterate, I don't use Ensembl, which is the source of your gene names.

On Genbank there is an entire "Structural" database for every structure that is on PDB. When Genbank hold the protein sequence that was used to generate the structure Genbank only use PDB codes. For example "5LBV\_A" is both a PDB code and a Genbank code. So this means PDB entry 5LBV chain A. Thus its easy to download the crystal because you simply ask PDB for 5LBV. To prove it,

```
efetch -db protein -format gpc -id "5LBV_A" | xtract -pattern INSDSeq -group INSDSeq -element INSDSeq_sequence

```

OUT ircigvsnrdfvegmsggtwvdvvlehggcvtvmaqdkptvdielvtttvsnmaevrsycyeasisdmasdsrcptqgeayldkqsdtqyvckrtlvdrgwgngcglfgkgslvtcakfacskkmtgksiqpenleyrimlsvhgsqhsgmivndtghetdenrakveitpnspraeatlggfgslgldceprtgldfsdlyyltmnnkhwlvhkewfhdiplpwhagadtgtphwnnkealvefkdahakrqtvvvlgsqegavhtalagaleaemdgakgrlssghlkcrlkmdklrlkgvsyslctaaftftkipaetlhgtvtvevqyagtdgpckvpaqmavdmqtltpvgrlitanpvitestenskmmleldppfgdsyivigvgekkithhwhrsgstiggpfeddddkagwshpqfekgggsgggsgggswshpqfek

What I cannot easily do, is go from a generic NCBI Genbank protein code to PDB code.

However, PDB uses Uniprot codes for proteins, so you may be able to use a Uniprot code to get hold of the PDB code, or Ensembl codes may, like Genbank, use PDB accessions.

So in summary perhaps you need to give more detail?

---
Tags: gene, protein-structure, pdb
---


# Question
Title: TPM or rlog(CPM) for comparing expression?
I want to see the expression of a gene in a group of patient amongst the entire cohort using my RNA-Seq data. While I can do a differential expression analysis with limma or DESeq2, I want to see how much each sample from my cohort expresses the gene.

The plan is to plot a waterfall plot (as defined in this paper). However, when I used TPM and rlog(CPM) they gave me very different looking graphs (Not just the shape, which is expected, but the ranking of expression). My concerns are, for TPM, the batch effect with two or three lots in one cohort, and, for rlog(CPM), the results don't match with prediction but is adjusted for batch effect. Which one would be a right or better choice?

Upper one is with plot with. Lower one is plot with voom transformed counts (log(CPM)). I didn't include all cases since there are too many. This look more or less what it is like in the one including every case. Target is the condition group in which I am interested. Differential expression analysis definitely showed that it is upregulated.

```
                     logFC       CI.L       CI.R    AveExpr         t      P.Value    adj.P.Val         B
GeneA  1.6576257  0.9143490  2.4009024  5.8713219  4.391796 1.646627e-05 0.0002581676  2.484413

```

Working steps TPM is basically from kallisto

log(CPM) calculation

```
#Import data from tximport
txi = tximport(files = read.hd5s,
                     type = "kallisto",
                     tx2gene = tx2gene,
                     countsFromAbundance = "lengthScaledTPM",
                     ignoreTxVersion = TRUE)
raw_data = DGEList(txi$counts)
# Filter low count genes
keep <- filterByExpr(raw_data)
raw_data <- raw_data[keep,,keep.lib.sizes=FALSE]
# Normalise counts
raw_data = calcNormFactors(raw_data)
# log(CPM) by voom
voom_data = voom(raw_data, design = NULL, normalize="none", plot = TRUE)
# Batch effect correction
voom_expr_psva = psva(as.matrix(voom_data$E),Batch)

```

Then I plot graph with TPM and log(CPM) respectively

# Answer
> 3 votes
For comparing the counts of different samples from DESeq2, Michael Love recommends using the variance-stabilized transform. It'd be great if you could provide some specific code examples in your question, but without that here's something that should work with the DESeq2 workflow as mentioned in the package documentation:

```
## [assuming "dds <- DESeq(dds)" has been run]
dds.counts <- assay(vst(dds, blind=FALSE));
## [This makes the scaled counts resemble more closely the actual read counts]
dds.quantile99 <- quantile(dds.counts[dds.counts > min(dds.counts)], 0.99);
dds.counts <- ((dds.counts - min(dds.counts)) / (dds.quantile99-min(dds.counts))) * (dds.quantile99);

```

The matrix *dds.counts* should now represent a model-normalised approximation of the log of counts, with column names as the sample names, and row names as the gene names. The default behaviour of the VST changes zero counts to an arbitrary positive number, which is why I've modified the output so that count values of zero will be changed to 0 in log space (i.e. a count of 1). To plot results for a single gene, subset based on the row name:

```
gene.VST <- sort(dds.counts["geneOfInterest",]);
plot(gene.VST);

```

# Answer
> 2 votes
To add to what @gringer, when you do use TPM, the normalization done is for both library size and gene length. When you use rlog, the normalization is done via median normalization (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2864565/). See below for a quick example about the rlog, and size factor normalization.

```
dds=makeExampleDESeqDataSet(betaSD=0.5,sizeFactors =runif(12))
rld = rlog(dds)
## you can see the data is variance stablized, i.e lower counts are converted
### see the size factor calculated
rld$sizeFactor
###
boxplot(assay(rld))
## look at unnormalized distributions
boxplot(log2(counts(dds)+1))
## now we normalize and see
dds = estimateSizeFactors(dds)
# the size factors should be similar
table(sizeFactors(dds) == rld$sizeFactor)
### here you can see the medians are similar but its a bit different from rlog
boxplot(log2(counts(dds,normalize=TRUE)+1))
###

```

To answer your question specifically, the sizeFactor is unaware of any batch effect. It assumes median expression of all genes to be similar and normalize using that. If you think there is a batch effect, it has to be present across many genes, as opposed to one gene. I would check it with a PCA plot (https://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/plotPCA) and then test for differential expression by including the batch effect in the model.

---
Tags: rna-seq, limma, batch-effects, tpm
---


# Question
Title: DMPfold: number of models produced using "10 500" feature in "run_dmpfold.sh"?
I ran this command:

```
run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./folder 10 500

```

From my previous discussion, I understand that I should run with default "run\_dmpfold.sh" parameters i.e. `3 50`; however, since I have already run the `10 500` I have some questions.

I see that after generating "sequence\_1.pdb to sequence\_500.pdb" files, these files are no more in the folder and, now I see "ensemble\_1.pdb". Will there be more ensemble files like "ensemble\_500.pdb" and "sequence\_n.pdb" files? I mean when can I expect "final\_1.pdb' file? Also, I see files like "CLUSTERS\_001.pdb".

# Answer
> 0 votes
The `sequence_n.pdb` files are temporary and are combined into an `ensemble.i.pdb` file after each iteration `i`. That file has many structures in. After the run is complete there will be a `final_1.pdb` file. Don't worry about the cluster files. If you are interested in what they do, consult the `run_dmpfold.sh` script and see what it is doing line by line.

As I have said multiple times elsewhere, you should run the script with default parameters and let it run *to completion*. Then use the `final_1.pdb` and any other `final_*.pdb` files as the output. If you want more structures, use the `ensemble.i.pdb` files.

---
Tags: protein-structure
---


# Question
Title: Is there a standard way to clean a PDB file and re-number its residues?
Is there a pre-existing tool which will tidy up the numbering of a PDB file?

Firstly, I would like to re-number the residues on inserts to make the icode an actual residue in the chain (by that I mean it's own number, shifting everything after it, so 45A becomes 46).

Secondly, I would also like to re-number so that each residue is one after another, as there are skips in the numbering, in spite of the two residues being physically next to each other.

I could do this in-house with a python script, but I thought I would ask specifically if there is a library or pre-existing tool?

Worst case scenario I have found this Perl script that does similar.

Thanks.

# Answer
Several years ago we have written a tool for this: http://kodomo.fbb.msu.ru/~avkitex/projects/pdbparser.html https://github.com/avkitex/pdbparserTools

The description is only in Russian but you can contact the author.
> 2 votes

# Answer
The pdb-tools suite has `pdb_delinsertion` and `pdb_reres` for this task. Something like:

```
pdb_delinsertion 1CTF.pdb | pdb_reres > out.pdb

```

should work.
> 1 votes

---
Tags: python, protein-structure, pdb
---


# Question
Title: 1000 genomes missing SNPs on chr X in GRCh38
I am trying to get a reference together for imputation using the 1000 genomes data set. The GRCh38 version of the data set was made available earlier this year, but seems to be missing a significant number of SNPs on the X chromosome (~96%).

Is there any resource that can explain why? I have tried searching on their home page with no luck.

The data sets in question are GRCh38, GRCh37.

# Answer
> 3 votes
As I told you on 1000 genomes helpdesk:
Because of the ploidy settings we used for our variant calling, we only called variants in the the pseudoautosomal region on the X chromosome. We are working on settings to correctly analyse X variants in both males and females and plan to analyse X properly soon.

# Answer
> 2 votes
I think these are biallelic snps:

> The set is restricted to biallelic SNVs and INDELs.

If you check the liftover version, the chromosome X is complete...

---
Tags: snp, impute2, 1000genomes
---


# Question
Title: Determining the quality of a publicly available genome assembly
I am looking for homologous genes in various protozoan species using BLAST. The genome sequences of these species are deposited in NCBI's WGS database. NCBI’s webpage includes the global statistics for each WGS project and it indicates the assembly level (scaffold, contig, etc.), the assembly statistics (number of contigs, N50, L50, etc) and (if available) the genome coverage. Can any / some of these parameters be used to evaluate the genome assembly’s "quality", especially if there is more than one genome assembly for the same species? Can these data be used to determine if the genome is very fragmented or not?

I understand that quality is a relative term and that there are several options for assessing genome assemblies (QUAST, BUSCO, etc.) but I am just looking for a "simple" way to compare several genome assemblies (for e.g., this post How can I improve a long-read assembly with a repetitive genome? states that the reference genome is highly fragmented, but there is no explanation regarding how the OP determined this). I believe that the N50 value can be helpful, but I need some context to understand what is "good" or "bad". I am at a loss of where to start, so any information will be greatly appreciated!

-Leah

# Answer
It is not so easy to determine the quality of an assembly - but sequence length, number of contigs, N50 and L50 can give you some guidance.

How much does the sequence length deviate from the expected length (or how much variation is among assemblies of the same species). Same for number of contigs/scaffolds. The sequence length, and number of contigs together with N50 and L50 give you an idea about the fragmentation. Very small N50 compared to expected genome size and chromosomes is a likely indication of high fragmentation, same for high L50 compared to expected number of chromosomes.

However, these number can be misleading. High N50 and low L50 can occur witho some very long contigs, so you may want to compare N75/L75 or better the Quast statistics
> 3 votes

# Answer
Carambakaracho's answer is good. A high N50 should usually be a good sign, as well as agreement between the overall genome assembly length and the expected genome size. The number of 'N's in the assembly may also be important to check, as Ns are sometimes used to indicate single unknown bases, and sometimes to indicate unknown numbers of unknown bases, and Ns may be used to connect scaffolds,which may lead to artificially higher N50. If you have access to a set of genome sequence reads (Illumina) for the same organism, then another more detailed check you can make is to compare the kmer spectrum of the reads with the kmer spectrum of the assembly, using the kat tool.

Finally, if you want to compare several assemblies, you might align them against each other, to see how they compare qualitatively (e.g. using mauve).
> 1 votes

---
Tags: genome, assembly
---


# Question
Title: How to predict stop codons in Illumina reads?
I have Illumina MiSeq paired-end reads from 150bp amplicons mapped to my reference genome (\> 1000X coverage). These reads have indels that may or may not induce frameshifts. If the indel induces a frameshift (i.e. indel is not a multiple of 3bp), I want to know if that leads to a stop codon later in the read. Is there any package/method that would allow me to do this?

Few extra details

* I might be wrong but I think this is different than gene prediction as I would want the tool to scan each read in the correct frame; not try different frames to detect ORFs etc.
* I expect mosaicism, so I do not want to compute a single assembly/consensus and predict stop codons in it.

Here is how I imagine such a script/tool would work:

1. Take the position of the first translated codon of the exon I am interested in (from Ensembl for instance?). *Eg. gene abc/exon1 starts at chr1:1000.*

Then for each read:

2. Look at which position the read is mapped to my reference. *Eg. read 1 is mapped to chr1:1050.*
3. Compute the position of the first full codon in the correct frame in my read. *Eg. start codon at chr1:1000 but my read starts at chr1:1050, the first full translated codon in my read should start at chr1:1052 (I think?).*
4. Split all my read into 3-mers starting at that position and build something like a 3-mer/position dictionnary.
5. Look for stop codons and report the positions.

Does it exist?

# Answer
I just recently discovered the ampliCan R package from the Valen lab after asking this question. I highly recommend it. It does all of what I was expecting & more.

To install: https://bioconductor.org/packages/release/bioc/html/amplican.html

To run (& other explanations): https://rdrr.io/bioc/amplican/f/vignettes/amplicanOverview.Rmd

This is the original paper & Do not forget to credit the people! https://genome.cshlp.org/content/early/2019/03/08/gr.244293.118.abstract
> 0 votes

# Answer
With 150bp reads, you might be able to blastx all the unique sequences you've got against the target proteins.
> 1 votes

---
Tags: sequence-alignment, sequence-annotation, illumina
---


# Question
Title: Find most similar genes to a set of genes of interest in time series RNA-seq data
I have gene expression data (RNA-seq) for 30 different time points (from 0 to 60 min each 2 min).

I have a set of 8 genes that behave similarly (although not identically) and I want to find the top X genes that behave most similarly to this set of genes in terms of their time-series expression profile.

I have tried so far to do unsupervised clustering (mfuzz) and check for genes that lie in the same cluster as my genes of interest. However, this approach relies on specifying a predefined number of clusters and since these genes are not exactly equal, sometimes they don't fall in the same cluster, depending on initial parameters.

Because of this, I have thought about doing it the other way around and look for genes similar to my set of genes. However, I haven't found any method/resource for doing that so far.

**So, is there a way to find genes with a time-series profile similar to a predefined set of genes?**

# Answer
> 1 votes
You can use GSVA to summarize the information of your 8 genes and look for the time-serie change using limma or DESeq2.

If you have enough samples for each time point you could also do a weighted gene co-expression network analysis (WGCNA), and see if those clusters in each time point are co-expressed similarly (But it can be tricky and you will need to define a lower bound of the clusters to 8 or so). You could also divide your time series in two point at the beginning and at the end, and see if there are differences or not.

If you want to get a list of genes similar to your set of genes you can use GOSemSim or BioCor (Disclaimer I am a contributor to the first, and author of the second), which can calculate similarities between genes based on the gene ontology and on pathways respectively. For this approach you'll need to calculate the similarity between all the genes and then select those that are closer to your genes of interest.

I'm not aware of any database/website with time-series profiles of genes (I assume that if there would be any would be on humans or some simple model organism).

# Answer
> 1 votes
DPGP (McDowell et al, 2018) uses a non-parametric clustering method, so should find clusters without you needing to specify cluster size or number of clusters. It is more or less explicitly designed for clustering timeseries gene expression data to find clusters of co-expressed genes.

https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005896

SplineCluster (Heard et al. 2006) uses a different method, but for a similar purpose and in a similarly non-parametric way, and worked well when we used it a few years ago with microarray data.

---
Tags: rna-seq, clustering, gene-expression
---


# Question
Title: Changing the axis limits of ggplot objects
By EnhancedVolcano R package I have this plot with this code

```
    > head(res1)
         baseMean log2FoldChange     lfcSE        stat     pvalue      padj
SAMD11   237.7975     0.03584334 1.3244737  0.02706233 0.97841002 0.9971691
NOC2L   4751.0381    -0.62090622 0.6500115 -0.95522347 0.33946470 0.9239759
KLHL17   735.0849    -0.95126291 0.8417444 -1.13010892 0.25843033 0.8926461
PLEKHN1 1117.5598    -2.11321359 1.0411108 -2.02976816 0.04238011 0.6350143
PERM1    409.0061    -2.32444927 1.1255813 -2.06511000 0.03891259 0.6285916
HES4    1088.7549     0.28395931 0.9091259  0.31234320 0.75477971 0.9866881
>
p= EnhancedVolcano(res1,
    lab = rownames(res1),
    x = 'log2FoldChange',
    y = 'pvalue',
    xlim = c(-5, 8))

```

I have stored my plot in an object

```
> class(p)
[1] "gg"     "ggplot"

```

I want to limit the axis length in this plot in the other word shrining my plot more but I don't know how to do that. `For example y axis be from 0 to 10 instead of 20`

Any help please

# Answer
Axes for `ggplot`s can be controlled via `xlim()` and `ylim()` functions. You can try:

```
your_ggplot_object +
  ylim(-5,20) # assuming that you want the gene labels to be readabe, not clear from your post

```
> 4 votes

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


# Question
Title: How many models should be generated by rosetta for denovo protein structure prediction?
How many models should be generated by rosetta for de novo protein structure prediction of a 382 long sequence? In the tutorial, it says 50,000 models but practically it takes a lot of time (weeks) to generate this number. For instance, I am currently running rosetta on a 382 long sequence and in a week it has generated 1500 models; I am using 16 cpus with Openmpi.

# Answer
Unfortunately you do in general increase your chance of success by going to the 50k range of models or higher. Projects like Rosetta@home exist because even large research groups can't access as much compute power as they'd like. In this case the long time arises from the combinatorial explosion of possible fragment combinations that need to be explored.

Success is rather protein dependent though - sometimes the right structure falls out after only a few models, other times your fragment library does not cover the protein adequately and you could run Rosetta forever without getting a good structure. Have a look at your 1,500 models and see how similar they look.

382 residues is rather long for Rosetta. Another factor to consider is whether you are using predicted contacts as part of the folding process (e.g. in the PconsFold protocol), which will greatly improve the chances of getting a good model and reduce the number of models you need. In fact the PconsFold group switched to CONFOLD for PconsFold2, I think mainly because it was much faster, so you are not alone in finding Rosetta runs long.
> 2 votes

---
Tags: protein-structure, computational-biochemistry, rosetta
---


# Question
Title: Parsing specific features from Genbank by label?
I'm trying to parse a Genbank file to find a specific feature. I can pull it out if I know the feature type (e.g. repeat\_region). E.g. if I'm looking for this feature:

```
 repeat_region   5623..5756
                 /label=5' ITR
                 /note="5' ITR"

```

I know that I can find it using:

```
for feature in reference.features:
    if feature.type == "repeat_region":
        print(feature.location)

```

But I don't trust that it will always be a repeat\_region. Instead, I'd like to look for it by label (5' ITR). I can seem to find a way to parse that from the feature object. Any suggestions?

# Answer
> 1 votes
The `/label=5'ITR` is called a **qualifier**. You can look for those qualifiers:

```
for feature in reference.features:
    for k, v in feature.qualifiers.items():
        if k == "label" and v == ["5' ITR"]:  # or "5' ITR" in v
            print(feature.location)

```

---
Tags: python, genbank
---


# Question
Title: Biopython Genbank writer not splitting long lines
I am parsing a csv file of annotated sequences and using Biopython to generate Genbank files for each. I want to add annotations of the sequence features.

My output file shows features listed without the correct line breaks. Other software is then unable to parse the names of the features

```
misc_feature    58..60  /note="Oxidation"
misc_feature    100..102  /note="Oxidation"
misc_feature    241..243  /note="Oxidation"

```

It should be formatted like this

```
misc_feature    58..60
                /note="Oxidation"
misc_feature    100..102
                /note="Oxidation"
misc_feature    241..243
                /note="Oxidation"

```

My section of code adds a space between the location and the /note="" but the Biopython writer doesn't add line breaks.

```
for i in range(len(start_list)):
    feature = Feature()
# Strict Genbank format calls non-standard things misc_feature
    feature.key = "misc_feature"
    liab_type = "{}".format(liability).split(":")[1]
#need to add 1 to start list since csv counting is off
    feature.location = "{}..{}".format(int(start_list[i]) +1, end_list[i]) + "  " + '/note=\"{}\"'.format(liab_type)
    container.features.append(feature)
    i +=1

```

The Biopython source is toggling the wrapping of the feature qualifiers on `no_space_key` and `Bio.GenBank._BaseGenBankConsumer.remove_space_keys`

How do I enable wrapping of text to format output correctly?

```
    def __str__(self):
    """Return feature as a GenBank format string."""
    output = Record.INTERNAL_FEATURE_FORMAT % self.key
    output += _wrapped_genbank(self.location, Record.GB_FEATURE_INDENT,
                               split_char=',')
    for qualifier in self.qualifiers:
        output += " " * Record.GB_FEATURE_INDENT

        # determine whether we can wrap on spaces
        space_wrap = 1
        for no_space_key in \
                Bio.GenBank._BaseGenBankConsumer.remove_space_keys:
            if no_space_key in qualifier.key:
                space_wrap = 0

        output += _wrapped_genbank(qualifier.key + qualifier.value,
                                   Record.GB_FEATURE_INDENT, space_wrap)
    return output

```

# Answer
I read the Biopython Genbank record source code carefully and realized I wasn’t using the Qualifier class correctly.

After reading the source code again, I realized I was using the data structure incorrectly. I need to use the Qualifier class to add a key and feature. This is formatted correctly and solves the problem.

```
feature = Feature()
# Strict Genbank format calls non-standard things misc_feature
feature.key = "misc_feature"
#need to add 1 to start list since csv counting is off
feature.location = "{}..{}".format(int(start_list[i]) +1, end_list[i])
liability_qualifier = Qualifier()
liability_qualifier.key = '/note="'
liability_qualifier.value = '\"{}\"'.format(liability).split(":")[1]
feature.qualifiers.append(liability_qualifier)

```
> 1 votes

---
Tags: python, genbank
---


# Question
Title: How to convert VCF format file into Fasta(.fa) with the Reference data
I need a Software or Python Program for Converting VCF(.vcf) to FASTA(.fa) format with the help of reference (.fa)

# Answer
> 2 votes
I used FastaAlternateReferenceMaker previously. You will need to download GATK first. Next you prepare your reference genome (reference.fasta) and your vcf file (input.vcf). You will call the GATK function like this:

```
 java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -V input.vcf \

```

output.fasta will contain the new fasta with snps inserted at sites specific by the vcf file

---
Tags: fasta, sequence-analysis
---


# Question
Title: How to load and split my two individual datasets when integrating datasets in Seurat?
I want to do comparative analyses of my two individual single cell RNA-Seq datasets in Seurat. And I load and list my data as follows:

```
Control.data <- Read10X("~/Documents/...  Control_filtered_feature_bc_matrix")

Control <- CreateSeuratObject(counts  = Control.data, project = "Control_summary", min.cells = 3, min.features = 200)

Mut.data <- Read10X("~/Documents/... /Mut_filtered_feature_bc_matrix")

Mut <- CreateSeuratObject(counts  = Mut.data, project = "Mut_summary", min.cells = 3, min.features = 200)

Object<- list(Control,Mut)

names(Object) <- colnames(Object)

```

Then the protocol https://satijalab.org/seurat/v3.1/immune\_alignment.html said, I need to use `SplitObject` function, but what do I split by?

```
Object_list <- SplitObject(Object, split.by = "?")

```

# Answer
`Seurat` is expecting individual datasets to be normalized separately prior to data integration. In that respect lists of objects corresponding to different datasets are handy to manipulate each object/dataset individually.

In your case you do not need to do any splitting, you have an object corresponding to each of your datasets.

If you would have started with a single count matrix of all datasets then you would need to split your data based on batches/experimental condition/...
> 4 votes

---
Tags: r, scrnaseq
---


# Question
Title: A program "PSICOV" "stopped" from a few days on Ubuntu 18.04?
I am running dmpfold, "seq2maps.sh" and from a few days the program is at the "PSICOV" stage:

```
psicov -z 15 -o -d 0.03 tail-fiber-sequence.aln

```

"ps aux" gives me this:

```
daniel     26693  490  2.5 2546588 2514320 ?     Rl   سبت20 25116:58 psicov -z 15 -o -d 0.03 tail-fiber-sequence.aln

```

Any idea on what could be the matter? Yes, I am runinng "rosetta" too on my 32 CPUs, so all my CPUs are covered to max 100%; howver, RAM is 93Gb and only 32 GB has been used. Thanks

# Answer
In summary just be patient. psicov doesn't parrallelise so a tenfold increase in CPU will not hurry it along.

Also psicov is not RAM intensive given the machine you are describing, so you haven't hit a RAM bottleneck.

I don't know the ps output format you are using I'm use to seeing %CPU and %MEM from `ps -aux` on Susie and OSX. You should also try,

```
top -c -u username

```

or for OSX `top -u`

You should see the psicov happily chuntering away consuming ~100% of a single process. There is a UNIX/Linux guru specialist on BioSE who could almost certainly provide better diagnostics.
> 0 votes

---
Tags: protein-structure
---


# Question
Title: Rosetta ab initio run gives 0.00 RMSD for generated 2000 models? (run still in progress)
I am running rosetta ab initio prediction using the command line:

```
"mpiexec -np 16 /main/source/bin/AbinitioRelax.mpi.linuxgccrelease -database /main/database/ @options -mpi_tracer_to_file log > log1.log &"

```

"Options" file looks like this: https://ibb.co/V9VrQHV

RMSD vs Energy plot is: https://ibb.co/1mjzwDj

There is no native structure as this is an ab initio prediction. Could you please let me know about your comments on the RMSD being 0.00 for all 2000 models generated till now? Is this normal and will produce results after more runs?

# Answer
RMSD is a measure of similarity between two structures. The relevant RMSD in protein structure prediction is between a model and the native structure. If you don't have a native structure you cannot have an RMSD in this sense. Ignore the fact it is 0.0, there is possibly some error that means you are calculating the RMSD between a model and itself.

And that's the difficulty of protein structure prediction... you don't know whether you are right or not! The structures you have generated are your results, and you'll need to find ways to filter or use them that are not reliant on RMSD. You could use estimation of model accuracy methods or clustering, for example.
> 1 votes

---
Tags: protein-structure
---


# Question
Title: split 1 column input into 5 column bed file
In `f1` below I am trying to split `$1` based on each line and create a bed file. If the line is a snp then the `:` is spilt and the text is `$1` the last digit is `$2` -1 and `$3` and the letter to the left of the `>` is `$4` and the letter to the right is `$5`. (lines 1 and 2)

If the line has `dup` then the `:` is spilt and the text is `$1` the last digit is `$2` and `$3` and the `$4` is `-` and the letter to the right is `$5`. (line 4)

If the line has `del` with nothing else then the `:` is spilt and the text is `$1` the last digit is `$2` and `$3` and the letter to the left of the `>` is `$4` and `$5` is `-`. (line 5)

I am not sure how to format for line 3 (the complex indel)... maybe split `$1` on the `:` and the text is `$1` the last digit is `$2` and `$3` and the letters to the left of the `>` is `$4` and the letters to the right is `$5`. Thank you :)

**f1**

```
chr7:140453145A>T
chr7:140453136A>T
chr7:140453135_140453136delCAinsTT
chr20:31022287dupA
chr19:13054614delG

```

**desired**

```
chr7    140453144   140453145   A   T
chr7    140453135   140453136   A   T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -

```

**awk**

```
awk 'BEGIN{OFS="\t"} {sub(/[^0-9]+$/, "", $1)}' f1

awk 'BEGIN { FS = "[ -]"; OFS="\t" } NF==3 { print $1, $2 - 1, $2, $3 } NF==4 { print $1, $2, $3, $4 }' f1

```

# Answer
> 4 votes
The following Python script seems to do the job.

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

from __future__ import print_function
import re
import sys

for line in sys.stdin:
    if '_' in line:
        match = re.search(r'(\S+):(\d+)_(\d+)del([ACGT]+)ins([ACGT]+)', line)
        print(match.group(1), match.group(2), match.group(3), match.group(5), match.group(4), sep='\t')
    elif 'dup' in line or 'del' in line:
        match = re.search(r'(\S+):(\d+)(del|dup)([ACGT]+)', line)
        if match.group(3) == 'del':
            print(match.group(1), match.group(2), match.group(2), match.group(4), '-', sep='\t')
        else:
            print(match.group(1), match.group(2), match.group(2), '-', match.group(4), sep='\t')
    else:
        match = re.search(r'(\S+):(\d+)([ACGT]+)>([ACGT]+)', line)
        coord = int(match.group(2))
        print(match.group(1), coord - 1, coord, match.group(3), match.group(4), sep='\t')

```

Invoked like so on the command line.

```
[standage@lappy ~] $ ./transform < f1
chr7 140453144 140453145 A T
chr7 140453135 140453136 A T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -

```

# Answer
> 3 votes
The Python script in my first response is heavy on regex matching, which is pretty clunky in Python. I like Python much better than Perl overall, but a throwaway script like this will be clearer and more concise in Perl.

```
#!/usr/bin/env perl
use strict;

while(<STDIN>) {
    if (m/(\S+):(\d+)_(\d+)del([ACGT]+)ins([ACGT]+)/) {
        print("$1\t$2\t$3\t$5\t$4\n");
    }
    elsif (m/(\S+):(\d+)(del|dup)([ACGT]+)/) {
        if ($3 eq "del") {
            print("$1\t$2\t$2\t$4\t-\n");
        }
        else {
            print("$1\t$2\t$2\t-\t$4\n");
        }
    }
    elsif (m/(\S+):(\d+)([ACGT]+)>([ACGT]+)/) {
        my $coord2 = $2;
        my $coord1 = $2 - 1;
        print("$1\t$coord1\t$coord2\t$3\t$4\n")
    }
}

```

It's been several years since I wrote Perl on a regular basis, so there are probably ways to make the script even more clear and concise. (Although I have a love/hate relationship with the default/hidden variables, which can make a script more concise but also harder for a newcomer to understand.)

Invoked like so on the command line.

```
[standage@lappy ~] $ ./transform < f1
chr7 140453144 140453145 A T
chr7 140453135 140453136 A T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -

```

# Answer
> 3 votes
Quick Perl solution (following on from Daniel's answer saying that this would be "clearer and more concise in Perl"). I was a bit confused with your BED output because you've treated the start/end locations of a single base deletion and insertion different to a SNP. In my example here, the start and end locations for those cases are the same to distinguish from the two-base change for TT-\>CA \[this can be changed\]:

```
cat 10454.txt | \
perl -ne 'if(/^(.*?):([0-9]+)(_([0-9]+))?(del|dup)?([ACGT]+)((ins|>)([ACGT]+))?/){
   print(join("\t",
         $1,
         $2,
         $4?$4:($2),
         ($8 eq "ins")?$9:($5 eq "dup")?"-":$6,
         (($5 eq "dup") || ($8 eq "ins"))?$6:($9)?$9:"-")."\n")
}';
chr7    140453145       140453145       A       T
chr7    140453136       140453136       A       T
chr7    140453135       140453136       TT      CA
chr20   31022287        31022287        -       A
chr19   13054614        13054614        G       -

```

**Regular Expression**

This uses enclosing groups to describe the important bits of the sequence. The groups are numbered from 1 upwards, with each number corresponding to a left bracket (in order), capturing the entirety of the sub-expression until the paired right bracket. I'll describe the groups, which are referenced in the subsequent part of the script as `$<group_number>`:

1. Any number of characters, up to (but not including) the first ':'
2. Any number of the characters '0' to '9' \[i.e. numbers\]
3. Optionally, a "\_", followed by any number of the characters '0' to '9'
4. The numbers from Group **\[3\]**, i.e. excluding the "\_"
5. Optionally, a "dup" or a "del"
6. Any number of characters selected from "A", "C", "G", "T" (i.e. bases)
7. Optionally, an "ins" or "\>", followed by any number of bases
8. The "ins" or "\>" from Group **\[7\]**
9. The bases from Group **\[7\]**

**Print Code**

The `print` statement prints out a tab-separated list, finishing off with a line break. I'll describe each list element:

1. Group **\[1\]**
2. Group **\[2\]**
3. If Group **\[4\]** exists, then Group **\[4\]**, otherwise Group **\[2\]**
4. If Group **\[8\]** is "ins", then Group **\[9\]**; \[otherwise\]
If Group **\[5\]** is "dup", then "-", otherwise Group **\[6\]**
5. If Group **\[5\]** is "dup", or group **\[8\]** is "ins", then Group **\[6\]**\[otherwise\]
If Group **\[9\]** exists, then Group **\[9\]**, otherwise "-"

# Answer
> 2 votes
Here is an `R` solution, would probably be slower than the suggested `Python` solution though:

```
bed <- readLines("bed.txt")

res <- data.frame()

for(i in bed) {

  temp <- strsplit(i, ":") %>% unlist()

  # case N>N

  if(grepl(">", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer() - 1
    end <- begin + 1
    ref <- str_extract(temp[[2]], pattern = "[A-Z]")
    alt <- str_extract(temp[[2]], pattern = "[A-Z]$")

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

  # case del_ins

  if(grepl("ins", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer()
    end <- str_extract(temp[[2]], pattern = "_[:digit:]+")
    end <- substring(end, 2, nchar(end)) %>% as.integer()
    ref <- str_extract(temp[[2]], pattern = "[A-Z]+$")
    alt <- str_extract(temp[[2]], pattern = "[A-Z]+")

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

  # case dup

  if(grepl("dup", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer()
    end <- begin
    ref <- "-"
    alt <- str_extract(temp[[2]], pattern = "[A-Z]$")

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

  # case del only

  if(grepl("del", temp[[2]]) & !grepl("ins", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer()
    end <- begin
    ref <- str_extract(temp[[2]], pattern = "[A-Z]$")
    alt <- "-"

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

}

print(res)

  temp..1..     begin       end ref alt
1      chr7 140453144 140453145   A   T
2      chr7 140453135 140453136   A   T
3      chr7 140453135 140453136  TT  CA
4     chr20  31022287  31022287   -   A
5     chr19  13054614  13054614   G   -

```

---
Tags: awk
---


# Question
Title: How to map selected genes to Metabolic pathway Maps
I have a selected Arabidopsis Genome Initiative (AGI) list for RNA seq and proteomics data, how can I map them to metabolic pathway maps to vilualize(e.g. TCA cycle / FA / Photosynthesis) in KEGG or AraCyc or any freely available pathway maps?

Could someone share an example if possible?

# Answer
KEGG has a webtool for this, KEGG Mapper.

You can just copy-paste your gene ids in there and see how your pathways are colored.
> 1 votes

---
Tags: r, rna-seq, pathway, proteomics, kegg
---


# Question
Title: How I deal with this kind of gene expression comparision
I have RNA-seq from two sequencing batches; Lab technician says that he has run the RNA expression quantification two times in bathes 1 and 2 for example tumor 1 in batch 1 and tumor 1 in batch 2 , normal 2 in batch1 and normal 2 in batch 2. PCA shows however shows two experimental runs are too close to each other and can be merged

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

Now my question is can I compare only `a tumour sample with its own matched normal sample in one patient` to get differentially expressed genes for this `patient` by using experimental runs (technical replicates )? Design would be `Tumour1_batch1` `Tumour1_batch2` ***versus*** `Normal1_batch1` `Normal1_bacth2`

```
> head(mycols)
              Condition
Normal1_batch1  N
Normal1_batch2  N
Tumour1_batch1  T
Tumour1_batch2  T

```

# Answer
> 3 votes
No. If the only difference between Tumour\_batch1 and Tumour\_batch2 was the library prep/sequencing run, then all this will tell you is if the particular sample of RNA taken from the tumour is different from the particular sample of RNA taken from normal. It almost certainly is.

Indeed if you did technical replicates of two RNA samples, both taken from normal, you would almost certainly find differences.

What you want to find is the differences in RNA level between tumour in general and normal in general. To do that you need to make an estimate of the mean and variance of RNA levels in the tumour and in the normal. To do that you need multiple samples from each.

If you need differences within a single patient, it would be more honest just to compute the fold changes without any implication of statistical certainty. You would then need to think about alternate schemes for measuring the confidence in inferences drawn from the data when you consider the actual biological questions.

---
Tags: r, rna-seq, differential-expression, deseq2, design
---


# Question
Title: What does final generated models with TM-score of DMPfold signify in "sequence.ema.txt" file?
So, my runs have finished and got 5 final models. I predicted TM-score of one model to understand the result. The output of "predict-tmscore.sh" for "final\_1.pdb" is:

```
1 models read from PDB file
1 models read from PDB file
nres is 377
neff is 15423
dpscore is -300.599030
dpscorenorm is -0.082424
Probability of each TM-score bin (cols 1-10), estimated TM-score (col 11) and prob TM > 0.5 (col 12) written to sequence.ema.txt
Estimated TM-score of ../final_1.pdb is 0.5

```

could you please help me understand what does "nres", "neff", "dpscore", and "dpscorenorm" below mean:

```
    nres is 377
neff is 15423
dpscore is -300.599030
dpscorenorm is    -0.082424

```

and it also generates a "sequence.ema.txt" file with below value:

```
2.0308916646172293e-05 2.0682286958617624e-06 0.04564645141363144 0.14857041835784912 0.27177685499191284 0.31814685463905334 0.19150513410568237 0.022412151098251343 0.001872337656095624 4.743059980683029e-05 0.5036169363265117 0.5339840054512024

```

Please let me know what these above values signify?

And,I would really appreciate your comments on what does "rawdistpred.1" file mean? There are total 9 such files. When predicitng TM-score, do I need to use "rawdistpred.1" file for all final models or I can use any "rawdistpred.\*" file?

thanks!

# Answer
> 2 votes
This is described briefly in the paper. The four quantities are inputs to the estimator of model accuracy:

* `nres` is the sequence length.
* `neff` is an effective sequence count, or how deep the alignment was.
* `dpscore` is how well the final model matches the initial distance prediction.
* `dpscorenorm` is a normalised version of the above.

Don't worry too much about directly interpreting these. They are only given for completion. The important line is `Estimated TM-score of ../final_1.pdb is 0.5`, indicating that DMPfold just about thinks this is a decent model.

The "sequence.ema.txt" is more raw output giving the probability of the TM-score being in each bin. "rawdistpred.1" is the distance predictions from the first iteration. Always use "rawdistpred.1" as input to the TM-score predictor as this is how it was trained.

---
Tags: protein-structure
---


title: Convert FPKM to RPKM

text:
# Question
Title: Convert FPKM to RPKM
I have my 24 samples analyzed by Hisat2\>Stringtie, from paired end sequenced data. I have got FPKM data and am planing to convert them to RPKM. to be able to compare my samples with another already published.

1- How can I convert from FPKM to RPKM?

2- Is it an acceptable method to compare samples using RPKM between two different studies or I need more or different procedures?

# Answer
I don't know if you have followed the discussions from the last years, but FPKM and RPKM are obsolete. Statisticians found that FPKM and RPKM are not the right values that can be used for proper statistics. It is a wrong way of normalizing RNA-seq data. In stead they suggest to use raw read counts, which can be used in tools like edgeR or DESeq2. Also, CPM would be better than FPKM or RPKM, for example you can use (log2) CPM values in limma with trend.

But to answer your questions.

1. You can't convert, you need to recalculate from raw data.
2. No not acceptable. To compare studies you need to correct for batch, which you can do with statistical programs such as edgeR and DESeq2.

The best advice I can give is to start with raw data, also from the already published data that you would like to compare your own experiment with. And continue with edgeR or DESeq2 or the like.
> 4 votes

---
Tags: rna-seq
---


title: Protein model quality assessment

text:
# Question
Title: Protein model quality assessment
How do I perform protein model quality assessment?

I will obtain predicted structures from DMPfold, I-Tasser, and Rosetta, what should be the next approach to select a structure from among this set of models?

Note, I have been looking at some, appropriate methods such as MetaMQAP, QA-RecobineIt, but their servers are not working.

# Answer
> 2 votes
Have a read of the Estimation of model accuracy in CASP13 paper to get the latest rankings.

Generally you would run the model quality software on the models you have generated and choose the best. Also have a look at the structures in PyMol, you learn a lot that way. See, for instance, how much variety there is, and where the models vary. Is there a conserved core in many of the models? Or is there just a blob of "spaghetti" that doesn't look like an ordered protein?

If you have many models (e.g. from a Rosetta run) you might want to look into clustering the models first to reduce the number you have to work with.

# Answer
> 0 votes
Co-variance with convolution is an extremely hot method, so DMPFold via psicov. Checking the reality of the folding is obviously important, convolutions can be missed, inaccuracies can happen.

---
Tags: phylogenetics, phylogeny, protein-structure
---


# Question
Title: How to allow "half windows" in biopython?
I am using Biopython to examine protein sequences including the Kyte and Doolittle hydrophobicity. However, for short sequences, this returns an empty value; the value is not allowed if the window does not have 20 values.

```
from Bio.SeqUtils.ProtParam import ProteinAnalysis

window_length=20
edge=1
kyte = {'A': 1.8, 'R': -4.5, 'N': -3.5, 'D': -3.5, 'C': 2.5, 'Q': -3.5, 'E': -3.5, 'G': -0.4, 'H': -3.2, 'I': 4.5,
        'L': 3.8, 'K': -3.9, 'M': 1.9, 'F': 2.8, 'P': -1.6, 'S': -0.8, 'T': -0.7, 'W': -0.9, 'Y': -1.3, 'V': 4.2}

full_sequence_analysis = ProteinAnalysis("CYGTGAVYLLGLGIGGFS")
kyte_window = full_sequence_analysis.protein_scale(kyte, window_length, edge)

```

Is there are a built in Biopython option in this for "allow half windows/short windows" for incomplete windows? There is nothing I can see to suggest so in the docs.

# Answer
> 3 votes
You can choose between a `window_size` of 20 for sequences longer than 20 AAs, and the length of the sequence for those shorter::

```
seq = "CYGTGAVYLLGLGIGGFS"
window_length = min(len(seq), 20)

full_sequence_analysis = ProteinAnalysis(seq)
kyte_window = full_sequence_analysis.protein_scale(kyte, window_length, edge)

```

---
Tags: phylogenetics, proteins, python, sequence-analysis, phylogeny
---


# Question
Title: Metagenomic shotgun data with internal control
For 24 human stool samples I have metagenomic shotgun reads from an Illumina platform. The reads were filtered and mapped against a bacterial species library and specific maps to species were kept and counted. Importantly, before DNA extraction from any sample it was spiked with a solution of a fixed DNA concentration of *Thermus thermophilus*, a species not occurring in the human microbiota. The reads from *T. thermophilus* allow me to convert relative abundance estimates back to abundance in the stool sample, to absolute abundance.

I want to test whether *T. thermophilus* is indeed a good internal control that removes confounding by converting to absolute abundance. To test this, I compare between-species correlations on the relative abundance scale with their correlations on absolute abundance scale. I would expect less correlation in the latter. I first take the raw read counts of a random set of species, I divide counts by a sample-specific size factor estimate (total number of reads r in sample i / geometric mean of r's) to adjust for library size, then I do arcsinh transform (to stabilise variance) and then I look at between species correlation. Next I take again the raw read counts, correct by the size factor, then correct for an estimate of the coverage ( estimated proportion of *T. thermophilus* in sample i ), and I do the acsinh transform. Again I look at between-species correlations.

I compare species-pair correlation before and after *T. thermophilus* correction in a scatterplot (see attached) to visualise the effect of the correction. I expected that the point cloud indicates a slope \<1 reflecting that previously strong correlations have weakened on the absolute scale. But that is not what my result shows. Indeed, some weakly positive correlations have become stronger. How does this make sense? What other ways are there to test validity of an internal control?

# Answer
> 1 votes
The problem with correlations in relative scale is know on the literature, it has been suggested to change to other metrics. So I wouldn't use it for internal control.

As you can see on figure 3 of the linked article the correlation on the absolute scale tend to be higher. So I would say this is within expectations.

As other way for internal controls are blanks (Just water to see if there are any contaminants), and technical replications, to make sure that the measure is correct not some problem with picking some amplifications instead of others...

Also you don't mention it but the tool and library used to impute the microorganism is also important. You can check with other tools to see if the estimation is similar.

---
Tags: quality-control, metagenome, coverage, correlation, validation
---


# Question
Title: Best way to "fish" long reads that match a query sequence
Very simple set-up: We have a PacBio long high-fidelity (HiFi) reads, genome sequencing, and some of those that contain a particular sequence. I want to find out which ones. The length of the query sequence is 7kb, though ideally I would use as query a sub-sequence of only 1.5kb.

The output can the just the read ID or the full read sequence.

We tried DAMAPPER and it gives some decent results but it might not be appropriate to search for the shorter sub-sequence of 1.5kb. It also feels like an overkill for this purpose.

Bear in mind that this is my first foray into long-read sequencing so my mindset is still very much one of short, low error-rate reads.

# Answer
I have seen the blasr aligner typically used for aligning PacBio reads to a reference. Is there a reason why blasr would not be a good choice? I guess the very short sequence you are trying to align to might be an issue.

Minimap2 is capable of mapping pacbio reads, and can also do pairwise mapping to see if reads overlap with each other, so this function might help you: From https://github.com/lh3/minimap2

```
minimap2 -x ava-pb  reads.fq query_seq.fq > ovlp.paf

```
> 3 votes

---
Tags: sequence-alignment, pacbio
---


title: how to generate an scRNA counts file

text:
# Question
Title: how to generate an scRNA counts file
I have barcodes.tsv, genes.tsv and matrix.mtx file from my RNA seq data and I want to generate a scRNA counts files using R. When I tried using the `readMM` function, it gave me the following error:

```
Error: cannot allocate vector of size 185.1 Gb.

```

This happens even when I run it as a job. How can I generate a combined csv file using R?

Thank you.

# Answer
> 2 votes
I use the DropletUtils package from Bioconductor. Getting counts begins like this:

```
library(DropletUtils)
filePath <- "path_to_data/raw_feature_bc_matrix/"
singleCellExperiment <- read10xCounts(filePath, col.names = TRUE)

```

# Answer
> 0 votes
You probably need to load the matrices part by part in order to avoid loading a large size vector. One way is to use features like `nlines` or `skip` of the `scan()` function. Also, this is suggested code by 10xgenomics:

```
library(Matrix)
matrix_dir = "/opt/sample345/outs/filtered_feature_bc_matrix/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
features.path <- paste0(matrix_dir, "features.tsv.gz")
matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")
mat <- readMM(file = matrix.path)
feature.names = read.delim(features.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)
barcode.names = read.delim(barcode.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)
colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

```

---
Tags: scrnaseq, seurat
---


# Question
Title: Getting species from ensembl sequence id
I have downloaded some orthologues from ensembl website in fasta format. The problem is that it uses the ensembl sequence id as name for the sequences. I want to know the species name or any other information about the sequence (chromosome, location...) For example for this sequence id: ENSMODP00000018933. How can I do that? I have been looking at the REST API but the most I can find is an API call for getting the sequence passing the id as parameter. https://rest.ensembl.org/documentation/info/sequence\_id I have also tried this API call but does not work for sequence ids: https://rest.ensembl.org/documentation/info/xref\_id

# Answer
> 4 votes
The lookup/id endpoint will get it for you. Of you can just look up the three letter species code, in this case MOD, on the list in the documentation.

---
Tags: ensembl, api, orthologues
---


# Question
Title: Are there databases or resources that contain expression of proteins or genes before and after cancer treatment?
Are there databases that contain information about gene or protein expression before and after targeted therapy or chemotherapy?

I am curious about different cancer cell lines or cancer types.

Any suggestion on a set of resources I should look into?

# Answer
> 2 votes
One such database/tool is `“[PharmacoGx][1]: an R package for analysis of large pharmacogenomic datasets.” Bioinformatics (Oxford, England)`.

They brought together a large compendium of experimental data (cell lines and the like if I remember correctly from a talk), standardized these and made available within a special object class. I think there is genomic and transcriptomic data, not sure if there is proteomic data though.

---
Tags: public-databases, cancer, data-retrieval, public-dataset
---


# Question
Title: Can DMPfold and I-Tasser work well for de novo structure prediction of 1135 amino acid long sequence?
I am trying to perform de novo structure prediction of a 1135 long protein sequence; can DMPfold and I-Tasser work well for de novo structure prediction of 1135 amino acid long sequence?

# Answer
> 1 votes
For DMPfold you might like to look at Figure 5B of the paper, which implies that performance starts to drop off after ~600 residues. Only proteins up to 500 residues were used in the training set, for example.

A protein of that length is almost certainly multi-domain. I would consider doing domain prediction; we have a tool DomPred, but many others are available. Then run the sequences of the domains individually and recombine them later.

I don't know about I-Tasser's ability to do this, perhaps see if they mention it in their papers.

# Answer
> 0 votes
Generally, but not always the I-Tasser algorithm is not in the same league as a covariant mutation-based analysis (psicov) and don't think many people would disagree with this. I don't do I-Tasser but I'm aware it is the algorithm of choice for rough and ready predictions.

Psicov based models are essentially de novo predictions, based on amino acid mutation patterns, which I assume are refined against known homologous structure in DMFold (inaccuracies can creep through). I suspect the starting point of I-Tasser is a known homologue and it's associated known structure.

Having said that you need alot of sequence data to make psicov covariance work (\>8000 aligned sequences) and (strictly in my opinion) it will also depend on the speed of mutation. Higher eukaryote gene alignments could require much higher numbers of taxa in the alignment. The covariance value will give you the strength of the signal and values \>0.7 should be used. I personally consider\>0.5 to 0.7 borderline, which is fair enough because those are the generic thresholds in statistics. I don't know if you can access the underlying covariance values in DMFold, or moreover set the threshold.

So I-Tasser could outperform DMFold in given situations, but the underlying output of covariance based approaches should make it clear how well the approach has performed. Obviously if it performs well I'd go with that.

I would follow@jgreener 's advice on the follow up program.

---
Tags: protein-structure
---


# Question
Title: How to edit the headers of multiple fasta files from multiple folders?
My directories are organized as follows: one main directory, in which I have multiple directories that end with a number ranging from 314 to 727, followed by .3 . For example, 'mgm4761314.3'. Within these directories, I have fasta files named after a similar pattern (for example, mgm4761314.3.299.screen.passed.fna is in the directory mgm4761314.3).

There is only one such file (.fna) per directory.

What I've been trying to do is to write a loop that would go through all these folders and fasta files, and add the short number associated with the folder/file name (eg. 314, 315...727) to the headers. For example, for my first folder/fasta file, instead of:

```
>seq1
>seq2

```

I want:

```
>314_seq1
>314_seq2

```

The output would be exported to a file called, for instance, 314.fna.

I've sort of made it work for one number at a time, but I have plenty of files, so I'd appreciate if anyone could help me with the loop part! This is what I tried (and it didn't work). Using ubuntu.

```
for i in {314..727}; do sed "s/>/>${i}_/" [...]/*${i}.3/mgm4761${i}.3.299.screen.passed.fna > [...]/*${i}.3/${i}.fna; done

```

# Answer
You were so close! This should do the trick:

```
for i in {314..727}; do
    sed -i.bak "s/^>/>${i}_/" mgm4761${i}.3/*.fna
done

```

Note that I am using `-i` which makes `sed` edit the original file. Because I use `-i.bak`, this will create a backup file with the original data. Once you are sure the command does what you need, you can delete all of these files with `rm */*fna.bak`.
> 3 votes

---
Tags: fasta, bash
---


# Question
Title: How to map multiple CDS sequences to exome data using bwa mem?
Normally, I use bwa-mem to map my fastq reads for a bunch of species using one protein-coding gene CDS from my reference species. This was fine when I had to try it for just a group of genes...I would obtain single CDS fastq files and use them as a reference. However, we now want to obtain protein-coding sequences for ALL GENES in my exome datasets. If I have one file containing all the CDS from my reference species, can that be used in bwa-mem? The end of my custom script produces one multiple sequence alignment file when I ran it using one gene only...my concern is that bwa-mem will somehow replace certain CDS or would I just end up with over 20000+ MSAs? Which is fine because I need them to be individual MSAs for each protein-coding gene.

I use a simple python script at the end to extract base information from the mpileup files generated to create my MSAs.

# Answer
> If I have one file containing all the CDS from my reference species, can that be used in bwa-mem?

That's fine, normal usage of bwa-mem entails giving it a fasta file with all chromosomes/contigs/transcript/etc. in an organism, so your usage is no different. How your custom "make an MSA" script handles that is unknown of course.
> 1 votes

---
Tags: bwa
---


# Question
Title: bedtools intersect on very large .bam file - -sorted confusion
I need to do a bedtools intersect operation on a *very* large .bam file. When I use the standard bedtools intersect operation, the process consumes all the memory on my system.

From the bedtools documentation page, I find the `-sorted` option to `bedtools intersect`, which sounds like it is the perfect solution to the problem, but the documentation seems incomplete. From the documentation, I need to ensure that the .bed file is sorted and I can use either `bedtools sort` or the Linux `sort` command to accomplish this:

```
bedtools sort -chrThenSizeA -i input.bed > sorted.bed

sort -k1,1 -k2,2n input.bed > sorted.bed

```

If I then execute the intersect operation on the .bam file, I get an error.

```
bedtools intersect -sorted -c -a sorted.bed -b input.bam > intersect_linuxsorted.out
ERROR: chromomsome sort ordering for file input.bam is inconsistent with other files.

```

Tracking this down I find that `samtools sort`, which was executed on the .bam file, "sorts by the leftmost coordinate, or the read name if -n is used." It would appear there is no way to put the .bam file in order by reference name, which is what bedtools instructs.

For a workaround, `bedtools sort` does have the `-g` option, which allows me to specify the ordering of the .bed file. To get this ordering, I execute:

```
samtools view input.bam | cut -f3 | unique > order.txt
bedtools sort -g order.txt > alt_sorted.bed

```

This puts the .bed file into the same order as the .bam file, and the `bedtools intersect -sorted` command now works, but is it working correctly? It seems that `bedtools intersect -sorted` wants the files in consistent ordering, which can be done, but is that consistent with the intersect algorithm itself? From the documentation's instructions to sort by the first 2 fields, chromosome name and left-most coordinate, it seems inconsistent to order it as the .bam file is ordered, by read name and not reference name.

I'm not sure how to check this assumption.

EDIT: Some extra details regarding sorting of the .bam file.

If the .bam file is sorted with `samtools sort input.bam > sorted.bam` I get the error, 'ERROR: chromomsome sort ordering for file input.bam is inconsistent with other files\`

If the .bam file is sorted with `samtools sort -n input.bam > sorted.bam` I get this error, `Error: Sorted input specified, but the file sorted.bam has the following out of order`

It would seem that if I sort without -n, it thinks the chromosomes are out of order, and if I sort with -n option it thinks the coordinate order is wrong. The only fix I have is forcing the order of the .bed file with the -g option, but again, I have no idea if that is consistent with the underlying algorithm's assumptions for the -sorted option.

# Answer
> 4 votes
It sounds like bedtools is behaving properly. The bam file is sorted by read name only if you used the `-n` option in `samtools sort`. Otherwise, if you used `samtools sort` without the name-sort option, it is 'ordered' by reference name (which can be arbitrarily ordered, as I'll explain below) and then by position (which is numerically ordered).

I think you are getting confused about 'ordering' the reference names. There are multiple ways to sort names; Some people might want numeric sorting (chr1, chr2, chr3...) and others might want lexicographic sorting (chr1, chr10, chr11...). Because of this, many tools allow you to specify the way you want to order the reference names.

There isn't necessary only one correct way to sort the reference names; what matters is that the order is the same between the files you're trying to intersect. This is exactly the purpose of `bedtools intersect -g`.

If you want to check the order of reference names in your bam, you can easily do so by checking the header via `samtools view -h input.bam`. It will say something like this:

```
@HD     VN:1.3  SO:coordinate
@SQ     SN:chr1     LN:248956422
@SQ     SN:chr10    LN:133797422
@SQ     SN:chr11    LN:135086622
@SQ     SN:chr12    LN:133275309

```

The `@SQ` lines tell you how your bam is sorted. My suspicion is that this order is different from the Linux `sort` or `bedtools sort` default, which is why your `bedtools intersect` was failing.

# Answer
> 1 votes
As per Lynn's answer, sorting is dependent on how chromosomes/contigs are named, how you want them to be sorted, etc.

`samtools sort` or `sambamba sort` will sort as per the SAM specification by the order of the contigs in the header of the file. If you've aligned these sequences to a human reference genome, typically it will be `chr1`, `chr2`, ... `chr10`, `chr11`, ...

To sort the BED file in the same way, you can specify the sort order with the `bedtools sort -g` option and providing the explicit order.

Alternatively, use `sort` with the `-V, --version-sort` option, that provides a "natural sort of (version) numbers within text".

```
sambamba sort -o output.bam input.bam
sort -k1,1 -V -k2,2n input.bed > sorted.bed
bedtools intersect -sorted -c -a sorted.bed -b input.bam > intersect_linuxsorted.out

```

---
Tags: samtools, bedtools
---