text:
RNA-seq: How to get new expression count after normalization
I've RNA seq, Human, Paired-end data, Sample size is <40. These are aligned using STAR, RSEM processed. With RSEM I've TPM and expected counts, that is two files columns as individual IDs and row as gene names.
I'm interested to normalize gene data. With edgeR tutorial (link in the end) and few other online resources I see that after following steps there's an R object that contains norm.factors (Page 15) value for each individual.
I'm unable to wrap my head around it, now. If I'm interested to get normalized gene counts, can I go ahead and multiple each individual's norm.factor into its gene counts? For example, the expected count for IND1 for Gene1 is 100 and it's norm.factor is 0.80, can I say that the normalized gene count is 100*0.80=80?
I'm not interested to perform voomm, or differential expression analysis.
Follow up question: If I'm to sample few individuals' data from these RNA-seq should I
a) normalize all together and extract interested individuals, or,
b) extract interested individuals and then perform normalization.
I think since normalization uses per person library size, each person's normalization is independent of the others.  Can someone please correct me if I'm wrong?
https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf

If I'm interested to get normalized gene counts, can I go ahead and
  multiple each individual's norm.factor into its gene counts? For
  example, the expected count for IND1 for Gene1 is 100 and it's
  norm.factor is 0.80, can I say that the normalized gene count is
  100*0.80=80?

No; you divide, but if you poke around, you can probably find a way to get edgeR to return the normalized counts.

I think since normalization uses per person library size, each
  person's normalization is independent of the others.

Not true for the default normalization used in edgeR, TMM.  CPM normalization is independently determined for each sample.
In general, normalize everything together, unless you think that the algorithm assumptions might be violated by including all of them together.

I was thinking to get the quantile normalization as on the question link:  https://bioinformatics.stackexchange.com/questions/2586/how-to-apply-upperquartile-normalization-on-rsem-expected-counts
Thank you for your kind reply.  :)

I got what I was looking for  thank you @swbarnes  :) poke around thing helped :D


text:
PDB file downloading: pymol automation vs. manual
I automated a PDB download using a Pymol script (below)
python
pdb_lists = ['3LZG', '6HJN'] # lots more pdbs
for x in pdb_lists:
    cmd.fetch(x)
    cmd.select(x)
    cmd.save(x + '.pdb', x)
    cmd.delete(x)
cmd.quit()
python end

When the script is run within Pymol it pulls down a collection of pdbs, however this flopped when I  plugged it through a Biopython PDB parser, but was okay when I manually downloaded the same pdb file from RCSB.org
A manual downloaded PDB was perfect, meaning the amino acids which cannot be observed by X-ray crystalography (e.g. they not sufficiently stationary) were recognised as being absent from the parser (resulting in a protein ~390 amino acids), but the Pymol scripted download contained all amino acid residues, i.e. even those which couldn't be observed in the X-ray (resulting in a protein of ~440 amino acid residues). Any ideas why the script resulted in a pdb that flopped and how to correct it?

In PyMol 1.8+ fetch downloads mmcif files by default. Your script must convert mmCIF to pdb, that's why the result differs from a pdb file downloaded directly from wwPDB.

To expand on the comment by marcin:

fetch downloads files in mmCIF format by default (https://pymolwiki.org/index.php/Fetch). Not all PDB entries have PDB format files, e.g. due to too many chains. Presumably this is why the change was made, though mmCIF files tend to be larger and hence download slower.
When calling save with the .pdb file extension the structure is converted to PDB format. Perhaps the converter uses the SEQRES records to define the sequence, and hence gives bad files when there are some residues not resolved.

To correct it you could try:

Saving as .cif and reading in with the Biopython mmCIF parser.
Providing the type=pdb argument to fetch, which will force PDB format downloads. Some PDB files may not be available.
Don't do it in PyMol and instead use the Biopython structure downloader (http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc187).

Also, if download speed is limiting you then consider using the MMTF file format, which is a binary format and hence smaller. Biopython can read in and write out MMTF files.

Thanks point 2 was exactly how I solved the issue (and it worked). Let me look at the Biopython solution. I looked in detail at MMTF but it that appeared more aimed at pdb->MMTF.

So this was the code I used and it worked (as far as I could see) and refers to point 2 of @jgreener's advice and was based on the comment by @marcin
python
pdb_lists = ['3LZG', '6HJN'] # lots more pdbs
for x in pdb_lists:
    cmd.fetch(x, type='pdb')
    cmd.select(x)
    cmd.save(x + '.pdb', x)
    cmd.delete(x)
cmd.quit()
python end

There are some Bio SE members who have serious skills in Pymol scripting, and constructive critique would be very welcome.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
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

extract all the hsps from each alignment
re-align them using MUSCLE
generate a distance matrix from the msa file
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?

You can do this directly in blast

Anyway you have the gist but load the MSA into Mega as an easy way to build a distance tree. I'll forward details later

@MichaelG. Right, that would be very helpful thanks. I'll take a look myself meanwhile.

Hi @DavidYoung I've not managed to answer this today. Promise a comprehensive answer late tomorrow.

@MichaelG. Sure, thanks.

@MichaelG. I should add that, the NCBI tree view is what I'm trying to get. But as I am using the NcbiblastnCommandline wrapper with my own custom nr database, I am trying to reproduce the tree view myself.

Specific questions

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.

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 ...

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Can I use GSEA to look for enrichment in data that is grouped into sets but are not necessarily genes?
I have yeast growth data in various experimental conditions that are classified by the strain type. Can I use the GSEA tool from the Broad Institute to test for enrichment of specific groups of strains within the total population (i.e. gene sets)? If not, would something like a hypergeomtric test be useful for this kind of analysis?

Yes, if you have a ranked list, you can use GSEA to calculate the enrichment of the strains along the ranked list.
Note that hypergeometric test is for count data, so you would test if there are more strains of one type that fit on certain conditions.
I would recommend this article comparing different enrichment methods to understand well what they do.


text:
snakemake with --conda-prefix, can I use a pre-built conda environment?
I have read this question (Running Snakemake in one single conda env) but still have some doubt.
I am using the '--conda-prefix /some/dir' option and I have a rule in my snakemake file with the directive:
rule xxx:
  ...
  conda:
    "envs/some.yaml"
  ...

The first call of that rule will create the conda environment according to the yaml file in /some/dir/envs with some random name (in my case 'bf55a79e'). If I run the workflow in a different directory, the same conda environment (/some/dir/envs/bf55a79e) gets activated directly instead of being built again. So this is what I expected.
My question is:

is there some way of skipping the building and specifying that rule xxx should use an already existing conda environment somewhere?
does anyone know where the mapping between rule and directory gets stored?

If the environment already exists it will simply be used. You cannot, however, tell it to use an environment with a normal name (e.g., "my_env").
I assume you mean, "where the heck does the random looking conda environment name come from?!?", which is a very good question. The answer to that is that it's a hash of the --conda-prefix setting and the contents of the environment yaml file you specify. The actual code for doing this (I'm "stealing" it from snakePipes) is:

