text:
Get & Annotation of isomiR (isoform of miRNA) profile from RNA-seq
I should get and analyse isomiR (isoform of miRNA) profile from RNA-seq with FASTA format. The best tool for this is "DeAnnIso". But here is a problem with this tool, it only allows you to upload FASTA format files  with a maximum size of 50MB. However, the size of my FASTA file is more than 1000MB. I tried to compress the file in different ways and uploading, but it doesn't work.
Does anybody know any way to solving this problem?
Or Does anybody know any other easy way to get isomiR profile from RNA-seq with FASTA format?
A sample of RNA-seq like this:
>RR42357.95441.1 I-D00105:274:H7MDXX:1:112:8124:865 length=51
TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACCTTGTAATCTCGTATGCC
>RR42357.95442.1 I-D00105:274:H7MDXX:1:112:8133:866 length=51
AAAGTGCTACTACTTTTGAGTCTTGGAATTCTCGGGTGCCAAGGAACTCCA
>RR42357.95443.1 I-D00105:274:H7MDXX:1:112:8144:8637 length=51

Thank you very very much!

Have you tried to compress your file into .tar.gz? What about the way to create a "tag" file like recommended here?

Have you tried cutting the sequence names down to integers instead of all the embedded information and then compressing them?

In fact it looks like the header lines are longer than the sequence lines so I would expect a significant reduction in the file size.

^ this is funny.. but its true. microRNA isn't huge in length. I'm not exactly sure how 50 bp of sequence can take up 100 MB of data.

The obvious answer is to break up the file size and upload each separately,
First I'd check the file size,
wc infile.fa

Output

XXXXX      YYYY     ZZZZZ infile.fa

The number you want is XXXX - the number of lines in the file.
You'd then decide how many files you wanted to break it up into ... in your case divide XXXXX by 3. Then .. for example
sed -n '1,5000000p' infile.fa > outfile.fa
sed -n '5000001,2000000p' infile.fa > outfile2.fa
sed -n '2000001,4000000p' infile.fa > outfile3.fa

I'm not certain this would be cool however. An miRNA expert could answer this question very quickly. I think is you should check your data for sequence duplication. Simply piping the fasta file through a simple distance matrix and screening for "0" values would confirm this.
Creating a frequency table of duplicate RNA must be a feature of Bioconductor, but in any case screening for duplication (technically duplicate records) is a formal SQL database command.
It can be done a Perl/Python but its a bit clunky for a data set this size. I'd do it by setting the sequence as the key within a hash, the inefficiency has been solved (but its a weird way to do it).
You would have a database for "headers" for a given sequence, which would take up space, but I would assess the potential level of duplication first.

The tool you want to use accepts compressed files, so the first thing you should do is compress your file. I created a 111M file of random, 51bp sequences, each with a random, 52 character header and a random 51 nt sequence:
$ ls -lh file.fa
-rw-r--r-- 1 terdon terdon 111M Feb  9 12:21 file.fa

$ head file.fa
>YZF8SGLVBPCKMY18GJ94XTRCS61BB4OZKUZ61D 1J6R7H3 EJQA5 length=51
ACGAAGGGCATTTGGCAACGAACCTATACTATATCACGAGCAATTAGGCAG
>SPR0INGWKMILAXZZPF73BKRV73ELL1ROZSHYO1E9JX6P5X6MUZ20 length=51
TTAAGATTCCGAGAGAGTCTGTGGTTGGGTCCCTCTCCATAATGCTTTCAT
>U CB82L2M3BWJ9ZZZG3V6DL0Y0FTW8BXIFTU0TK7SL VN5QE42J6 length=51
TGGCTGGGAGGGATAGTAGCGAGATTATCCAGTCCGGTGTTCACAAGCCGG
>TGH62IM2TJK2QJTIVWGZR IVKOW94X7D52KGM7DJP61KKF81KPYA length=51
GCCGCAATACCAAAAAGATGGTGTGATTAGGAAATAAATAGGATCCAACAA
>L5VS9AVEF3UTEO3341STET092YEZ43ESVH32W7GI6Q27OIR NTXF length=51
CTATAGACGACTACGAACGACTGCAGGTGTGCCAGTTCTCGGTGAAACGGA

If we now compress that file:
$ tar cvzf file.fa.tar.gz file.fa
file.fa
$ ls -lh file.fa.tar.gz
-rw-r--r-- 1 terdon terdon 56M Feb  9 12:22 file.tar.gz

You will likely get a significantly better compression ratio since your headers won't be random and will likely be very repetitive.
If that isn't enough, you should try simplifying your headers. For instance, just use an increasing number:
$ perl -pe 's/^>.*/">" . ++$i/e' file.fa > file.short.fa
$ ls -lh file.short.fa
-rw-r--r-- 1 terdon terdon 58M Feb  9 12:23 file.short.fa

See what a huge difference that one change made? And if we now compress the file:
$ tar cvzf file.short.tar.gz file.short.fa
file.short.fa
$ ls -lh file.short.tar.gz
-rw-r--r-- 1 terdon terdon 19M Feb  9 12:24 file.short.tar.gz

Well below the threshold!

If compression isn't enough, you should make sure you have no duplicate sequences in the file. Use the FastaToTbl and TblToFasta scripts I have posted here and:
FastaToTbl file.fa | sort | uniq | TblToFasta > file.uniq.fa

Finally, if none of these helps enough, split your file into multiple smaller ones:
split -l 100000 file.fa

That will create multiple files of 100000 lines (so 50000 sequences, assuming all your sequences are one line) each.


text:
What is Feature Barcoding technology?
I heard the term Feature Barcoding a lot so I find resources to understand it. I have read docs of 10xgenomics and I see that Feature Barcoding technology is a method for adding extra channels of information to cells by running single-cell gene expression in parallel with other assays, this technology also do something with t-SNE and something with matrices.
I can't imagine the whole picture because I don't know what are those adding extra channels of information to cells, single-cell gene expression and parallel with other assays mean. And what is the input and output of this technology? Is there any paper go along with this technology?

Welcome to the site.  Have you looked at the 10X genomics website?  Usually they post all the papers and materials.  How you done pubmed or google scholar search for papers?

I did but I only found this paper relates to the term "Feature Barcoding" when some people tweeted about https://www.biorxiv.org/content/10.1101/237693v1 , and this paper in Feature Barcoding Glossary https://www.nature.com/articles/nmeth.4380 . But I cannot find the whole match "Feature Barcoding" in that paper which make me unclear how it relates to each other.

"Feature Barcoding" refers to using antibodies linked to specific oligos. The antibodies can then bind to cell surface markers of interest and the oligos they're conjugated with turned into normal scRNA-seq reads. This effectively means you can get a FACs profile for each cell without needing to run it through a FACs sorter and individually label it.
The "feature barcoding" nomenclature comes from 10X I think, but there are a few similar protocols (maybe they're exactly the same in the end, I'm not sure about the details) such as CITE-seq and REAP-seq. You can see a talk from 10X on this here.

Thank you so much. Now it's more clear for me. May you suggest any paper come with it? I found some paper but don't know which one is better for understanding: https://www.nature.com/articles/nmeth.4380 or https://www.biorxiv.org/content/10.1101/237693v1

I've updated my answer, the nature methods paper you linked to at least is related (I mentioned it in my update). The cell hashing paper you mentioned is also very related.

Thanks. I think this is it. I didn't know there is a talk about it. Thank you again for your answer. I will dig into it.

source: Stack Exchange

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

sort: votes


text:
Installing DESeq2 in Ubuntu
I am trying to install DESeq2 in my Ubuntu with R version 3.5.1. According to the package repository in Bioconductor the version should be 3.5.
> R.version

platform       x86_64-pc-linux-gnu
arch           x86_64
os             linux-gnu
system         x86_64, linux-gnu
status
major          3
minor          5.1
year           2018
month          07
day            02
svn rev        74947
language       R
version.string R version 3.5.1 (2018-07-02)
nickname       Feather Spray