import hashlib
md5hash = hashlib.md5()
md5hash.update("what you gave to --conda-prefix".encode())
f = open("your environment yaml file", 'rb')
md5hash.update(f.read())
f.close()
h = md5hash.hexdigest()

h is then the directory name. There are actually two lengths of those, the shorter 8 character version that snakeMake uses by default or the full length version that the code above will produce. snakeMake will actually use either of these. If you're wondering, no I don't think any of this is documented anywhere, you have to dig through the snakeMake code to find it.

Ok, that answers the core of my question!
So if I change the way I refer to the conda-prefix (via symlinks or relative/absolute paths), the environment will get rebuilt. Which is perfectly acceptable as long as you understand what is happening.

Exactly, I'll suggest to Johannes that he document this somewhere :)

and just to be absolutely clear, if you don't specify a conda-prefix, the hash uses  the yaml file content only, right?

I actually don't know that off-hand, but I presume so. The conda environments end up under .snakemake in the output directory, so it's possible that it uses that path somewhere.

upon the first workflow run I had the conda directives in all rules, so that the virtual env was created (relative path to Snakefile is ../outSnakemake_ngs_bngs05b/.snakemake/conda/a87f68b4). In order to suppress the creation of another (identical) environment I commented all the conda directives and called snakemake -p -s Snakefile_v4_ngs_bngs05b --cluster "qsub -q researshq" -j 5 --use-conda --conda-prefix ../outSnakemake_ngs_bngs05b/.snakemake/conda/a87f68b4, though the jobs failed. Is this the expected behaviour i.e. can I suppress the conda directives by specifying the env?

If the --conda-prefix doesn't change then an identical env won't be made (it'd have the same hash anyway).


text:
How to take into account alternative bwa mem mapping when computing coverage
when mapping short reads with bwa mem, if a read has alternative mapping positions they are reported by bwa mem in the X0 and XA tags.
Now, let's say I want to compute the coverage of my bam file. How would I do for the alternative hits to be taken into accounts. What I mean is that, if a read QUERY has a selected best HIT position and 3 alternate HIT positions, I would like that QUERY read to be counted as 1 unit of coverage for the best HIT position, but also as 1 unit of coverage for each of the alternate HIT position.
My understanding is that samtools depth will only consider the best hit (which makes sense in 99.99% of the cases).
I guess I could edit the bam file with custom scripts so that the alternative hits are also reported as primary hits but I was wondering if there was an already existing tool (as I am not that confident with my ability to mess around with sam files).
Thanks a lot.

You would want to instead use the -a option with bwa mem to output all alignments. You would be best off filtering these with a custom script to only keep those with a close enough score to the primary hit. samtools depth can't be told to include secondary alignments (it's hard-coded), so you'll need to use something else like bamCoverage or samtools mpileup.


text:
Understanding Single Cell RNAseq Plots
I have just started some single cell RNAseq analysis with 500 cells. I have used Seurat for analysis. Following the guideline I have generated few plots. Now I am trying to understand how to read those plots. Can the community provide their input in the following plots. To give a brief background of the experiment. This data was processed in MinION. Its from human brain sample from our collaborator. In their shortseq data they identified 15 clusters. However, both the libraries were processed independently. The main question I am trying to answer is celltype specific isoforms. but in my dataset I can find only 6 Cell clusters. So, I am trying to understand why I am getting such a small cluster? Is it because I stared with only 500 Cells? Also I want to know what these two PCA plot is indicating. Thanks

What are your specific questions?

@haci am trying to understand whether the data looks good or not. This data is from long read RNAseq and I had never worked on even short read scRNA seq.

The data looking good or not is not a specific question, it is subjective. For example how many cell types or states do you expect and does it make sense (regarding your research question) to end up with 7 clusters in the UMAP plot? Unless you provide some info on your experimental set-up and research questions you cannot get useful feedback.

@haci Thanks for your input. I am editing my question

How many PCs did you use? And how many cells are there in the dataset where 15 clusters were identified?

Identifying 15 clusters from 500 cells is challenging. It will depend on the relative abundance of each cell type within the tissue. You can use the cell type frequencies found in the shortSeq data set to see what to expect from a 500 cell dataset using the howmanycells tool
About the pc plots:
The first is an elbowplot which explains how many variation is explained by each PC. The knee basically defines when components become less informative.
EDIT: I noticed that the y-axis of the elbowplot starts at 4 suggesting that the dimensionality of your data is most likely ~10 dimensions. But I expect many of the clusters that don't show up on your Umap plot to be present at very low frequencies.
The second plot is a jackstraw plot and provides a visualization of a statistical test that calculates the significance of your principal components. This one shows that all your components are significant. You can combine this with the elbowplot to estimate the dimensionality (10-15 dimensions) of your data and decide which components to use for downstream analysis. In general the more they dropoff to the dotted line the less significant they become.
While these plots are useful as a first impression, always rely on biological information and common sense to decide what information (here: components) to include in downstream analysis such as Umap, clustering, etc. For example: usually you don't want to include components mainly defined by ribosomal, mitochondrial or hsp genes.

Whats the importance of second plot? What can I read from this plot?

@user3377241 updated my comment

The first two plots are described on the Seurat PBMC tutorial. Specifically:

‘Elbow plot’: a ranking of principle components based on the
  percentage of variance explained by each one (ElbowPlot function).

Based on the tutorial, you should use at least 10 PCs.

The JackStrawPlot function provides a visualization tool for comparing
  the distribution of p-values for each PC with a uniform distribution
  (dashed line). ‘Significant’ PCs will show a strong enrichment of
  features with low p-values (solid curve above the dashed line).

Based on the tutorial, you should use more than 15 PCs.
No one can tell you the right number of clusters. If you already identified 15 clusters in a different dataset, then you should have cluster markers for those genes. You can check those genes in your dataset.

The first two plots are used in order to estimate the number of PCs to be used in later stages, for example clustering (this was split in to two functions in Seurat v3). The number of PCs selected would have an impact on the number of clusters obtained. The more PCs the more information for downstream applications. It would be a good idea to go over the genes within each PC, a trained eye can tell what kind of biology (or technical artefacts) contribute to the PCs.
Moreover, FindClusters() has a resolution parameter, which enables the number of communities (clusters) that is generated. In Seurat v3, there is a wide selection of different clustering algorithms but recommending one over the other would be beyond my expertise.
If you have a comparable dataset generated with some other technology/method, you can give a try with "integrating" the two datasets. An alternative non-Seurat approach that aims to do the same is MNNCorrect from the scran package.


text:
how to Know the longest ORF in Protein?
I have a protein Seq. and I want to know

how to Know the longest ORF in this Protein?

and if there is any tutorial for this process ??

Hi, you need to supply example data. The question has been asked and answered several times already and referring to past posts would help get you started.

If you have a protein sequence, you already have the ORF. I assume you mean you have a DNA sequence, not a protein sequence, right? In what species? Is there splicing? Is this a bacterium or a eukaryote? Please [edit] your question and give us some context. Add the actual sequence and the organism.

Proteins don't have ORFs, DNA sequences do.  What DNA sequences do you have?  If cDNA (i.e. transcript/mRNA) then an ORF finder will give multiple ORFs and the coordinates on the full-length input DNA which indicate it's length.

can you give me example of ORF finder.

EMBOSS has ORF finder, but it will find all ORFs, which may not be what you want (e.g. you may have to post-process to pick the longest of the overlapping ORFs).  Usually folks use a gene-finder instead, which includes an ORF finder but also other features, such as promoter/transcriptional start site and may include special handling for intron/exon boundaries (i.e. for eukaryotes).  I like Prodigal for microbiome data but there are several choices including genemark, fgenesh, glimmer/glimmerhmm, genscan although particular genefinders are usually suited for either prokaryotes or eukaryotes.


text:
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.

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.

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.


text:
DMPfold error while running protein structure prediction
I am using DMPfold to predict protein structure of a 382 residue protein sequence. All goes wel, and even the files are produced but during the "seq2maps.csh", I get this error in between:
Running PSIPRED & SOLVPRED
Running PSI-BLAST with sequence c600m2-tail-before-lysin-sequence.temp.fasta ...
Segmentation fault (core dumped)
[blastpgp] WARNING: Unable to open BLOSUM62
[blastpgp] WARNING: BlastScoreBlkMatFill returned non-zero status
[blastpgp] WARNING: SetUpBlastSearch failed.
BLASTP 2.2.17 [Aug-26-2007]

Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.

Reference for composition-based statistics:
Schaffer, Alejandro A., L. Aravind, Thomas L. Madden,
Sergei Shavirin, John L. Spouge, Yuri I. Wolf,
Eugene V. Koonin, and Stephen F. Altschul (2001),
"Improving the accuracy of PSI-BLAST protein database searches with
composition-based statistics and other refinements",  Nucleic Acids Res. 29:2994-3005.

Query= C600M2_00100 tail protein [Salmonella phage BPS17S6]
         (377 letters)

Database: c600m2-tail-before-lysin-sequence.a3m
           46,978 sequences; 11,323,713 total letters

SearchingFATAL: Error whilst running blastpgp - script terminated!
Running PSICOV
Running FreeContact
Running CCMpred
 _____ _____ _____               _
|     |     |     |___ ___ ___ _| |
|   --|   --| | | | . |  _| -_| . |
|_____|_____|_|_|_|  _|_| |___|___|
                  |_|

using CPU (10 thread(s))

Reweighted 30404 sequences with threshold 0.8 to Beff=22917.9 weight mean=0.753779, min=0.0232558, max=1

Will optimize 62686429 32-bit variables

iter    eval    f(x)        ║x║         ║g║         step
1       1       1.50338e+07 152357      3.021629e+11    1.55e-06
2       1       1.47369e+07 152356      2.1447226e+11   1.18e-06
3       1       1.44292e+07 152355      1.4646483e+11   1.22e-06
4       1       1.41268e+07 152353      1.0388224e+11   1.5e-06
5       1       1.38295e+07 152358      7.337667e+10    1.92e-06

Can you please let me know if I need to rerun seq2maps.csh because of this error? And, how this error can be solved? Thank you!

The blastpgp command threw a segfault, so either it has a bug or you ran out of memory. Yes, you need to rerun things, but do so on a different computer (with more memory) or using a different blast installation.