then I attempt to install DESeq2, following the instructions in Bioconductor:
> if (!requireNamespace("BiocManager", quietly = TRUE))
       install.packages("BiocManager")
   BiocManager::install("DESeq2", version = "3.8"

and at the end I get:
Warning messages:
1: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘RCurl’ had non-zero exit status
2: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘GenomeInfoDb’ had non-zero exit status
3: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘annotate’ had non-zero exit status
4: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘GenomicRanges’ had non-zero exit status
5: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘genefilter’ had non-zero exit status
6: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘geneplotter’ had non-zero exit status
7: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘SummarizedExperiment’ had non-zero exit status
8: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘DESeq2’ had non-zero exit status

I have visited $2*10^{9999}$ discussions on this issue and none of them helped me solving the problem. In addition, I didn't have any problem installing other packages like limma and edgeR, so I assume that this might be specific to DESeq2 or one of it's dependencies.
Thank you for the help!

I used to have an issue where I had to run R-studio with sudo or my OS wouldn't allow me to install packages.

@Pallie, thanks! I had similar issues with other packages but in this case I run rstudio with sudo and it didn't help.

First I installed this libraries:
sudo apt-get install libcurl4-openssl-dev libxml2-dev

I then installed RCurl:
BiocManager::install("RCurl")

and lastly:
BiocManager::install("DESeq2", version = "3.8")

Have you tried just typing,
install.packages("BiocManager")

Without all the proceeding if statement? If this doesn't install correctly nothing else will.  Personally I would be tempted install each of the named libraries independently. If that works okay, e.g. its not a problem with the mirrors then I'd contact the developer michaelisaiahlove at gmail.com

BiocManager is successfully installed. The problem is with DESeq2...

Thanks it worked! I will post the answer for others with the same problem

Cool, happy to hear! Systems admin is the bane of every bioinformatician.

I noticed that some packages sometimes failed to compile depending on the version of gcc present on the system. So far, most of my problems were solved using a minimal R 3.5 installation, gcc 4.9 and trying to install everything using Bioconductor.

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


text:
UniProt flat file. Interpreting format of SUBCELLULAR LOCATION. What does flag mean?
In creating a parser for the UniProt flat file uniprot_sprot.dat (directory) in the manual for the section SUBCELLULAR LOCATION is listed

The format of SUBCELLULAR LOCATION is:
CC   -!- SUBCELLULAR LOCATION:(( Molecule:)?( Location.)+)?(
  Note=Free_text( Flag)?.)?
Where:

Molecule: Isoform, chain or peptide name
Location = Subcellular_location( Flag)?(; Topology( Flag)?)?(; Orientation( Flag)?)?

Subcellular_location: SL-line of subcell.txt ID-record
Topology: SL-line of subcell.txt IT-record
Orientation: SL-line of subcell.txt IO-record

What does Flag mean?
My best guess at present is that these are the text surrounded by parenthesis, e.g.
 Note=Not associated with microtubules (PubMed:21558460).

but that does not help in understanding why it is named flag.
While I can use my parser to verify this is correct when this section is working, it would be nicer to have a better understanding of what the word flag means in this context.
I also checked the XML schema but no help there.

Per request.
More examples from file.
These where gathered using grep.
grep -A 3 -n 'CC   -!- SUBCELLULAR LOCATION'  /mnt/.../uniprot_sprot.dat

71:CC   -!- SUBCELLULAR LOCATION: Host membrane {ECO:0000305}; Single-pass
72-CC       membrane protein {ECO:0000305}.

2163:CC   -!- SUBCELLULAR LOCATION: Membrane {ECO:0000305}; Single-pass membrane
2164-CC       protein {ECO:0000305}.

6003:CC   -!- SUBCELLULAR LOCATION: Secreted {ECO:0000305}.

7039:CC   -!- SUBCELLULAR LOCATION: Cell membrane {ECO:0000305}; Lipid-anchor,
7040-CC       GPI-anchor {ECO:0000305}. Note=In microneme/rhoptry complexes.
7041-CC       {ECO:0000250}.

7409:CC   -!- SUBCELLULAR LOCATION: Secreted {ECO:0000305}.

8677:CC   -!- SUBCELLULAR LOCATION: Virion {ECO:0000269|PubMed:9499098}. Host
8678-CC       endoplasmic reticulum-Golgi intermediate compartment
8679-CC       {ECO:0000269|PubMed:15016891, ECO:0000269|PubMed:9499098}.

12871:CC   -!- SUBCELLULAR LOCATION: Nucleus {ECO:0000250|UniProtKB:P48349}.
12872-CC       Cytoplasm {ECO:0000250|UniProtKB:P48349}. Note=Translocates from
12873-CC       the cytosol to the nucleus when phosphorylated.
12874-CC       {ECO:0000250|UniProtKB:P48349}.

grep -A 6 -n -m 200 'CC   -!- SUBCELLULAR LOCATION.*(.*)'  /mnt/.../uniprot_sprot.dat

1560049:CC   -!- SUBCELLULAR LOCATION: Secreted. Note=Dense core vesicles (DCV).

2452012:CC   -!- SUBCELLULAR LOCATION: DNA-(apurinic or apyrimidinic site) lyase,
2452013-CC       mitochondrial: Mitochondrion. Note=Translocation from the
2452014-CC       cytoplasm to the mitochondria is mediated by ROS signaling and
2452015-CC       cleavage mediated by granzyme A. Tom20-dependent translocated
2452016-CC       mitochondrial APEX1 level is significantly increased after
2452017-CC       genotoxic stress (By similarity). The cleaved APEX2 is only
2452018-CC       detected in mitochondria. {ECO:0000250}.

9243083:CC   -!- SUBCELLULAR LOCATION: Isoform Soluble copper chaperone CopA(Z):
9243084-CC       Cytoplasm {ECO:0000250|UniProtKB:Q59385}.

38197671:CC   -!- SUBCELLULAR LOCATION: E(rns) glycoprotein: Host cell membrane
38197672-CC       {ECO:0000269|PubMed:28290554}; Peripheral membrane protein. Virion
38197673-CC       membrane {ECO:0000269|PubMed:28290554}; Peripheral membrane
38197674-CC       protein {ECO:0000305}. Note=The C-terminus membrane anchor of Erns
38197675-CC       represents an amphipathic helix embedded in plane into the
38197676-CC       membrane.

4924705:CC   -!- SUBCELLULAR LOCATION: Note=(Microbial infection) In the presence
4924706-CC       of M.tuberculosis EsxA-EsxB complex decreased amounts of B2M are
4924707-CC       found on the cell surface (PubMed:25356553).
4924708-CC       {ECO:0000269|PubMed:25356553}.

02/18/2019
For the DE line it is noted

A block of DE lines may further contain multiple Includes: and/or Contains: sections and a separate field Flags: to indicate whether the protein sequence is a precursor or a fragment:

At present this seems to be related to the answer. When I have the file correctly parsed, then I can check for a correspondence as verification.

Of interest: UniProt subcell.txt - Controlled vocabulary of subcellular locations and membrane topologies and orientations.

Of interest: UniProt Evidence - Most information in UniProtKB has one or several “evidence tags” which describe the source of the information, e.g. an experiment that has been published in the scientific literature, an orthologous protein, a record from another database, etc.

Thanks. And sorry, I forgot to delete my comment after I added the example. Have you looked at the Swissknife perl module? This is produced by Uniprot precisely to deal with these files. That should at least point you in the right direction.

Ah, maybe. Perhaps to specify that the precursor is in one subcellular location while the mature product is in another. I don't know. And I couldn't find an example of these "flags" in an actual CC section of the uniprot_sprot.dat I had lying around. But I admit I just looked at a few, I don't really know how to grep for it.

Please don't update the question with what you learned, but post an answer instead! That way, the question will be marked as answered and can help future users. And thank you, I am very familiar with grep, I've been using it every day of my life for almost a couple of decades now.  :)

Thank you for asking this question. "Flag" is actually an obsolete notation which predates the introduction of proper evidence attribution in UniProtKB, and which escaped us in our efforts to update documentation to use the correct terminology. The correct  term would now be "Evidence", and the next version of the user manual will have this corrected (release 2019_06 due on July 3rd, 2019).
Before the current way of evidencing annotations was introduced in UniProtKB, the following "non-experimental qualifiers" or "flags" were used: By similarity, Potential, Probable. Absence of such a flag would generally imply that the information is based on an experiment.
When in doubt, please don't hesitate to contact the UniProt helpdesk.


text:
Getting the licenses of all conda installed packages
I would like to get a list of the licenses of all packages I installed using Conda. pip has what I need in the form of pip-licenses, is there a similar thing for conda?