The comment by Devon Ryan is correct - looks like your blastpgp command is not working properly. For a protein of that length it's more likely to be a problem with the BLAST executable running on your system than a memory issue.
Try running the blastpgp command (https://github.com/psipred/DMPfold/blob/master/bin/runpsipredandsolvwithdb#L26) on its own and see what happens:
~/path/to/blastpgp -a 12 -b 0 -j 3 -h 0.001 -d yourfile.a3m -i yourfile.fasta

Thanks, I did this but got this erro:
blastpgp -a 12 -b 0 -j 3 -h 0.001 -d c600m2-tail-before-lysin-sequence.a3m -i c600m2-tail-before-lysin-sequence.temp.fasta
[blastpgp] WARNING: Unable to open BLOSUM62
[blastpgp] WARNING: BlastScoreBlkMatFill returned non-zero status
[blastpgp] WARNING: SetUpBlastSearch failed.
BLASTP 2.2.26 [Sep-21-2011]

Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.

Reference for compositional score matrix adjustment: Altschul, Stephen F.,
John C. Wootton, E. Michael Gertz, Richa Agarwala, Aleksandr Morgulis,
Alejandro A. Schaffer, and Yi-Kuo Yu (2005) "Protein database searches
using compositionally adjusted substitution matrices", FEBS J. 272:5101-5109.

Reference for composition-based statistics starting in round 2:
Schaffer, Alejandro A., L. Aravind, Thomas L. Madden,
Sergei Shavirin, John L. Spouge, Yuri I. Wolf,
Eugene V. Koonin, and Stephen F. Altschul (2001),
"Improving the accuracy of PSI-BLAST protein database searches with
composition-based statistics and other refinements",  Nucleic Acids Res. 29:2994-3005.

Query= C600M2_00100 tail protein [Salmonella phage BPS17S6]
         (377 letters)

Database: c600m2-tail-before-lysin-sequence.a3m
           46,978 sequences; 11,323,713 total letters

SearchingSegmentation fault (core dumped)

Could you please let me know how this could be solved? Thanks

Ok, I got it resolved by changing the correct blastpgp (legacy blast) executable. Thanks!


text:
Does this FASTQ data contain single or paired end calls?
I have this fastq data from GEO:
zcat SRR1658526.fastq.gz | head -n 20
@SRR1658526.1 HWI-ST398:296:C1MP4ACXX:1:1101:1093:2094 length=102
GATCTCTATTACTTTTTGAAGGATTNNNNNNNNNNAANTTTTGAATCANNNNNNNNNNNNNNNNNNNNNNNNNNNNNTNNNNNNNNNNNNNNNNNNNNNNNN
+SRR1658526.1 HWI-ST398:296:C1MP4ACXX:1:1101:1093:2094 length=102
@<@FFDFDHHFHHIIIIGHBGGGIG##########10#0:BGDDHHII######################################################
@SRR1658526.2 HWI-ST398:296:C1MP4ACXX:1:1101:1167:2107 length=102
ATAATATTGTAGATATAAATGTTATCTAATCTTATCTGATCAGCTTGCTNNATANNNNNNNNNNNNNNNNNNACNTATGNNNNNNNNNNNNNNNNNNNNNCC
+SRR1658526.2 HWI-ST398:296:C1MP4ACXX:1:1101:1167:2107 length=102
CCCFFFFFHHFHHJJIIIIJIJJHJHJJJJIJJIJJJJIJHIJJJJGII#####################################################

It is supposed to be paired-end sequences. Is the prefix @ and + are the R1 and R2? What's the convention here?

Entries in a fastq file occupy 4 lines each and R1 and R2 are typically in different files. Since that SRA project is 2x51 you seem to have run fastq-dump without the --split-3 option, so both R1 and R2 are merged together. Make your life easier and never use SRA, but instead just download the individual files from ENA.

what's ENA? and how what would be the right fastq-dump command?

ENA is the European Nucleotide Archive. Don't use fastq-dump, don't use SRA at all.

What would be the equivalent command?

wget or curl or just right clicking on the web interface and selecting save.

I see, so ENA is a European mirror of GEO?

Use aspera, just as with SRA.

the link will be here: https://www.ebi.ac.uk/ena/data/view/SRR1658526

The data shown in the fastq files are paired end data.
In fastQ files @ and + do not correspond to R1 and R2.
The wikipedia page for fastq files describe the following conventions:

A FASTQ file normally uses four lines per sequence.
Line 1 begins with a '@' character and is followed by a sequence identifier and an optional description (like a FASTA title line).
Line 2 is the raw sequence letters.
Line 3 begins with a '+' character and is optionally followed by the same sequence identifier (and any description) again.
Line 4 encodes the quality values for the sequence in Line 2, and must contain the same number of symbols as letters in the sequence.

When using fastq-dump
".1" and ".2" stand for read suffices assuming u used the -I flag.
In this case you can see that R1 (".1") and R2 (".2") are within the same file.


text:
MiXCR: only create a single export file for all clonotypes
I am using following command from MiXCR to both align, assemble and export my input files:
mixcr analyze amplicon \
    -s hsa \
    --starting-material rna \
    --5-end no-v-primers \
    --3-end c-primers \
    --adapters adapters-present \
    --export "-aaFeature CDR1 -nFeature CDR1 -aaFeature CDR2 \
              -nFeature CDR2 -aaFeature CDR3 -nFeature CDR3 \
              -aaFeature VGene -cGenes -count" \
    --assemble "-OassemblingFeatures=[CDR1,CDR2,CDR3]" \
    INPUT_FILE.fasta analysis_name

Unfortunately, this creates multiple output file (analysis_name.clonotypes.ALL.txt, analysis_name.clonotypes.IGH.txt, etc.), which stands in contrast to the output of the mixcr exportAlignments command.
How can I export a single txt file instead, while still using the mixcr analyze command?

This is from the MiXCR documentation:
Under the hood analyze amplicon is equivalent to the execution of the following
MiXCR actions options:
> mixcr align -s hs -p kAligner2 input_R1.fastq input_R2.fastq alignments.vdjca
... Building alignments
> mixcr assemble alignments.vdjca clones.clns
... Assembling clones
> mixcr exportClones --chains IGH clones.clns clones.txt
... Exporting clones to tab-delimited file

So when using analyze amplicon, you are using mixcr exportClones under the hood and not exportAlignments, hence the difference in output.
I believe what you are experiencing results from a bug, please see this issue. If this is the case, all these files should be the same according to the bug report. Fortunately the developers have already addressed the bug in version 3.0.3.

As I was already running version 3.0.9 of mixcr, the bug you found (which corresponds to my observation) should have been fixed. However, although running mixcr exportClones prefix.clna clones.txt returns only one file as expected, the same is not true for me when using the analyze amplicon command. So, even though those commands are supposed to be equivalent, their output suggests otherwise for me.

Currently in MiXCR SERIES 4 you can use --dont-split-files  parameter here.
Please check our new documentation portal:
https://docs.milaboratories.com/mixcr/reference/mixcr-export/#clonotype-tables
Also now we support presets that allow to perform analysis for a lot of standard protocols using a one line command.

After finding the correct entry in the MiXCR documentation (thanks to haci), I learned that the parameter --receptor-type is responsible for which sequences will be returned.
The default for this parameter when using mixcr analyze is --receptor-type xcr, and stands for 'all chains', thus returning such a vast output. Changing this parameter to e.g. BCR solved my issue.

Where can I find the default value for the --receptor-type parameter? I cannot find it in the manual, can you please point out the relevant section?

I found it in the analyze section of the MiXCR docs, under the optional parameters.


text:
blast nt vs Genbank
What is the difference between these two databases?
ftp://ftp.ncbi.nlm.nih.gov/blast/db/
 and
ftp://ftp.ncbi.nih.gov/genbank/
I understand that BLAST database has both nt for nucleotide and np for protein. But genbank is also for genomic database. So what is the difference between nt and Genbank databases?
Thanks

They are the same data, but the interegation of the database is very different. Blast uses a similarity algorithm but Genbank uses record referencing. It is better to limit the quesiton towards a specific biological problem for a more detailed answer

I think the BLAST folder contains pre-formatted databases for BLAST whereas the GenBank folder contains sequences.


text:
I have a list of protein targets. How to get structural data (atomic coordinate) from PDB?
I have a list of around 500 protein targets (human). I only have their protein symbols. As an example:
EGFR
AR
FOXA1
CXCR3
...
I want to get structural data (atomic coordinate) from PDB (Protein Data Bank) at once, either programmatically (Python or R) or with a web api. Do you know a clear step by step tutorial for this task or code example?
If that is not available, at least I need to find the PDB IDs of my list of protein symbols at once.
thanks

Biopython PDB is a good start together with associated documents. It is better to describe the list of proteins for a more detailed answer

What do you mean by a "protein symbol"?

I mean they are like gene symbols. For example, CXCR3, EGFR and so on.

I think I should have mentioned you @jgreener

Please [edit] your question and add a few lines of your input file with the list of proteins. That way, we can use that to test our answers and we can know what kind of IDs you have. Also, if all you have are things like CXCR3, so gene names and not protein IDs, we will need to know the species.

There is not a simple way and there are factor that you also need to consider before starting.

Are you using gene names or gene accession numbers/id (NCBI/Uniprot)?
If not, are the names the current standard ones (e.g. scraped names from old papers is a definite no)
PDB residue index is rarely the same as protein index.
Do you want a perfect match or do you want say 80% homology?
Which PDBs do you want for a protein?

Personal example
For a website I have developed I have a route that gets all the PDBs for a given gene. (enter link description here). The way I did it is that I parsed Uniprot using a custom parser and made a dictionary/DB for each species with many synonyms for a gene going to a Uniprot ID and then from there I have another dictionary/DB that goes from Uniprot to PDBs if any. I went for Uniprot as I needed other data too (e.g. features) to present to the user. Including a renumbering option.
All this is way overkill for what you need. However, if you specify the species I can easily share a JSON file online.
SIFTs
If point 3 is important and you have Uniprot IDs I would suggest going with SIFTs. It has the offsets and has mapping for Uniprot IDs (pdb_chain_uniprot.tsv).
Gene names
In terms of mapping gene names to NCBI/Uniprot/Ensembl IDs, there are mapping tables everywhere. But if you may have obscure synonyms, using the NCBI API is a good option.
Uniprot data is better for protein and has PDB codes and chains, but the API gives a XML reply which is a bit annoying (I wrote a parser for it, but I think there was going to be one in Biopython).
Homologues
If you would like homologues, say your organism is human and you are okay with mouse, you can either use the threaded SwissModel models or use the blast query in RCSB PDB (not the PDBe, which has better data normally) or NCBI set to PDB database.
Swissmodel actually contains also PDB structures, so that may be a good starting point if your organism is one of those where precompiled exists.
Parts
Some proteins, especially mammalian are formed by separate domains linked together, especially the docking scaffolds. In which case you may care for a specific model within the gene, hence the indexing issue and the potential getting features.
If you want a single PDB structure for each gene and do not care about the technical I have some Pymol-using py3 code to fuse models together with the N and C termini aligned on line. But I strongly discourage that and would argue for protein feature based routes.

Thanks for the suggestions. I get some more information but I think I could not get the specific answer. My main problem is this: I have my human protein targets (around 500) as gene names (like EGFR and so on). I also do not mind info from other species as I can input this species info into model as well. I want to develop ML model and predict compound hits for each protein if possible. For this, I need to have structural data for each protein. I want to be able to get it with an R or Python script of possible. If you know of such software or tutorial over a script, I would be glad to hear.

In which case. Download the Swiss-model collection and see if it suits. It runs off Uniprot IDs, so here is a JSON file, containing a dictionary with key=names and value=uniprot id. As mentioned the Swiss-model file contains both the renumbered structures and threaded structures (=close homologue structures converted to that species basically).

In terms of your application, I am not sure what you mean by "compound hits". You mean epistatic mutations (i.e. one residue change alleviates the strain from another mutation) or is it a typo of compound het? I am assuming the former as the latter makes less sense. To calculate the ΔΔG (the energy of being folding) of a mutation you need forcefield calculations (I'd recommend looking into Rosetta pmut_scan as that has double mutations too. Not sure how ML fits in, but that is your remit. Although it does sound like Tom Blundell's SDM.

I think the JSON file helps answering my question, thanks, and I will accept that. Regarding second comment, there may be a confusion. I am trying to find ligands that bind to my target proteins.

Ha. That's funny and awkward for me. I talk to folk doing Fragment based screening most days and they say compound hit all the time and I didn't click that it was that. Sorry.


text:
ATAC-seq macs2 peak splitting in sliding windows
This question has also been asked on Biostars
I used macs2 to call peaks for atac-seq data. now my goal is to split the peaks into 50 bp windows with 25 bp steps and then calculate the Tn5 integration frequency in each window.
How should I proceed with that?

import numpy as np
np.random.seed(0)
import pyranges as pr

gr = pr.random()
gr.Score = np.random.randint(100, size=len(gr))

gr = gr.slack(25) # make data wider for this example

print(gr)

t1 = gr.tile(50)

def increase_by_25(df):
    df = df.copy()
    df.Start += 25
    df.End += 25
    return df

t2 = t1.apply(increase_by_25)

tiled = pr.concat([t1, t2]).sort()

print(tiled)
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | Score     |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 5205300   | 5205350   | +            | 33        |
# | chr1         | 5205325   | 5205375   | +            | 33        |
# | chr1         | 5205350   | 5205400   | +            | 33        |
# | chr1         | 5205375   | 5205425   | +            | 33        |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 41326450  | 41326500  | -            | 3         |
# | chrY         | 41326475  | 41326525  | -            | 3         |
# | chrY         | 41326500  | 41326550  | -            | 3         |
# | chrY         | 41326525  | 41326575  | -            | 3         |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 7,964 rows and 5 columns from 24 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

This is as far as I can get without example data and a clearer explanation of what you need.


text:
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...

So what is the actual question here? How you can root the tree?

The question is how to make a robust amino acid phylogeny with adequate robustness support. This is a generic question that many people will want answering and something I know a lot about. In my opinion, this has been a successful investigation given that the final tree is not a "rough and ready tree", but a formal likelihood-based tree. I see this as the end of the investigation and final response

I just have not got around to answering at this moment. I hope this is ok and we'll know in future when a rolling investigation is in-play either we  section the questions to specific topics or find an alternative forum

So @leah for forum rules it might be better to delete this and ask again OR change the title and specifically what you need is state of the art tree building to finish and conclude your study. Amino phylogeny is complicated and not trivial so in my opinion is a valid question. The idea is that someone in years in the future can understand from this post how to address a specific question. Stack exchange aims to get 90% of traffic via Google questions. Thus the question format is off-topic.

Exactly, thanks @MichaelG.! I agree that you two have had a very interesting discussion, it's really a shame it happened here since this site is not a good place to store discussions. It would be really, really great if you could combine all of this into one comprehensive Q&A though since there is a lot of great info you've both come up with.

Ok.. @leah I will try and set up a Q&A, I can't remember how to do this, but I'll give it a whirl. Essentially the site has a chat feature.

Thank you for pointing out the problem with this thread. As @Michael G. suggested, I changed the title so that it is better suited to the format of the forum. Hopefully these posts can be moved to the Q&A format so that the information is available for other people later on.  As for the question itself, I agree with Michael G. that we are almost done, and I am very much looking forward to his input to round things off and get some final conclusions on this investigation!  Again, many thanks for all your help, I have certainly learned a lot through these discussions.

I cannot comment... which seems odd while I can write answers.
Could you share with me the data you are using for the tree? I am developing a software for phylogenetic tree generation. I'd like to compare the tree and see if it can help you.
You can also try yourself if you like: https://www.dnaservic.es/
I find that the root it finds makes usually far more sense than with other algorithms.

Thank you @juanjo75es. I will give your software a try and I will let you know how it goes. I hope the results will helpwith the software you are working on.

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.

The mod is asking for a Q&A, when I find out what this is, this will need to be done.

Thank your for the detailed instructions @Michael G. I will give it a try and let you know how it goes. And thank you for the heads up about the program taking a long time to complete the calculation... I was not expecting that!

Hi @Michael G. Just a quick note to let you know that I am still working on the phylogenetic tree, but since I have no experience compiling, it is taking me a little longer than I expected. I will post the results as soon as possible.

Hi @leah you can always post it as a separate question

That's a good idea @Michael G. Will do!

Well, it has taken me some time but I have finally managed to install RAxML and FigTree in Linux (I could not compile RAxML on OSX). I took my alignment file, saved it in Phylip format and I tried to run the first command line to build the robust maximum likelihood tree (raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup-p 10000 -s ~/Desktop/alignment_file.phy), but I got the following error: Error, the outgroup name "outgroup-p" you specified can not be found in the alignment, exiting ....  and I don't know how to fix it. Could you please help me out? Thank you very much!!

Can you put this as a separate question please?

Certainly. The question is posted in this thread: https://bioinformatics.stackexchange.com/questions/10534/how-to-define-an-outgroup-to-build-a-robust-amino-acid-tree. Thank you @Michael G.


text:
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

You need raw count data, and let limma normalize it (with e.g., voom). Read the manual, it is pretty well written.

@benn But I have the normalized data already, I want to get single expression values

Try to get raw counts.

@benn Ok , and then, I should normalize with voom, after which how should I get DE values for each sample separately?

Please go over the limma manual (https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf)  and edit your question with the potential problems you will encounter. Be careful to use the right section given the technology used (microarray or RNA-seq).

@haci Hey, I have added the microarray tag, I am sorry for the confusion

Okay I assumed it was RNAseq... Just read the manual...

@benn I am sorry if I have confused you, I tried to find this in the manual but I could not find any instance of comparing single sample vs a population

If you do statistics you compare groups, not separate samples. Separate samples do not have information about variance which is needed for statistics.

@benn But what I am doing is One Sample t Test. It may not have meaning as such, but I wanted to try it out with DE Values

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.

Thanks for the help!


text:
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!

Please [edit] your question and show us line 177 of /media/Deepak/NPRP10/softwares/dmpfold/DMPfold-master/run_dmpfold.sh. As you can see in the error message, the script is looking for files named c600m2-tail-before-lysin-sequence_[1-9]*.pdb and it is doing so on line 177. What is this script? Where did you get it? Could it have hard-coded input files? All of your errors are about missing files.

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!

Ok, I cant ask more questions on this thread as stated by the person below. So, I will be asking next questions in another thread with title "DMPfold: is it suitable to get 5000 models to compare with rosetta".

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

Thanks for posting your solution! But please ask a new question for the follow-up question. This site is not a forum, so we don't do discussions, only questions and answers. You might want to take the [tour] or visit our [help] to understand more. In the mean time, I have removed the question from this answer, but please feel free to ask it as a new question.


text:
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!

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.

This would be a very good alternative. Thanks a lot @terdon for the answer.


text:
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.

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


text:
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!

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.


text:
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.

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!

Just one question: "output files before clustering" and "ensemble files": are these also 500 + 500 in number?

Yes, after a run has finished the sequence_n.pdb files will be 500 in number and correspond to the 500 structures in ensemble.n.pdb. Other ensemble files contain structures from previous iterations.

I noticed a behavior of the program I have: I left "run_dmpfold.sh" to run for two days without touching the system using the command line I mentioned, at that time dmpfold had produced  "sequence_71.pdb" file, and after 2 days when I check, there were no more number pdb files than "sequence_71.pdb". When I logged into the system as it was on sleep mode, the number of "sequence_71.pdb" stared increasing to "sequence_72.pdb"..etc. Hence, actually 2 days were lost where dmpfold was not producing files. Can you please let me know why this might have happened?

It sounds like your machine went to sleep, which will have paused all running processes. This isn't a DMPfold issue, it's more an OS issue. Have a look at tools like https://mac-how-to.gadgethacks.com/how-to/take-control-your-macs-sleep-functions-with-these-commands-0168109/ to sort this out.

yes, I corrected it by "aborting suspend" in my system settings. Thanks

Another query to clarify, I see that after generating "sequence_n.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"? I mean when can I expect "final_1.pdb' file? Also, I see files like "CLUSTERS_001.pdb". thanks!

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.


text:
Is there a JSON-based genomic feature format?
There's a variety of formats for storing genomic features (e.g. Genbank, GFF, GTF, Bedfile) but all of the ones that I'm familiar with use either a custom format (Genbank) or a CSV with defined column conventions and space for additional data (GFF, GTF, Bedfile).
Is there a feature format that is JSON based and adheres to a defined schema?

It is a nice question, but why do you expect JSON schema for this type of data What would be its advantages over other file formats?

Easy to validate, easy to understand, enforcement of conventions. I find that with bed or gff files I'm constantly either going back to the definition to find out what the columns are or trying to figure them out based on their content. This is especially annoying with a multi-column tsv file that has missing values for columns (as can happen in the genePred format).
Some file types (e.g. GFF) allow for additional fields whose names are not enforced. So one file could have geneName and another could have gene_name.

Bear in mind that this is a very subjective issue. For example, I find JSON an absolute pain to work with and very hard to understand, while simple, tab-delimited text files are a pleasure to work with. Note that if a field is empty, you have F1\tF2\t\tF4, so empty fields are also easily dealt with.

@terdon Completely agree on the subjectivity of the issue. I'm not trying to convert anybody. Just curious if such a format exists.

Would you be interested in solutions that somehow convert gff to json? I'm sure there are tools or Bio[Python|Perl|R] modules for this. It would help if you could clarify exactly what genomic features you're after though. Since it's such a vague term, you kinda need vague formats for it :/.

@terdon Yes, I would be interested in those as well. Genomic features is indeed a broad term but there are a number of concepts that reused repeatedly across formats. contig, start, end, gene, cds, mrna, trna come immediately to mind but I'm sure there's plenty of others.

So you are only interested in 1-dimensional, linear features? No trans-splicing, secondary, tertiary or quaternary structure etc? Should be trivial to convert that to JSON if that's what you want, yes.

Yes, only 1d linear features. No doubt about it being trivial to convert, I'm just curious if there's a format / schema people already use.

JVCF is a JSON schema for describing genetic variants. BioJSON is a JSON schema for multiple sequence alignments. The output of mygene.info's gene annotation web service is JSON-formatted.
You might build one of your own schemas from any or all of these, as inspiration for structure and content.
Consider that most genomics formats are historically tabular, because they are not hierarchical, because they originated prior to the advent of JavaScript as a common programming language (and subsequent use of JSON to serialize data), and because a lot of stream-oriented UNIX tools read in data and process it in ways that are natural to tabular structure.
Some tabular results can be hierarchical, but use inefficiencies in field values to apply hierarchy. JSON is inefficient for storing tabular data, but it is absolutely wonderful for hierarchical data. This difference is perhaps worth considering.
One could think of storing a gene as a master object containing various haplotypes, exons, introns, features, associated variants, associated TF binding sites, gene ontologies, etc. Ordering or prioritization of such items would be possible with the use of sorted lists within a gene object. The format is as extensible as it is a matter of adding a property, list, or object to the schema.
To gain adoption or "buy-in", build tools which publish it to users, and tools that consume what users get in order to do things, along with conversion and validation utilities.

Thanks for the detailed answer! I think you hit the nail on the head with tabular vs. hierarchical data. Formats like GFF, Genbank and genePred store hierarchical data in a tabular format which makes it difficult to parse. Contigs, genes, mRNAs, CDSs, introns and exons have a very natural organization that may be well suited to a JSON-like structure.

BioJS has capacity to read and write GFF as streams, it might have something of interest such as this: https://github.com/biojs-io/biojs-io-gff

Oh neat! Looks like this implementation converts GFF to its own internal JSON representation.

I've written a Python library called PyReference which converts a GTF / GFF (RefSeq and Ensembl) into a gzipped JSON file.
You can just use this if you like but there's also a Python wrapper around the JSON, which allows you to write genomics code more naturally.
Python isn't known for being fast but the library function for reading a JSON file is highly optimised - the following takes less than 4 seconds on my laptop:
import numpy as np
import pyreference

reference = pyreference.Reference()

my_gene_ids = ["MSN", "GATA2", "ZEB1"]
for gene in reference[my_gene_ids]:
    average_length = np.mean([t.length for t in gene.transcripts])
    print("%s average length = %.2f" % (gene, average_length))
    print(gene.iv)
    for transcript in gene.transcripts:
        if transcript.is_coding:
            threep_utr = transcript.get_3putr_sequence()
            print("%s end of 3putr: %s" % (transcript.get_id(), threep_utr[-20:]))

This looks promising, thank you! It would be even better if there was an API reference linked to in the GitHub repo ;-)


text:
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

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))