I proposed grep '"license":' ${CONDA_PREFIX}/conda-meta/*.json in a similar SO query.

To my knowledge there's no single listing of the licenses, but you can:
conda list | awk '{if(NR>3) printf("%s=%s", $1, $2)}' | xargs conda info | grep license

Allegedly conda info --license should work, but at least in conda 4.6.3 it doesn't show any license information. It unfortunate that there's no convenient way to get a tabular list of package:version:channel:license information from the command line.

It doesn't work for me, I get InvalidVersionSpec: Invalid version '0.1.0_license=1.1_r-mutex=1.0.0_tflow_180_select=3.0_tflow_select=2.1.0absl-py=0.2.2ad3=2.2.1alabaster=0.7.10anaconda-client=1.6.14anaconda-navigator=1.8.7anaconda-project=0.8.2anyqt=0.0.8appdirs=1.4.3argcomplete=1.9.4asn1crypto=0.24.0astor=0.6.2astroid=1.6.3astropy=3.0.2atomicwrites=1.1.5attrs=18.1.0awscli=1.15.57babel=2.5.3backcall=0.1.0backports=1.0backports.sh ... etc BTW, I can edit your response, isn't that strange? Conda version 4.6.3. Conda list does not seem to output any licenses.

Try adding a \n in the printf ("%s=%s\n"). Regarding the edit, it would become a suggestion and I'd get a notice for approval.

Thanx, this does produce a list of licences indeed! However, I don't see what they belong to. I also don't understand how your command works otherwise I would distill from it something that prints a packages name and a license only. What does the command do? There are no licences in the conda list output so where do they come from? How would I get the license for one package?

This one uses JQ for parsing and querying conda generated json output (https://stedolan.github.io/jq/)
conda list | awk '{if(NR>3) printf("%s=%s\n", $1, $2)}' | xargs -I{} conda search --info --json "{}" | jq --raw-output '.[][0] | "\(.name)\t\(.version)\t\(.license)"'

builds the list, asks conda search for all info on package as json, feeds json to JQ, which selects the first item and outputs its name, version and license as tsv.

I like this solution, although it does take a bit of time, since conda search queries aren't that fast. The only problem I ran into is if you have a later version of conda installed where the command line conda is actually a function, not an executable file, you'll need to changes ... xargs -I{} conda .. to ... xargs -I{} $CONDA_EXE ....

I can't comment, but here's an improvement for the top answer to list more information:
conda list | awk '{if(NR>3) printf("%s=%s\n", $1, $2)}' | xargs conda info | grep 'name\|version\|license'


text:
Error in as.vector(x) : no method for coercing this S4 class to a vector
I tried to run the following code in R studio. Everything worked fine, except at the last step [write.table(mdat, "recount_mdat.csv")] when I tried to export the 'mdat', I got the following error:
Error in as.vector(x) : no method for coercing this S4 class to a vector

My code:
 library('dplyr')
 library('recount')
 library('magrittr')
 library('ffpe')
 library('RSkittleBrewer')
 library('SummarizedExperiment')
 library('devtools')
 trop <- RSkittleBrewer::RSkittleBrewer('tropical')
 colon_proj <- c('SRP025982', 'SRP42161')
 if(any(!file.exists(file.path(colon_proj, 'rse_gene.Rdata')))) {
 sapply(colon_proj, download_study)
 }

 blood_proj <- c('SRP066834', 'SRP041736', 'SRP060416')
 if(any(!file.exists(file.path(blood_proj, 'rse_gene.Rdata')))) {
 sapply(blood_proj, download_study)
 }
 proj <- c(colon_proj,blood_proj)

 dat <- lapply(proj, function(x) {
 load(file.path(x, 'rse_gene.Rdata'))
 return(rse_gene)
 })
 proj

 sapply(dat, dim)

 metadata <- all_metadata('sra')

 write.table(metadata, "recount_metadata.csv")

 mdat <- do.call(cbind, dat)

 write.table(mdat, "recount_mdat.csv")

is(mdat)
## "RangedSummarizedExperiment" "SummarizedExperiment" "Vector" "Annotated"

What is the class of mdat? Did you try is(mdat)? What did you get?

@IIrs Hi I got this after this step
           [1] "RangedSummarizedExperiment" "SummarizedExperiment"
           "Vector"
   [4] "Annotated"

Great, see how fast you got your answer when you provided enough details, remember to include all the relevant details next time! Glad you found how to solve the problem

You can't cbind a bunch of obscure object types.
If you want merged count tables you should do this:
mdat <- do.call(cbind,lapply(dat,assay))

Where row.names are Ensembl gene IDs and col.names are the SRR accessions.
Then run your table writing command.
-
If you want the coordinates of your genes then do this to make a bed with the genomic locations:
library(rtracklayer)
export.bed(rowRanges(dat[[1]]),'gene_coords.bed')


text:
How can I get plant Genome size, total number of genes, EST, and chromosomes number
How can I get plant Genome size, total number of genes, EST and chromosome number from phytozom and plantDB?

Have you checked the URLs from the image you pasted? The genome size should be trivial to calculate if you download the files. Also, please add the source of the image you are showing.

Both resources do not offer a well documented "programmatical" access option - either REST or something equivalent, so I believe it is not possible to (trivially) do it with Python or R.

It looks like you've got to download the files. A user defined script can do all the stuff you want.

Ensembl plants has a REST API, you might be able to get all these details (except possibly EST numbers) there: http://ensemblgenomes.org/info/access/rest
For example this query gives you the Arabidopsis genome size: http://rest.ensemblgenomes.org/info/genomes/arabidopsis_thaliana?content-type=application/json


text:
Using preprocessing/alignment functions on the server
I am new to bash and the processes behind cluster computing in general and need some help with understanding some basics.
After looking all over the internet and this forum (+ askUbuntu) I found nothing that addresses this issue.
I have a series of raw RNA sequencing files: .bam and wish to utilize samtools (or others) to begin the data preprocessing steps defined in many workflows.
I have installed samtools and many other programs necessary for the workflow, however, once I am on the remote server in bash, and try and run a function over the files of interest, it gives me the error: -bash-samtools: command not found
When I try to install the package, say samtools, by running sudo apt install samtools, I get an error stating I do not have permission to do such.
What do I do? Must I bring up this issue to the server provider or is there a way around it?
Thank you for your time,

The short answer is to use conda. In the bioconda channel we have most of the tools used in bioinformatics, such as samtools. You do not need administrator permissions to install conda or packages with it, so your lack of sudo ability is not an issue.

I have installed samtools and many other programs necessary for the workflow, however, once I am on the remote server in bash, and try and run a function over the files of interest, it gives me the error

Installing packages on your local machine won't affect the remote server unless they share a file system, so this is expected.
As an aside, BAM files aren't really "raw", they're considered processed data (fastq/fast5 files are considered truly "raw" data in the sequencing world).

At a very wild guess
-bash-samtools

This looks shaky. For a start you need a space between the command and the application in all Linux/Unix/Mac OS X applications (its complicated to explain right now [see below] just accept it).

Wait you are using sudo and thats how you are installing on a local machine. sudo on a remote machine will never work, because you are not the systems admin. sudo on a local machine is not cool (see below)

Anyway I need to do some next-gen and have 6 MiSeq kits waiting for action so I installed it from here:
http://www.htslib.org/download/
Okay so here's my tcsh history
 4  21:44   cd sam* # like I was in the Desktop
 7  21:44   more README # I needed to know how they wanted the installation
 8  21:45   ./configure # okay there are setting the paths for me - saves me a job
11  21:46   make # like fairytale stuff for a load of C scripts
16  21:49   ls | grep samtools # I wanna know if its compiled and I see samtools* ... yep
17  21:49   ln -s ./samtools* ~/bin # I'm gonna pop a link in my bin so after I type
18  21:49   rehash # 'cause I still insist on tcsh rather than bash (don't do it in bash)

However, the following is a bit creepy because I get loads of perl scripts in my /usr/local/bin ...
make install

Anyway,
samtools

It took 4 mins on an OS X machine, using tsch, after I just quaffed 2 pints at my  pub (10pm London time). So probability is set against me ....
OUTPUT
Program: samtools (Tools for alignments in the SAM format)
Version: 0.1.18 (r982:295)

Usage:   samtools <command> [options]

Command: view        SAM<->BAM conversion
         sort        sort alignment file
         mpileup     multi-way pileup
         depth       compute the depth
         faidx       index/extract FASTA
         tview       text alignment viewer
         index       index alignment
         idxstats    BAM index stats (r595 or later)

etc ...
BTW seriously, fun over, avoid "sudo" if you can, for lots of reasons (which I'm not gonna explain). Its a last resort and usually means your paths are incorrectly set.
sudo will be refused on a remote machine, without question. You COULD ask systems admin to do this for you. HOWEVER, I'd try and install yourself without sudo and your bioinformatics will benefit.

Program: samtools (Tools for alignments in the SAM format)
Version: 1.9 (using htslib 1.9)

Usage:   samtools <command> [options]

Commands:
  -- Indexing
     dict           create a sequence dictionary file
     faidx          index/extract FASTA
     fqidx          index/extract FASTQ
     index          index alignment

Output following make install

0.1.18 is a truly ancient version of samtools and is unlikely to be what you compiled.

Yes you are right, currently 1.9 (above)


text:
How to quantify similarity of genomes and find differences in set of S aureus genomes?
I have around 500 annotated proteomes of different bacterial strains and would like to quantify their similarity (or difference). I found gt genomediff from genometools gives me some scores that I can use to generate nice clusters, but I am not sure whether that tools really works. The fasta-files that I use contain multiple sequences.
I ran some tests and it looks ok.
1_reference.fna:
>1
TAAGTTACT
>2
TAAGTTACA

2_eq_to_ref.fna:
>1
TAAGTTACT
>2
TAAGTTACA

3_tags_diff.fna:
>1asdfadf
TAAGTTACT
>2asdfasdfffa
TAAGTTACA

4_orde_diff.fna:
>2
TAAGTTACA
>1
TAAGTTACT

5_add_subse.fna:
>1
TAAGTTACT
>2
TAAGTTACA
>2
TTACA

6_point_mut.fna:
>1
AAAGTTACT
>2
TAAGTTACA

7_different.fna:
>1
TAAGTTACT
ATTACCTAA
>2
AAAAAAAAA

Then:
gt genomediff --indexname test *fna
7
1_reference.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
2_eq_to_ref.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
3_tags_diff.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
4_orde_diff.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
5_add_subse.fna 0.206969        0.206969        0.206969        0.206969        0.000000        0.212596        1.349206
6_point_mut.fna 0.199527        0.199527        0.199527        0.199527        0.212596        0.000000        0.569180
7_different.fna 0.794496        0.794496        0.794496        0.794496        1.349206        0.569180        0.000000

I am running the clustering on the difference matrix calculated with genomediff: "These distances are Jukes-Cantor corrected divergence between the pairs of genomes, that is, the number of mutations per base between them."
Currently, we are studying S aureus. We genomes are assembled genomes (three different methods). My guess is that the sequences from the plasmids are present. Furthermore, we do have drug-resistance measured in culture. So, we will be able to compare the genomes and the resistances.

Can  help here (I think), but need more info. Could you explain why you talk about "proteins" but present nucleotide data. What do the fasta sequences represent and what is your goal of analysis? Also which bacteria, pathogens?

Arguably, depending on what the goal of using said differences is for, it would be far more useful to infer a phylogeny of the bacterial strains you're analyzing and consider differences within the context of their evolutionary relationships.

Thats what I was thinking too. In pathogens its quite common to track microbial drug-resistance and that tends to be presence/absence stuff.

The talk about proteins was an error, but I also have annotated proteomes generated from the dna sequence files. From here I would do hierarchical clustering. What tool would you use to generate a phylogenetic tree? First of all I want to quantify the difference between the strains to see whether drug-resistant or other classes of strains cluster together.

Ah! right now we have a question. Its a while since I've done anything on microbial drug resistance so I'll need to collect my thoughts.

In a second step, I would like to evaluate why they are (are not) clustering.

... stuff of original post deleted.
On second thoughts what you might be doing is a template assembly of your genomes. It is a possible interpretation of the 1.2. fasta sequences above (i.e. 1 is template).
A microbial genomics professional would advise also performing a de novo assembly, particularly if you are interested in presence and absence. The reason is that if a gene is present in your query and absent in your template, you will miss it. Again it all depends on what bacteria are being assessed, some are more prone to "genetic islands" than others.
You need a collaborator beyond that.
For the onlookers here bacterial genetic behaviour is very different to eukaryotes and what they get up, switching DNA etc.., to would appear bizzare from a higher eukaryotic world.

You mentioned a comparative analysis of 500 strains: I have worked on this in MRSA. Note, the bacteria is important considering the approach you will adopt.
Anyway, you want a single aligned file to produce a phylogeny, Bayes or likelihood. This is a model of point mutations. Bootstrapping doesn't really help because of SNP differences between isolates.
There is a complex problem. Generic phylogeny can hit the buffers in my opinion because often:

low numbers of SNP differences between isolates.
The error across the genome is poorly defined
Multi-clonal infection ignored and your isolates not cloned.
The other issue is the QC of the genome, I've been fairly astonished by the variation (it affects the tree)

Put it all together and you can get a nice tree but the topology of a given MLST is in my opinion unlikely to be correct. A referee may not bother with this, but they might.
You then map the phenotype against the phylogeny (the tips of the tree) and look for clusters. HOWEVER, there are problems:

Drug-resistant e.g. mec cassette (methicillin-resistant ) are often on plasmids and these are lost during isolation. Albeit they can integrate onto the genome. So an isolate can be drug resistance but you fail to find the gene.
Drug-resistance transmission is unlikely to cluster, its too quick, so you don't see nice tight clusters against a theoretically perfect tree
The best approach is drug-resistance in culture.

Hierarchical clustering, has many meanings. It is used in presence/absence data, I demonstrated the method failed (on bacteria) based on cluster analysis. It is currently revised and used in unsupervised deep (or machine) learning prior a training method.  The question then is what is what are you modelling? Presence/absence (of genes), point mutations, epidemiological data, drug-resistance/non-drug resistance?
My assessment is that given the potential complexity of these biological scenarios and the quality of the question you should seek formal collaboration and this is before considering what bacteria are being assessed. Some of epidemic drug-resistent bacteria rip up this rule book.

I am running the clustering on the difference matrix calculated with genomediff.  "These distances are Jukes-Cantor corrected divergence
between the pairs of genomes, that is, the number of mutations per base between them." Currently, we are studying S aureus. We do have drug-resistance measured in culture. And we do have assembled genomes. My guess is that the sequences from the plasmids are present.

JC is a primitive algorithm, particularly using cluster analysis ... a definite no no. Acually its quite rubbish and if I reviewed the ms I'd politely return it. On other issue is that I was authoritively informed that plasmids for MRSA (i.e. staph, i.e. S aureus) could not be guarrenteed and can be lost in isolation. E.g. the mec cassette may be present, but could be lost. I'm not a microbial biologist so wouldn't have a clue about plasmid retention, but it was the clear advise I got from exactly the same question, i.e. from people who are supposed to know.

Anyway I think we're done here ... at least get a negative rep for this thread (at least so far :-O )  ... Oh I see "S aureus" is mentioned in the title, my bad, ... the good thing is that you can use phylogeny on this bacteria, i.e. the theoretical model (which is what I do) matches the biology of the bacterial genetics ...

Michael, do you have a publication for your SA project? It seems to be very relevant for us.

I do, but its not an area I work on and it didn't float my boat, so I prefer to remain anonymous. There are lots of publications on MRSA you can follow. I personally don't think the area is very well developed in terms of modeling brilliance, but  one of Peacock's studies stands out.


text:
Error in PSI-BLAST when running psipred locally
I am very new to bioinformatics and Linux and am trying to run psipred locally on Ubuntu 18.04 but cannot get PSI-BLAST to work. I have downloaded BLAST legacy version 2.2.26 and unpacked it, uniref90 database and psipred version 4.01
I installed psipred according to instructions:
cd psipred/src
make
make install

And this seems to work fine
I used the following code to create a database, which proceeded without errors:
makeblastdb -in uniref90.fasta -dbtype prot -out uniref90

In runpsipred I changed the directories:
# The name of the BLAST data bank
set dbname = /home/hilda/uniref90

# Where the NCBI programs have been installed
set ncbidir = /home/hilda/blast/bin

# Where the PSIPRED V4 programs have been installed
set execdir = ./bin

# Where the PSIPRED V4 data files have been installed
set datadir = ./data

I then try to run psipred:
tsch ./runpsipred example/example.fasta

And get the error
Running PSI-BLAST with sequence example/example.fasta ...
FATAL: Error whilst running blastpgp - script terminated!

I am not sure how to solve this as I do not understand what is going wrong. As I said I am new to Linux. I use the command tsch to run psipred (otherwise I get the error "command not found"), but not to install the database, could this be a problem? Or have I missed something in the setup process? Running psipred does create a copy of the example.fasta in the psipred folder named psitmp4286007f010.fasta, but that is all that seems to happen.
I would greatly appreciate any suggestions on how to make it work!
Edit:
I tried to run with BLAST+ (which seemed to work without tsch) and get the following error:
BLAST Database error: No alias or index file found for protein database [/home/hilda/uniref90] in search path [/home/hilda/psipred::]

I have also tried setting dbname to home/hilda/uniref90/uniref90 and got the same results. there is a uniref90.pal file in the uniref90 directory, which I understand is the alias file? Other forum answers seem to use pfilt to index (?) the database, but there is no pfilt in my psipred directory. Any ideas how to fix the problem?

Welcome to the site Hilda. I think that the supported version is now blast+ or something alike. Is there any special reason why you need the legacy 2.2.26? It will get easier to get support and trouble shot the problems if you install the supported version. In addition, it could be that the legacy version does not work with Ubuntu 18.04...

I use tcsh, but I would make a new attempt using bash.

I picked the BLAST version based on what seemed to be the most commonly used in forums and similar questions, and that I got the impression that psipred preferes legacy to +. Will definitely try downloading and running BLAST+ instead. Thank you!

I'm having the exact same problem. I'm using BLAST+ though and runpsipredplus script. Did you manage to make it run and how did you do it?

This error could be due to the fact that you are using legacy blast, have you tried using BLAST+ instead?
If you wish to use PSIPRED with BLAST+, then use the runpsipredplus script
here, rather than the normal runpsipred script.


text:
How to create Phylogenetic Trees from fasta files in Python or R?
I have around a hundred Fasta files (and will collect several thousand) with DNA sequences and +50x coverage. What is a recommended method to construct a phylogenetic tree? Solutions in Python or R are sought.
I found Phylo from Biopython only handles already calculated trees.

Would you be open to using a program designed to do this or do you need to roll your own in a programming language?

I would not look for a package for this, but instead build a small pipeline calling external tools with something like the following workflow:

Cluster the ~100 sequences with CD-HIT-EST/PSI-CD-HIT or many other options
Take all the sequences that form one individual cluster and build a multiple sequence alignment (MSA) with MAFFT/ClustalOmega or similar
Take the MSA and build a phylogenetic tree with a Maximum-Likelihood approach like IQ-TREE or similar
Visualize the tree file with Jalview or similar

Of course this is rather general and depending on exactly what you're doing you may want a different workflow and/or different tools. You should also explore the parameter space, do not assume the defaults are necessarily good choices

Yes, I agree with this and thats what I do. R is has its uses however in phylogeny

Jalview can also make the tree for you. Just load the alignment into jalview and tell it to build a tree from it.

@terdon not recommended though, I think jalview only implements NJ

Oh sure, I am not saying it's the best tool. I just know that ~10 years ago when I last used it, it was perfectly capable of producing a decent tree. But then I was using about a dozen sequences, all of which were homologous. So not a very hard tree to build.

@terdon 10 years ago is about right - pre-history and even Perl programmers existed then ... I think you critiqued me for 'samtools' so touche!

I agree with the first two steps. For the tree building I might consider RAxML. I might then load the finished tree into R with the ape/phangorn package for manipulation/visualisation.  For stand-along visualisation figtree is great.

The obvious single answer is R "ape". This will give you access to PhylML for tree building and Clustal/Muscle for alignment building. The paths to the binarys are important. There are several distance methods in there such as NJ and BIONJ. Its distance approaches however don't look mainstream, but I could be wrong.
There are functions within ape which are cool, the tree sorting is very cool and I need to read through this with much greater care. Personally I wouldn't perform a core phylogenetic analysis within R, because the standalones are sufficient and the analysis is intensive.
https://cran.r-project.org/web/packages/ape/ape.pdf

I agree with Chris Rands that a reasonable approach would be to call external tools.
However, if you really want to do the phylogeny from within Python, you could use the P4 package, which is a bit complicated to handle but gives you lots of options in the way to build MCMC-based bayesian phylogenies:
https://github.com/pgfoster/p4-phylogenetics
You would still need something else to align the sequences before.
To visualize the tree using python, you could use the ete toolkit, which is likely more powerful than what you can find in Biopython: http://etetoolkit.org/

(if I understand your situation correctly) https://www.rdocumentation.org/packages/seqinr/versions/3.6-1/topics/read.alignment shows how to use the function read.alignment which can take fasta msf etc. The docs provide the example' read.alignment(file = system.file("sequences/LTPs128_SSU_aligned_First_Two.fasta",  package = "seqinr"), format = "fasta", whole.header = TRUE) but you can use this code below (assumes those files are aligned) to go from reading the tree to getting the distances, producing the neighbor joining phylogenetic tree, and then plotting the tree.
library("Biostrings")
library("seqinr")
library("ape")
library(phylogram)
library("dendextend")

fasta.res <- read.alignment(file = "geneticAlignment.msf", format = "fasta")
fasta.res.dist.alignment = dist.alignment(msf.res, matrix = "identity")
fasta.res.dist.alignment.nj = nj(fasta.res.dist.alignment)
plot(fasta.res.dist.alignment.nj, main = "from fasta files")

Within the R/Bioconductor ecosystem you could use the DECIPHER package (https://bioconductor.org/packages/release/bioc/html/DECIPHER.html) to align your sequences and then continue with either ape or phangorn  (https://klausvigo.github.io/phangorn/) packages to compute a tree. You can then visualize your tree using for example the ggtree package (https://bioconductor.org/packages/release/bioc/html/ggtree.html).


text:
Targeted NGS, up to 99% of reads have been marked as duplicates
Currently I'm performing whole analysis (pipeline from *.fastq to *.vcf) of 41  samples (targeted NGS). I rely on GATK best practices, however with some modifications. I decided to use the following tools:
#mapping
bwa mem (with mem - alternate haplotypes)

#preprocessing
RevertSam (Picard)
MergeBamAlignment (Picard)
MarkDuplicates (Picard)
BaseRecalibrator (GATK)
ApplyBQSR (GATK)
AnalyzeCovariates (GATK)
SortSam (Picard)

#call, filter and annotate variants
HaplotypeCaller (GATK)
GenotypeGVCFs (GATK)
snpEff, SnpSift

Nevertheless, after marking the duplicates with MarkDuplicates about 99% of reads have been marked as duplicated, hence the coverage decreased from several hundreds to literally up to five. So I decided to omit this step.
Have you ever encountered a similar problem? What is your opinion about this?
Best Regards,
Adam
EDIT
The Agilent Haloplex custom target enrichment kit (HaloPlex Target Enrichment for Illumina Kit) has been used to prepare the libraries, sequencing done using Illumina HiSeq 2500. Targeted NGS of all exons (only) of 90 genes.
Here I provide how the final bam look like. This is the result without duplicate marking, otherwise there're only several reads per exon...

can you clarify what you mean by targeted NGS? What are you targeting/capturing and how were the sequencing libraries prepared?

I believe, using haloplex target capture, you expect reads to have identical begin and end coordinates, so imho in this case you should skip duplicate marking/removal.

Wouter is right, Haloplex is PCR based, consequently you expect high duplication. Do not deduplicate this data prior to analysis.

Thanks for help! By the way, what do you think about GATK as a tool for such analysis. For me it's nice but a bit too complicated.

Indeed after MarkDuplicates the coverage was up to 4. It was hard to distinguish between actual alternations and mismatched bases.

Most duplicate marking programs do it by looking at the start and end location of hits. In this case, if two reads have identical start and end locations for both reads, then they will be considered [PCR] duplicates of each other. When running through a duplicate filter, you will only get one read for each unique pair.
This is much more likely for targeted sequencing, because the length covered by the template is much smaller. It's also more likely that you actually do have lots of PCR duplicates for the same reason.

If you're doing target capture, especially of a small region, it is entirely possible that there are only five unique reads. Target capture is molecularly finicky and sometimes the entire region is not captured.


text:
Retrieve RNA sequencing data for human p53 colon cancer cell lines
I want to download the next generation RNA sequencing data in different cell lines. I want the data (alongwith normal) for p53 wild type, p53 knock-out and p53 null cell lines (human colon cancer). From where do I download this data? In fastq format if possible.
Currently, I am going through SRA, to look for the datasets, using their search engine. But the process is extremely tedious.
Also can someone please hint about how to download the data from SRA, in case that's the only way to get data. I'm getting very confused.

Welcome to the site. Could you help us explaining if you could manage to download any dataset? Also, what kind of sequencing data do you want: single-cell RNa-seq, long read or bulk RNA-seq? Also if you could clarify if you have some preference for a programming language to automatize the process it would help too

Hello, What exactly you find tidious? If you search for "p53 AND "Homo sapiens"[orgn:__txid9606]" you get plenty of hits for NULL and WT. For example reads liked to this project. If you want fastq files, you can pull them from ftp server once you know IDs, here is an example: ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR835/SRR835937/

Hello, thank you all for such an immediate response. Ilrs, I want RNA seq data obtained from p53 WT, NULL and KO colon cancer cell lines, without any other modifications in the cell. I'm not bothered about programming language for further analysis now. I am actually new to this, and I don't know if I am collecting the right kind of data or not.  Kamil, these are shPFKFB4 silenced. I want normal vs experiment data from same experiment if possible. I might be making very silly mistakes, I apologise for the same. Thanks for your help.

Downloading from ENA is simpler: https://www.ebi.ac.uk/ena

I suggest you to try difference resources, including ArrayExpress.
A simple search at ArrayExpress returns multiple experiments which might be of your interest. Just a pointer, ArrayExpress is also synced with GEO.

I agree with Kamil on the downloading method.
Moreover, SRA will give you your file in the raw format which is usually a Binary Alignment Map (BAM) or Sequence Alignment Map (SAM)
Download SRA toolkit available online, to convert your raw reads to fastq format. Follow the steps using command line. It is a very simple command to convert your SRA file to fastq which uses the fastq-dump command.
Further reading: https://www.ncbi.nlm.nih.gov/books/NBK158900/
Hope it helps!

Hello Pawan, my main concern here it to pinpoint the data I require from the database. I am not able to find the data I require. I'm aware about fastq-dump. Thank you for your efforts.

I understand! Since you had asked if it was possible to download in fastq format, i had answered accordingly. Cheers

Thanks. Actually, when I'm using fastq-dump, all I see upon execution of command is a bunch of .cache files.Also, as huge files in GB scale take time to download, I would like to know if there is a way where we can track the progress of download?

Can you specify the size of your file?

Last file I tried downloading  was roughly 12GB. Max download speed on my network is 550kbps.

Well, downloading at this speed will take over 6 hours.


text:
How can I obtain the percentage gene expression per identity class in Seurat as further processible numbers (e.g. matrix)?
I am analysing my single cell RNA seq data with the Seurat package. I want to know if there is a possibilty to obtain the percentage expression of a list of genes per identity class, as actual numbers (e.g. a matrix) which I can write out to say an excel file. For example I want to be able to say gene X is expressed in 30% of cell A, 10% of cell B etc.
I can generate something comparable with:
DotPlot(object = my_object, genes.plot = "my_gene")

However the results are only graphic and I wish to have further processible numbers.
Furthermore:
AverageExpression(object, genes.use = "my_gene")

Produces expression values which I cannot transform to percentages.
I will be very grateful on any hints.

Did you read the tutorials of seurat? Why you cannot transform to percentage expression values?

@IIrs I will be very grateful for any hint. I am not able to figure a way how I can convert expression values for one cell type into percentage expression for that cell type. To say X% of my cell type C expresses this gene G.

This can be solved like this:
library(Seurat)

my_genes <- c("gene1", "gene2", "gene3")

exp <- FetchData(object, my_genes)

matrix <- as.matrix(colMeans(exp  > 0))*100

Glad you found your way! Perhaps explaining how did you find this solution (what did you look up and where) will help your future self or other people in the same situation. Also are you sure that FetchData retrieve the expession data? (Also which one, the raw data, the scaled data or the imputed data?)

You can get the table that is used to make the dot plot if you modify the DotPlot function to return it instead of the ggplot, and use the argument do.return=T.
To edit the function, the command is:
trace("DotPlot",edit=TRUE)

Then replace the last line "return(p)" by "return(data.to.plot)" and save the edit.
You can call the function and store the return in a variable.
The editing doesn't affect the source code, so the function will be untouched at your next R session.
Best,

Also potentially useful is the solution proposed here: https://github.com/satijalab/seurat/issues/1888
Basically, save DotPlot to an object, and your percentages are found within.


text:
How to extract metadata from NCBI's short read archive (SRA) for a few runs?
I wish to extract metadata from a list of runs on NCBI's short read archive. For instance, I'd like to extract the library name ("HS0798") from the following run info:
https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=run_browser&run=SRR390728
I have to extract metadata for a few hundred run numbers. Is there a existing command-line tool (sratoolkit, python, perl, bash, etc) that will extract the metadata about the run?
I could grep it out of a webpage dump, but I notice there are XML schemas provided on the site, which to me indicates that there's some sort of XML API available somewhere.

using xslt, see an old answer I wrote for biostar: https://www.biostars.org/p/69974/#69995

Entrez Direct is perfect for this sort of thing. You can do something like this:
$ esearch -db sra -query 'SRR390728' \
  | efetch -format runinfo \
  | cut -f1,12 -d,
Run,LibraryName
SRR292241,HS0798
SRR390728,HS0798

The runinfo table format has 47 fields including experiment accession, project accessions, bioproject and biosample identifiers and such. You may want to explore that a little bit more if you are interested in downloading metadata in addition to the LibraryName.
Also, if you prefer to parse XML output, you can append -mode xml to the efetch command to retrieve the runinfo table in XML format. xtract tool that comes with the Entrez Direct package can then be used to convert that into a table.

Definitely better than my screen grab. I installed the conda package for sratools, which doesn't have that esearch binary with it. I'll have to look what package comes with that.

Follow the installation instructions here: Entrez Direct. There is a package in Ubuntu repositories but it is quite behind from the latest version.

But it seems that the runinfo is not the same as the RunTable found at this page: https://www.ncbi.nlm.nih.gov/Traces/study/?. For example, go to the link and search for PRJNA301892. Under the "select" section (blue flag), there is a Metadata button. Is there a way to get that done?

The -q option is not valid with the version I just installed. -query should be used instead.

Thanks, I updated the answer to use -query now.

Here's a very crude script version, which should work without downloading fancy software tools. It will only work as long as the website stays in its current layout.
It extracts the particular <td> from the webpage, based on a nearby table link ID. and prints it out.
Note: This is brittle. It will only work for library names that fit in a single-line of the HTML, and will ignore some webserver errors. The content is also not html-decoded (so its output might show things like & if the field contained html-special characters).
#!/bin/bash

set -e -o pipefail

# pass the run as a single argument to the script
run="$1"

# extract library name field. second cell in visual order.
curl -o - "https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=run_browser&run=$run" 2>/dev/null | \
   grep 'a id="exp-accession"' -A 7 | \
   tail -n +3 | head -n 1 | \
   sed -e 's@^<td>@@' -e 's@</td>$@@'

Its got to be the SRA Toolkit
The way I would do it (a traditional Perl digger) doesn't work on the Tracer NCBI site, at least not readily and dumping the XML is likewise unlikely to succeed.

I've been told by the support team that there's a way to pull it out in csv form by GETting the following URL, replacing <SRRNUMBER>, with, e.g. SRR000001:
https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?save=efetch&db=sra&rettype=runinfo&term=<SRRNUMBER>

So this would be:
# download the file attachment
wget -nv -O foo.csv  "https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?save=efetch&db=sra&rettype=runinfo&term=SRR000001"

# or just get the output on stdout
curl -v "https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?save=efetch&db=sra&rettype=runinfo&term=SRR000001"

As an example, the output looks like this for SRR000001. It namely contains additional info on other runs from the same SRX:
Run,ReleaseDate,LoadDate,spots,bases,spots_with_mates,avgLength,size_MB,AssemblyName,download_path,Experiment,LibraryName,LibraryStrategy,LibrarySelection,LibrarySource,LibraryLayout,InsertSize,InsertDev,Platform,Model,SRAStudy,BioProject,Study_Pubmed_id,ProjectID,Sample,BioSample,SampleType,TaxID,ScientificName,SampleName,g1k_pop_code,source,g1k_analysis_group,Subject_ID,Sex,Disease,Tumor,Affection_Status,Analyte_Type,Histological_Type,Body_Site,CenterName,Submission,dbgap_study_accession,Consent,RunHash,ReadHash
SRR000001,2008-04-04 14:36:13,2015-04-07 21:54:15,470985,129515135,236041,274,298,,https://sra-download.ncbi.nlm.nih.gov/traces/sra27/SRR/000000/SRR000001,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,25D46C95842BCD2F235DDD48070341C4,5FC1B56C9A5332BD21A6FCC6713C3159
SRR000013,2008-04-04 15:15:44,2012-01-19 15:14:06,478304,131410350,316908,274,307,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000013,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,A576F73351E329617F91C92845ADE3BE,0EFE27254FB8790C80E4A8A3A3CCBF1D
SRR000019,2008-04-04 15:20:15,2012-01-19 15:14:06,537808,147928455,358612,275,347,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000019,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,5E42E66B2972A28C3C1B714499B3CE1E,A44F729DF562E11BDA4F950F171085B3
SRR000033,2008-04-04 15:26:19,2012-01-19 15:14:16,476737,131230950,316780,275,307,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000033,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,8D775F0FB1E78F2818166B4E6C7BC320,6C794108B0172AE98323C9A809BE892E
SRR000047,2008-04-04 15:34:57,2012-01-19 15:14:25,322279,88301883,214744,273,208,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000047,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,B4D7BC3DAD22AE58EADC18B2F05B2732,154C9532706F955A63CF685C4D71D9E6
SRR000049,2008-04-04 15:35:58,2012-01-19 15:14:34,506408,139029601,328595,274,336,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000049,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,60D02918CD7686DA084607BFD9D8E6A5,E9A732F2A8F36B8EB1CD38CA41BA8FB5
SRR000061,2008-04-04 15:41:48,2012-01-19 15:14:37,532884,146330433,351541,274,348,,https://sra-download.ncbi.nlm.nih.gov/traces/sra4/SRR/000000/SRR000061,SRX000007,SID2748,WGS,RANDOM,GENOMIC,PAIRED,3000,0,LS454,454 GS FLX,SRP000001,PRJNA33627,1,33627,SRS000100,SAMN00001583,simple,9606,Homo sapiens,NA18505,,,,,female,,no,,,,,454MSC,SRA000197,,public,8F6E1CA597BCAA19A852036775FA1FAB,2C035CCBD62808A01CEF6EB13BAEB40F

Providing an invalid SRR number will generate an empty HTTP 200 response.

Not a command-line tool but also useful is redesigned SRA Run Selector. Insert comma separated accessions and you can filter/download based on metadata.

You can get metadata for an entire project using pysradb:
> pysradb metadata SRP063732 --detailed

study_accession experiment_accession experiment_title                                                       experiment_desc                                                        organism_taxid  organism_name   library_strategy library_source library_selection sample_accession sample_title instrument           total_spots total_size   run_accession run_total_spots run_total_bases run_alias             sra_url                                                                                experiment_alias isolate                 dev_stage sex   tissue             treatment     BioSampleModel
 SRP063732       SRX1249861           Bisulfite-seq of rhesus macaque 26168_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 26168_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071080                    Illumina HiSeq 2000  89051187    11830170462  SRR2420581    89051187        17988339774     26168_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420581/SRR2420581.1                   26168_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249860           Bisulfite-seq of rhesus macaque 26148_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 26148_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071081                    Illumina HiSeq 2000  87335402    11763248173  SRR2420580    87335402        17641751204     26148_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420580/SRR2420580.1                   26148_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249859           Bisulfite-seq of rhesus macaque 25822_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 25822_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071082                    Illumina HiSeq 2000  83544130    9481342835   SRR2420579    83544130        13784781450     25822_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420579/SRR2420579.1                   25822_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249858           Bisulfite-seq of rhesus macaque 25154_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 25154_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071083                    Illumina HiSeq 2000  115476876   13003165573  SRR2420578    115476876       19053684540     25154_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420578/SRR2420578.1                   25154_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249857           Bisulfite-seq of rhesus macaque 23784_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 23784_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071084                    Illumina HiSeq 2000  101539017   13697691917  SRR2420577    101539017       20510881434     23784_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420577/SRR2420577.1                   23784_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249856           Bisulfite-seq of rhesus macaque 23779_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 23779_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071085                    Illumina HiSeq 2000  98394582    13326784419  SRR2420576    98394582        19875705564     23779_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420576/SRR2420576.1                   23779_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249855           Bisulfite-seq of rhesus macaque 23773_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 23773_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071086                    Illumina HiSeq 2000  84442355    11218164664  SRR2420575    84442355        17057355710     23773_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420575/SRR2420575.1                   23773_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249854           Bisulfite-seq of rhesus macaque 23764_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 23764_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071087                    Illumina HiSeq 2000  87197006    11291883405  SRR2420574    87197006        17613795212     23764_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420574/SRR2420574.1                   23764_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249853           Bisulfite-seq of rhesus macaque 23762_HeavyDrinker: Nucleus accumbens  Bisulfite-seq of rhesus macaque 23762_HeavyDrinker: Nucleus accumbens  9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071088                    Illumina HiSeq 2000  89455608    11917710767  SRR2420573    89455608        18070032816     23762_NaC_HD_files    https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420573/SRR2420573.1                   23762_NucleusAccumbens  Adult     male  Nucleus accumbens  HeavyEthanol  Model organism or animal
 SRP063732       SRX1249852           Bisulfite-seq of rhesus macaque 25811_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25811_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071089                    Illumina HiSeq 2000  88216999    10285962574  SRR2420572    88216999        14555804835     25811_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420572/SRR2420572.1                   25811_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249851           Bisulfite-seq of rhesus macaque 25790_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25790_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071090                    Illumina HiSeq 2000  94731241    10748979142  SRR2420571    94731241        15630654765     25790_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420571/SRR2420571.1                   25790_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249850           Bisulfite-seq of rhesus macaque 25742_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25742_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071091                    Illumina HiSeq 2000  113512952   17943763654  SRR2420570    113512952       22929616304     25742_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420570/SRR2420570.1                   25742_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249849           Bisulfite-seq of rhesus macaque 25240_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25240_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071092                    Illumina HiSeq 2000  81002663    9232986327   SRR2420569    81002663        13365439395     25240_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420569/SRR2420569.1                   25240_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249848           Bisulfite-seq of rhesus macaque 25184_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25184_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071093                    Illumina HiSeq 2000  110228135   12594727479  SRR2420568    110228135       18187642275     25184_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420568/SRR2420568.1                   25184_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249847           Bisulfite-seq of rhesus macaque 25157_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 25157_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071094                    Illumina HiSeq 2000  97249804    11026452312  SRR2420567    97249804        16046217660     25157_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420567/SRR2420567.1                   25157_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249846           Bisulfite-seq of rhesus macaque 21119_LowDrinker: Nucleus accumbens    Bisulfite-seq of rhesus macaque 21119_LowDrinker: Nucleus accumbens    9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071095                    Illumina HiSeq 2000  68899204    9241244635   SRR2420566    68899204        13917639208     21119_NaC_LMD _files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420566/SRR2420566.1                   21119_NucleusAccumbens  Adult     male  Nucleus accumbens  LightEthanol  Model organism or animal
 SRP063732       SRX1249845           Bisulfite-seq of rhesus macaque 26104_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 26104_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071096                    Illumina HiSeq 2000  84859783    9604170375   SRR2420565    84859783        14001864195     26104_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420565/SRR2420565.1                   26104_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249844           Bisulfite-seq of rhesus macaque 26082_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 26082_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071097                    Illumina HiSeq 2000  73346886    8248399265   SRR2420564    73346886        12102236190     26082_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420564/SRR2420564.1                   26082_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249843           Bisulfite-seq of rhesus macaque 25207_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 25207_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071098                    Illumina HiSeq 2000  77306973    12235827725  SRR2420563    77306973        15616008546     25207_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420563/SRR2420563.1                   25207_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249842           Bisulfite-seq of rhesus macaque 25407_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 25407_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071099                    Illumina HiSeq 2000  83985968    13355795218  SRR2420562    83985968        16965165536     25407_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420562/SRR2420562.1                   25407_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249841           Bisulfite-seq of rhesus macaque 24818_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 24818_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071100                    Illumina HiSeq 2000  72782928    8484613123   SRR2420561    72782928        12009183120     24818_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420561/SRR2420561.1                   24818_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249840           Bisulfite-seq of rhesus macaque 24698_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 24698_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071101                    Illumina HiSeq 2000  79414016    12627204820  SRR2420560    79414016        16041631232     24698_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420560/SRR2420560.1                   24698_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249839           Bisulfite-seq of rhesus macaque 24340_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 24340_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071102                    Illumina HiSeq 2000  84772576    13436053136  SRR2420559    84772576        17124060352     24340_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420559/SRR2420559.1                   24340_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal
 SRP063732       SRX1249838           Bisulfite-seq of rhesus macaque 22770_Control: Nucleus accumbens       Bisulfite-seq of rhesus macaque 22770_Control: Nucleus accumbens       9544            Macaca mulatta  Bisulfite-Seq    GENOMIC        other             SRS1071103                    Illumina HiSeq 2000  83310630    13224887149  SRR2420558    83310630        16828747260     22770_NaC_Ctrl_files  https://sra-download.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-3/SRR2420558/SRR2420558.1                   22770_NucleusAccumbens  Adult     male  Nucleus accumbens  Control       Model organism or animal


text:
wtdbg2: practical implications of k-mer fsize and psize choice
I am using wtdbg2 2.3 to assemble a human genome (sequenced on PromethION from a cell line). I filtered out reads with low average quality, and now I am trying to determine the parameters that will optimize overall assembly quality.
The wtdbg2 assembler has two k-mer related settings, documented as follows.
 -k <int>    Kmer fsize, 0 <= k <= 25, [0]
 -p <int>    Kmer psize, 0 <= p <= 25, [21]
             k + p <= 25, seed is <k-mer>+<p-homopolymer-compressed>

The requirement that k + p is <= 25** suggests that these parameters can be set simultaneously. However, in all of the presets for the various long read sequencing technologies, one of the two settings is always set to 0.
My questions are as follows.

Are k and p ever intended to be used simultaneously?
If so, how does the behavior change with -k 2 -p 19 versus -k 19 -p 2?
I'm familiar with some common implications of choosing k size for k-mer based analyses of Illumina reads. (Smaller k-sizes are more sensitive but also more subject to spurious matches. Larger k-sizes are more specific but also more likely to include a sequencing error.) Do similar implications apply for -k and -p with wtdbg2?
BONUS: Is there an explanation anywhere of how -k and -p are used? It looks to me like -k is related to k-mers in the traditional sense, and -p is a homopolymer compressed k-mer. Are seeds made by concatenating the regular k-mer and the compressed k-mer? It's possible I missed some important details in the preprint but I haven't been able to make the connection yet.

**Incidentally, the documentation here is wrong. During runtime wtdbg2 enforces a constraint of <= 23, not 25. This matches the description in the bioRxiv preprint stating that 46 bits (2 bits per base) are used to store each k-mer.

There is some description of homocompressed polymers (-p kmer) in this paper: https://arxiv.org/pdf/1708.01492.pdf


text:
How to do protein quantification with Peptideshaker output (TMT labeled, MS2)?
I have TMT tagged MS2 mass spec run processed with Peptideshaker. I saved the project in a peptideshaker format cpsx. Now, I would like to do the quantification with Reporter. I can load the file with Reporter GUI. Though Reporter is still in development and cannot be run from commandline yet. Are there other tools available that can run from Peptideshaker exported reports for downstream analysis? I tried to run Reporter from commandline version 0.7.16. I will try 0.7.17 which was just released.

It seems to me there should be something, but without knowing what you have tried we could repeat the same search you did. What search did you perform? Do you need it to be accessible via command line?


text:
Randomness in BLAST
So if you see the BLAST parameters it says

The Expected value E  is a parameter that describes the number of hits one can "expect" to get by chance when searching a database of particular size. It decreases exponentially as the score (S) increases. Essential, the E value describes the randomness background noise. For example an E value of 1 assigned to a hit can be interpreted as meaning that a database of the current size one might expect to see 1 match with a similar score simply by chance. Lower the E value the better."

How could there be something "BY CHANCE" it's a software right ?

The Blast E is the expected frequency of obtaining false positive. This will depend on your query size (number of nucleotide or amino acid residues) and the size of the database. With a short query and a large database you are more likely to have a sequence in your database that matches your query by simple chance. In summary, the smaller the query and the larger the database the greater the chance of a spurious result and that is what E is measuring.

Sure, so you start out with what is called a bitscore, which is a normalized to the score calculated from the alignment between the 2 seqs which depends on the following equation. It is independent of database size:

(lambda * S - ln(k))/(ln)2

Then, the p-value of a local blast is just:

1/2^bitscore

So if your bitscore is 15, you need 1/32768 alignments before you will get a score as good or better (highly similar sequences) BY CHANCE ALONE. This addresses your "BY CHANCE" question.
Earlier I said the bitscore was independent of database size. The E value is just the p-value above normalized to the database size (so it is dependent on the db size) by the following equation:

(query length * database length * p-value)

which simplifies a bit to:

E = (query length* (db length/2^bitscore)

I hope this helps!

Calculation of E for Blast is a good question.
Two methods:

Use a Poison or binomial distribution.
Use randomization

I think its no. 1 because no. 2 would use too much computational power.  I did know but someone borrowed by O'Reilly book on the subject.
One of the contributors here will know, because I suspect their genomics developers.


text:
Questions regarding Nanopore sequencing analysis
I am new to Nanopore sequencing analysis. I have a couple of questions regarding it which are as follows:

How do I know if my fast5 file is multiread or single read file?
Is there a way to combine all the fast5 files into a single fast5 file?
How can I extract the flowcell type and sequencing kit information from the fast5 file? Is there a command that I can use to extract it? I know I can use the following command to retrieve it:

h5dump -a /UniqueGlobalKey/context_tags/sequencing_kit -a /UniqueGlobalKey/context_tags/flowcell_type <fast5 file>

However, it works well for single read file and not for multiread files. Please correct me if I am wrong. Thanks!

[Note: As of 2024, ONT now uses pod5 as the standard output format from their MinKNOW sequencing program, and all pod5 files are "multi" pod5 files.]
1. The quickest way to know if it's a single or multi fast5 is the file name. If it mentions channels and read numbers (e.g. ..._read_9_ch_471_..., then it's a single read file. Alternatively, the multi-read files have only read IDs at the top level (i.e. read_<32-character-string-with-additional-dashes>). Example:
$ h5ls <read_file_name>.fast5 | head
read_0115573a-042a-42a3-b4e2-f0338b4a6c66 Group
read_015b870d-8a03-4f22-a73e-fb9f833b94c7 Group
read_0167393a-dede-4a27-91af-494622b63261 Group
read_01e658de-853b-4f58-93ad-87fbd735254c Group
read_0202acc2-33f4-4d29-9dac-c7f80c758c86 Group
read_022d68f0-6f8e-4124-bfb0-11262efa4a67 Group
read_0258a3c9-4ffe-42ab-a42d-89466fcefbf5 Group
read_032da6db-e84a-43ec-8c25-bcaead992066 Group
read_03383754-be43-420b-87fc-60c2eb833963 Group
read_036b4824-bb0a-42dd-ade0-29d89112ff07 Group

2. ONT provide a single_to_multi program for converting to the multi-read format. Details can be found on this github page. Here's a quick example of how it could be run:
$ single_to_multi_fast5 --input_path fast5 --save_path fast5_multi \
    --filename_base multi_file_myRunID --batch_size 4000 --recursive

3a. For multi-fast5 files, the kit/cell attributes are in the context_tags group of the file (because a multi-fast5 file could contain fast5 elements from multiple different sequencing runs with different run parameters):
$ h5dump -g /read_<readID>/context_tags multi_file_myRunID_0.fast5 | \
  grep -A 11 -e 'sequencing_kit' -e 'flowcell_type'
   ATTRIBUTE "flowcell_type" {
      DATATYPE  H5T_STRING {
         STRSIZE 11;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "flo-min106"
      }
   }
--
   ATTRIBUTE "sequencing_kit" {
      DATATYPE  H5T_STRING {
         STRSIZE 11;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "sqk-pbk004"
      }
   }

3b. For single fast5 files (at least, in the most recent version of them), they are found in the UniqueGlobalKey group:
$ h5dump -g /UniqueGlobalKey single_file_read_X_ch_Y_strand.fast5 | \
  grep -A 11 -e 'sequencing_kit' -e 'flowcell_type'
      ATTRIBUTE "flowcell_type" {
         DATATYPE  H5T_STRING {
            STRSIZE 11;
            STRPAD H5T_STR_NULLTERM;
            CSET H5T_CSET_ASCII;
            CTYPE H5T_C_S1;
         }
         DATASPACE  SCALAR
         DATA {
         (0): "flo-min106"
         }
      }
--
      ATTRIBUTE "sequencing_kit" {
         DATATYPE  H5T_STRING {
            STRSIZE 11;
            STRPAD H5T_STR_NULLTERM;
            CSET H5T_CSET_ASCII;
            CTYPE H5T_C_S1;
         }
         DATASPACE  SCALAR
         DATA {
         (0): "sqk-pcs108"
         }
      }

[note that this file structure may change in the future]


text:
R studio failed to install the latest development version of hdf5r from Github (win10)
I install hdf5r following https://github.com/hhoeflin/hdf5r， it shows two methods:

When I use the second method to install from github,it showed many errors：

The downloaded source packages are in
‘C:\Users\dell\AppData\Local\Temp\RtmpI7TP9f\downloaded_packages’
√ checking for file 'C:\Users\dell\AppData\Local\Temp\RtmpI7TP9f\remotes8c87faaf3d\hhoeflin-hdf5r-cc278c1/DESCRIPTION'

* preparing 'hdf5r': (4.9s)
  √  checking DESCRIPTION meta-information ...
* cleaning src
* checking for LF line-endings in source and make files and shell scripts (676ms)
* checking for empty or unneeded directories (1.5s)
* building 'hdf5r_1.0.1.tar.gz'
  Warning: file 'hdf5r/cleanup' did not have execute permissions: corrected
  Warning: file 'hdf5r/configure' did not have execute permissions: corrected

* installing _source_ package 'hdf5r' ...
  checking hdf5 headers and libraries
  Environment variable LIB_HDF5 is not set. Attempting to download HDF5 files from repository.
  试开URL’https://github.com/mannau/h5-libwin/archive/master.zip'
  Content type 'application/zip' length unknown
  ** libs
  错误: (由警告转换成)this package has a non-empty 'configure.win' file,
  so building only the main architecture
* removing 'D:/soft ware/R-3.5.2/library/hdf5r'
  In R CMD INSTALL
  Error in i.p(...) :
  (converted from warning) installation of package ‘C:/Users/dell/AppData/Local/Temp/RtmpI7TP9f/file8c863b6795/hdf5r_1.0.1.tar.gz’ had non-zero exit status
  In addition: Warning messages:
  1: In untar2(tarfile, files, list, exdir) :
  skipping pax global extended headers
  2: In untar2(tarfile, files, list, exdir) :
  skipping pax global extended headers

But when I use the first command:install.packages("hdf5r"), it worked:

Can anyone help me? I search many times, but I could't find the answer.

Usually installing a package from github means installing the developing version, so it might have errors, you probably should install the stable release stored in CRAN by using the install.packages("hdf5r"). Also it is incompatible to have two versions of the same package at the same localization. Decide either if you need both version and where do you want to install them if you want both versions. Are you sure you need the development version of this package?

Try running this:
source("https://bioconductor.org/biocLite.R")
biocLite("hdf5r")

Just a different method of installing packages besides github. It is more centred around bioinformatics packages as well so it is pretty useful.

Bioconductor packages now use 'install.packages("BiocManager")' / 'BiocManager::install("")'


text:
failed to install hdf5r
I failed to install hdf5r using this command provided by https://satijalab.org/seurat/mca_loom.html: devtools::install_github(repo = 'hhoeflin/hdf5r'),below is the code,can you help me ?


text:
How to get fasta alignment file from SAM/BAM file?
I am not talking about consensus sequence, I know how to get consensus sequence using mpileup in samtools/bcftools. As I understand , SAM/BAM files are basically sequence alignment format so it's  natural to expect a straightforward way of conversion between them. Sadly I have yet to find one. Vcf-kit allows for conversion between vcf file and fasta alignment, but the fasta alignment contains only the SNPs. I need the whole alignment.  Does samtools/bcftools provide option for such conversion? If not, is there any other tool?

I am not sure what you mean by "fasta alignment file". If you mean a multi-sequence alignment (MSA) in the fasta format, you can't get that because SAM keeps pairwise alignments only and doesn't align inserted sequences. Even if you don't care about inserted sequences, a MSA in fasta is far to big to be practical. Alternatively, by "fasta alignment file", you could mean pairwise alignment, but it is still impractical to output every pairwise alignment in a separate fasta file.
There are tools to convert SAM to BLAST-like format if that is what you want. The following shows one of them. I believe there are other tools easier to use, but I forget what are they.
Anyway, to install:
wget https://raw.githubusercontent.com/lh3/minimap2/master/misc/paftools.js
curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` k8

Then you can convert your ALN.sam with:
./k8 paftools.js sam2paf -L ALN.sam | ./k8 paftools.js view -l0 - | less -S

Or generate a SAM stream from ALN.bam and pipe it with:
samtools view -h ALN.bam | ./k8 paftools.js sam2paf -L - | ./k8 paftools.js view -l0 -

Here is an example output:
>MT_orang   16499   0   16025   +   MT_human    16569   576 16569   13773   16095   60  tp:A:P  mm:i:2150   gn:i:172    go:i:101
Ref+:        577 GTTTATGTAGCTTACCTCCT---CAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTAAGATTACACATGCAAGCATCCCCGTTCCAGTGAGTTCACCCTCTAAATCACCACGATCAAAAGGAACAAGCATCA
                 ||||||||||||||  | ||   ||||||||| |||||||||||| | ||||||| |||| | |||||||||||||||||||||||||||||||||||||||||||||||| || |||||||||||||||||||||  |||||||| || ||||| || ||||   ||  || |||| |||||||||
Qry+:          1 GTTTATGTAGCTTA--TTCTATCCAAAGCAATGCACTGAAAATGTCTCGACGGGCCCACA-CGCCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTGAGGTTACACATGCAAGCATCCCCGCCCCAGTGAG-TCGCCCTCCAAGTCACTCTGACTAAGAGGAGCAAGCATCA
//

You can change -l0 to -l80 to print alignments in multiple lines.

I meant "aligned fasta format" which I think is MSA. Beside general curiosity, my immediate need for the format is as the input file for gubbins.

gubbins seems to be a program for the fiel of phylogenetics. So the input file must be an alignment file between your species and not the alignments of the reads for one sample.

@AhmedAbdullah you should probably ask a new question, explaining your final objective. Going from reads aligned to a genome (bam) to a multiple sequence alignment you want to use for phylogenetic analysis makes no sense. If you explain what you are trying to achieve in a new question, we should be able to help.

@terdon I am trying whether this works or not. I have 48 SAM files corresponding to 48 samples. I converted them to 48 ordered BAM file. and the BAM fiels to a single vcf file using : bcftools mpileup  -d 1000 -Bf R.fa *.bam  | bcftools call --skip-variants indels --multiallelic-caller --variants-only --samples-file samples.txt -O z -o aln.vcf.gz. vcf-kit can generate aligned fasta for only SNPs from this vcf file (vk phylo fast . The problem is alignment is only of SNPs, I was wondering whether one can generate full alignment from the (SAM files or the subsequent vcf file)

@AhmedAbdullah you are almost certainly approaching the problem in the wrong way. Please ask a new question explaining exactly what you need. What you're describing here makes no sense.

cool...did not know that paftools can print out ref vs query like this. another tool that can do similar is https://github.com/mlafave/sam2pairwise

InIGV you can visuals the alignments of individual reads to the reference as contained in the .bam file. I guess what you want to is get that information into a file containing all the (paired) read sequences aligned to a particular region, with padding for the missing data. That way you could inspect all SNPs or sequencing errors, check possible haplotypes across a few SNPs, use these for phylogenetic analysis etc.
I am also looking for such tool and best I could find so far is an answer on BioConductor. It refers to some code in R that can be used for this purpose:  genomicalignments::stackStringsFromBam().

MSA fasta file can be generated from SAM,
but it is almost impractical to generate this HUGE file.
If you still insist to generate the file, then please check this.
https://sourceforge.net/projects/sam2fasta/files/


text:
How to cluster the human genes by pathways/system-biology/metabolic properties?
I would like to make a chord plot from my data. I have a list of genes that according to my experiments are divided into 64 clusters of enrichment pattern.
I would combine my 64 clusters with a second clustering. A clustering of all the human genes into subgroups of "interacting" genes (e.g. both are pathway, co-expressed) will allow me to combine the two clustering data on a chord plot.
In the analogy of the following graph the continents will be my 64 clusters of the genes and the edges/chords will be assigned based on the biological function  clustering. I guess that such a functional clustering can be done ore has been done using platforms such as: KEGG, genescf, or reactome.org.

Is there a biological pathway clustering of genes for lymphoblastoid or breast cell lines?
Are there any suggestions how to build the clustering between genes whose functions are strongly correlated in the cell?

Edit: the question is about getting the clustering from a system-biological perspective not how to plot chord diagrams.
Is there such a metric to quantify how likely two genes are correlated in their enrichment, function etc. For example, using STRING we can see that PIK3CA and PTEN are more co-functioning than PIK3CA and SF3B1. How can I do something like String on genome wide scale.

Its difficult to understand the experiment, so not easy to recommend something other than generic correlation/covariance

@MichaelG. Would like something that measure the distance between two genes on the multigraph of pathways.

To measure if two genes are functional similar I developed the BioCor package. It calculates a similarity score between genes by the measuring the amount of pathways  shared. However, it doesn't take into account the strength of the connection between two proteins. You could multiply the similarity by a value between 0 and 1 according to the number of evidences of the relation between the two proteins.
It is optimized for large amounts of data, you could even perform them in parallel, but I doubt you find information about pathways for most of the genes (In my informal test, I think that I could calculate similarities for about 1600 genes not the ~20000 of the human genome).
When I am interested in co-expression data I calculate the co-expression by correlation and add the similarity of the genes calculated with BioCor. The resulting distance matrix can be used for clustering. Maybe the section about using weighted co-expression network analysis and BioCor can be helpful.
To have a "biological pathway clustering of genes for lymphoblastoid or breast cell lines", you would need to have a specific signature/pathway of the cell lines. Currently I am not aware of accepted standard signature database. The single cell atlas could be useful to find expression of your cells, but I don't know if they provide a signature for cell line

Is there a python version for it?

No, but you can use it in python by using rpy2

Production of a chord diagram is described in detail here. You use the R package with either,
library(chorddiag)
library (circulize)

The first level clusters are predefined and the second layer of clustering is what you need to work out. The obvious stat is correlation, secondly (better approach IMO) you might look at covariance, although I don't exactly know what this "interaction" constitutes, which makes it difficult.
My personal view is the diagramatic representation doesn't really matter, its the robustness of the stats behind the diagram that counts.
I was a bit surprised to learn ggplot2 doesn't do circular diagrams :-(

Lymphoblastoid Cell lines

I think you need to clarify the experiment on these cells lines, I assume its a time course, i.e. time series.

At a guess this is RNA expression, if you are using a chip then its between predetermined genes (i.e. spotted on the chip). If it is RNA seq then it can be anything. I suspect you are treating the cells and monitoring changes in RNA expression as a function of their biological (immunological?) activity. If the basic measure is time, i.e. a time series - thus two genes being expressed in the same unit time the best opening statistic is covariance. If all you want to do is make a chord diagram out of it, I think thats fine, this will form a 2x2 variance matrix for which any standard clustering will work (Jaccard distances [there will be better ones]) or direct plotting from the matrix. If you looking at unique expression behaviour to attempt to unlock some sort of pathway then you are into more serious regression time series (autoregression that sort of thing - this is complicated). If its oscillations then its a separate branch of stats based on sinosoids (sine waves) called periodicity analysis. HOWEVER, if the interaction is genetic similarity (gene families) e.g. mRNA then you want a phylogenetics measure, which if its clustering neighbor-joining is fine.

My personal view if it is the experiment I think, I would,

Construct a matrix of the covariance of the (I assume) time series
Make sure the covariance is standarised between -1 to 1 (Pearson
approach), R-cran.
Assign a tight cut-off (tricky because there's no "hard cut-off") but say 0.7 (you need to experiment a bit)
Plot the raw covariance on the chord diagram.

Easy, simple and highly visual.

For my experiment, I know how to cluster it. My question is independent! I would to find a metric that assign a distance to two genes according to their functions. If they happen to function together in the same context/pathways they will be close. For example, the proteins of NRAS and PIK3CA will be closer than NRAS and SF3B1 (http://version10.5.string-db.org/cgi/network.pl?taskId=KbO2wnULZxKQ)
I am looking for a resource to be able to build such thing on a genome wide scale.