text:
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?

I'm not an expert but hierarchical clustering could be an alternative start

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).

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.

I already have a tool for doing the clustering (mFuzz). What I am looking for is to take 8 genes and look for the top X most similar genes to those genes. But I will try this DPGP, seems pretty good for the data I have!


text:
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")

Answered on GitHub here: https://github.com/satijalab/seurat/issues/2106

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.

Hi haci, I'm using Seurat v3. And I do not precisely know what you mean. Can you send me the code to have tried? Thank you.


text:
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

Why do you want to collapse replicates? It's really good that you have replicates. Treat them as such.

She's talking about technical replicates.

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.

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                             fastq2_lane1_batch1 fastq2_lane1_batch2 fastq2_lane2_batch1   fastq2_lane2_batch2  fastq2_lane3_batch1 fastq2_lane3_batch2 fastq2_lane4_batch1 fastq2_lane4_batch2

@Fereshteh That's a different question, and one you can find answers to by searching this site or biostars. Please do not stretch people offering to help by piling on more questions.

Sorry @Devon I think this question relates to this post; Can I use technical replicates for gene differential expression analysis? For example I am only interested in tumour A and normal A , can I compare them by using batch1 and batch2 as my replications?

I think you know the answer to that is "no".

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.

Thank you @swbarnes2, now my question is can I compare only one tumour with its own matched normal sample to get differentially expressed genes for this pair by using experimental runs (technical replicates )? Design would be Tumour1_batch1 Tumour1_batch2 versus Normal1_batch1 Normal1_bacth2

@Fereshteh As you're well aware, new questions should be posted as such.


text:
Coverage required
I was came across a problem during an exercise in a book and I don't really know how to solve it. I feel like something's missing.
"coverage, c = $NL/G$ (N=number of reads, L=read length, G=genome length)
total expected gap length = $G(e^{-c})$
total number of gaps = $N(e^{-c})$
Question: You want to sequence a 4 Mb genome by the shotgun method, by assembling random fragments with read length 500. What coverage would you require, to expect no more than four gaps, assuming no complications arising from repetitive sequences or far-from equimolar base composition?"
The answer from the book is basically just: Solve $N(e^{-c})=4$ for $c$, but the question doesn't give what $N$ is?

What do you mean? Isn't N the number of reads? That's what you say in your question.

yeah, what I meant was that the question doesn't give me the value of N, without which I can't solve for c

Solving for N and solving for C are the same thing (you know L and G in the first equation after all). It also doesn't give you a nice answer because the actual value (11,926,375 reads, which is ~14.9X coverage) isn't really a nice pretty number.

According to the book you mentionned Introduction to Genomics SecondEdition by Arthur MLesk p.113
(a) What fraction of a genome could you expect to assemble from eightfold coverage?
(b) What total gap length would you expect in an assembly of a 2 Mb target genome size from eightfold coverage?
(c) How many gaps would you expect in an assembly of a 2 Mb target genome size from an eightfold coverage of fragments with a read length of 500?
(d) You want to sequence a 4 Mb genome by the shotgun method, by assembling random fragments with read length 500. What coverage would you require, to expect no more than four gaps, assuming no complications arising from repetitive sequences or far-from-equimolar base composition?

I'm agreed that whitout knowing $ N $ it's not obvious to answer this last question so perhaps they're a mistake in the statement of this exercise.
As Devon Ryan suggested, let's assume eightfold coverage means the coverage $ C $ = 8X
So we can calculate the number of reads:
$ N = \frac{GC}{L} $
$ N = \frac{4000000\times{8}}{500} = 64000 $ reads
Then we can calculate the required coverage to expect no more than four gaps:
$ N( e^{-C}) = 4 $
$ e^{C} = \frac{N}{4} $
$ C = log{\frac{N}{4}} $
$ C $ = 10X
(PS: I would like to thank @LordOfRacoons to help answer this question.)


text:
Module preservation analysis WGCNA
I have done a WGCNA (weighted gene co-expression network analysis) analysis for various brain disorders. So along with control, I have 3 different disorders such as BPD, MDD and SCZ. Now I have to find out what sort of modules are preserved between the control and disease condition.
I have 54 samples in control and 50 samples in MDD as one of the condition. The other conditions are similar, somewhere between 50-60 samples each.
I have cited a paper where they have done module preservation analysis. Now I have checked the WGCNA or website they do have given example of mouse male and female liver module preservation.
I am stuck at this do I need both the expression data together to do module preservation analysis? Or do I first have to run WCGNA on both the data sets and then look for?
For example if I want to look at what are the modules preserved in CTRL and BPD condition.
I have followed this tutorial section 12.5 Module preservation between female and male mice
where it reads both male and female data together which is bit confusing
maleData = read.csv(file.path(data.dir, "LiverMaleFromLiverFemale3600.csv"))

I have also followed this tutorial:
https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/ModulePreservation/Tutorials/
Any suggestion would be really helpful.
Weighted gene network coexpression analysis (WGCNA) reveals stage-specific modules (gene network) in neonatal left and right ventricles


text:
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)
    }

}

Welcome to the site! Could you please post the code on how are you doing this ? (This could clarify things about the order and if it is correctly done)

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).

Thanks @devon that I understand (R and loops) but I am not able to confirm if assessing FPKM from normalized read counts make sense.

It makes more sense then using straight FPKMs. Inputting normalized values helps to negate the otherwise intrinsic issues with FPKMs (namely, that they're heavily affected by outlier genes).

To expand on what @DevonRyan is saying, if you calculate FPKM neat (without sizefactor adjustment), you are normalizing by library size. If you use DESeq2 to normalize, it is median normalization. You have to choose one of these, depending on what makes sense. Then using the normalized counts, you divide by transcript length etc, same value across all samples to get the final normalized values. You don't need to loop.

Thanks
So then if I apply these two normalization approaches in this order- median normalization using DESeq2 and then on normalized read counts library size normalization using FPKM; it is okay to do so

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.

Thanks @haci
But this is not what I am looking for. I have now added more information to my question.


text:
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.

Recent 1D data is also about 95% accurate.

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.


text:
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

in addition to TPM, I think kallisto also produces count-like data as output, which would be better for downstream analysis (presumably the default when tximport is used in 'kallisto' mode). The count-level data gives a better idea of uncertainty based on shot noise (i.e. the likelihood of missing a gene because its expression levels are low).

I think there might be some misunderstanding. I actually used the count data as shown in this line: raw_data = DGEList(txi$counts), which is to generate a DGEList object from the count-like data you mentioned.

You're using the countsFromAbundance="lengthScaledTPM" argument on tximport, which looks suspiciously like a form of normalisation to me ("Whether to generate estimated counts using abundance estimates"). The import default for that, as used in the kallisto import example is countsFromAbundance="no".

Thanks for pointing that out. I didn't notice since it was something shown to me by a bioinformatician when I first started learning R. Yes it looks like the reads are scaled to the library size and average transcript length as show in this tutorial. I guess raw_data = calcNormFactors(raw_data) was intended to renormalise to library size after filtering due to keep.lib.sizes = FALSE. Still trying to work this around my head.

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);

Thanks! Sorry for not providing an example. I thought it was more about the practice of the community than a technical question so I omitted it. I use voom transformation followed by psva() in sva package but I assume a model-normalisation can also be done. The reason why I have this question is that using TPM shows a result closer to our expectation, while limma actually uses the voom transformation counts to detect differentially expressed genes. I am just not sure if I am misleading the reader if I use TPM just to show a better looking graph.

Specific examples are helpful in bioinformatics because there are a huge number of solutions available that are context-dependent. It'd be great if you could add this information to your question to improve the quality of the answers you get.

Hi! I have added the graphs that I get and the results. I am not sure if that's an example that would help answering the question, but that's basically the problem I am facing: will using TPM misleading/exaggerating the upregulation and trick people in interpreting the data wrong? And is it better to use voom since they are variance stabled and batch effect corrected so that reader can better interpret the data e.g. observe any trend in clinical data rather than randomness that might lead to the difference in gene A expression level within the group?

It looks like you're doing a double normalization (i.e. using kallisto TPM values in addition to calcNormFactors(raw_data)). I'm not familiar with voom, so don't know if voom is attempting additional correction despite the use of the normalize="none" option.

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.


text:
From ensembl id to protein structure
Forgive me if this question is too trivial.
I have a list of ~700 Ensembl ID. I need to extract the protein structures of each Ensembl id. How can I do this?
Is there any script(python, R) so that we can download the structure from the website(like RCSB)?

Welcome, parse the pdb and then run a Pymol script. Otherwise Biopython will have a fetch which will likely pull down the pdb. Be careful of the format.

Thanks @MichaelG. As you understand, I am very novice in this field, can you be more specific with the answer. Like any link or blog.

Hi I've a lot on right now and I don't do Ensembl (I just use Genbank) so I'm not sure how to parse Ensembl seq -> pdb codes. Emily from Ensembl is on the forum, and might help.

Again you have to be super careful of the pdb formatting. Finally you might  want to split this into 2 questions, i.e. I think this involves 1. parsing 2. downloading.

Ensembl to PDB
So the dataset you need is from SIFT. There is one with Ensembl IDs. Majorly this also includes what chain is your peptide and what is the offset. About the latter, I actually wrote a blog post about it recently.
PDB download
Using Python, the Biopython PDB module is very good. Although the PyMol module is more powerful. But as you want only coordinates (right?), you'd best use the former, which has a better documentation.


text:
Find indels between two short sequences
I have two sequences, say AAAGCTCGAGG and AAAGCGAGG. I need a convenient tool which shows me insertions or deletions between these, i.e. in this case something like
AAAGCTCGAGG
|||||  ||||
AAAGC  GAGG
2bp deletion starting at base pair 5

and also insertions such as
AAAGCTC GAGG
||||||| ||||
AAAGCTCTGAGG
1bp insertion at base pair 7

Does something like that exist? Blastn finds the alignment of course, but only outputs the number of gaps, not their position. I guess I could parse a pileup, but do not want to create a bam file from one of the sequences first, so samtools seems to be out. Any further ideas?

IMO you need to write custom code. I just done this (again) via pandas, albeit elegant its a lot of code. Its quite easy to do via Perl5.

Smith-Waterman is guaranteed to give you the best alignment given your choice of substitution matrix.  For a few such tiny sequences, it's what I would use.  Many implementations exist.

Perhaps you could use something like SearchIO to parse the output of BLAST, to get the desired alignment attributes: https://biopython.org/DIST/docs/api/Bio.SearchIO.BlastIO-module.html

blastn should be able to output this format of alignment result, using the --outfmt 0 option

Is there any particular reason why you don't want a BAM output format? If the sequences are sufficiently long, then the CIGAR field from the SAM output of minimap2 could be used:
minimap2 -a seq1.fa seq2.fa | samtools view | awk '{print $6}'


text:
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.

Please alert users to a cross-post to show us that you do not take volunteer effort for granted. https://www.biostars.org/p/398893/

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?