text:
Unable to Install PGDSpider using Ubuntu and wine
I'm running wine v 2.0.3 using ubuntu v. 16.04.2 on a Mac. I'm trying to install a program called PGDSpider2 from here. Essentially the program is supposed to convert different file formats used for bioinformatics.
I'm running it in a 32-bit wine server and I have the Java V.9 and the Microsoft.net packages downloaded on that same server. I've been getting this error.
I have tried the installation on a clean bottle but get the same following error, any thoughts?
parallels@ubuntu:~/.wine32/drive_c/Program Files/PGDSpider_2.1.1.3$ WINEPREFIX=~/.wine32 wine PGDSpider2.exe
fixme:heap:RtlSetHeapInformation (nil) 1 (nil) 0 stub
fixme:process:SetProcessDEPPolicy (1): stub
fixme:heap:RtlSetHeapInformation (nil) 1 (nil) 0 stub
fixme:process:SetProcessShutdownParameters (00000380, 00000000): partial stub.
fixme:process:GetProcessWorkingSetSize (0xffffffff,0x73e94c,0x73e950): stub
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e6fc,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
parallels@ubuntu:~/.wine32/drive_c/Program Files/PGDSpider_2.1.1.3$ fixme:win:EnumDisplayDevicesW ((null),0,0x39edcf0,0x00000000), stub!
fixme:win:EnumDisplayDevicesW ((null),1,0x39edcf0,0x00000000), stub!
X Error of failed request:  BadMatch (invalid parameter attributes)
  Major opcode of failed request:  156 (GLX)
  Minor opcode of failed request:  5 (X_GLXMakeCurrent)
  Serial number of failed request:  241
  Current serial number in output stream:  241
k
fixme:advapi:RegisterEventSourceW ((null),L".NET Runtime Optimization Service"): stub
fixme:advapi:DeregisterEventSource (0xcafe4242) stub

According to this forum it is likely a problem with the wine configuration of the video card drivers or opengl. I would suggest posting on a forum for wine as I don't think this is an issue with PGDSpider.

@gringer the community consensus is that questions asking for help installing bioinformatics software are on topic since they're often non-trivial to install. I agree that this is a complex case, and it might be easier to find an answer on sites like [apple.se], [ubuntu.se] or [su], but that does not make it off topic here.

Mollie, while this is on topic here as I explain above, you might want to consider deleting it and asking instead on either [apple.se] or [ubuntu.se]. In either case, please [edit] to explain whether your Ubuntu is actually running on your mac hardware or if it is running using the parallels virtualization tool. It looks like the latter and that is relevant.

Also, why in the world are you trying to run the Windows version of this? The website you linked to provides both a MacOS and a Linux version. Either of those would be simpler to use than trying to run it through Wine. Don't make your life harder than it needs to be!

Don't use Wine at all! You seem to be running an Ubuntu installation inside a MacOS system using parallels. And then a Windows emulation layer (wine) on top of that. This is a very complex system.
However, the software you are trying to install is written in Java and can run perfectly well on any of the operating systems you have available. So there is absolutely no reason to even try to run it via Wine. Just run it natively on either your MacOS or your Ubuntu. As explained in the page you linked to:
Execute PGDSpider GUI:

Windows:
  execute the file PGDSpider2.exe to start the program
Linux:
  execute the command ./PGDSpider2.sh to start the program
Mac and others:
  execute the command java -Xmx1024m -Xms512m -jar
  PGDSpider2.jar to start the program

So ignore the Windows instructions and pick either one of the others.


text:
Canu assembly not making a single consensus?
I've downloaded reads from this BioProject. Using canu with default parameters (no correction), I've got 4 contigs, none of which really look like the reference plasmid here.
The command I used was:
canu -p ip40a -d ip40a_assembly -useGrid=false -maxThreads=6 -genomeSize=175k -nanopore-raw reads.fastq

Granted, the original publication for this sequence used minimap/miniasm but my questions are:
A. Why do I get multiple contigs? Doesn't canu generate a single assembled sequence?
B. This is a plasmid but canu says suggestedCircular=false. What's going on here?

Looking that the link you provided only 30% of the reads appear to be e. coli.  I would suggest filtering the reads for those that align to e. coli and see where that gets you.  Did the original publication describe the method they used for data cleaning?
A:
https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR3385327

Unidentified reads: 68.82%
Identified reads: 31.18%
cellular organisms: 30%
Bacteria: 30%
Proteobacteria: 28.09%
Gammaproteobacteria: 19.7%
Enterobacterales: 13.3%
Enterobacteriaceae: 12.54%
Vibrionales: 0.01%
Viruses: 1.17%

B:
Without a good assembly canu won't be able to determine if the assembled genome is circular.  I would ignore this until you have a good example.

Sadly the original publication only seems to say "assembled by miniasm/minimap" (along those lines) and I'm trying to figure out the pipeline as we speak -- there's no indication on whether the reads were error corrected, etc.

Also, how would I filter the E. coli-identified reads?

I would align the reads to the plasmid reference and see how many align and what the coverage of the reference looks like.  Are there any regions of the reference that are uncovered.  how even is the coverage across the reference.  These are things that can help you understand what to expect from the assembly.

It seems like it's actually the opposite that @thestatnoob needs to do: exclude E. coli reads, instead of selecting for them

Correct he has the plasmid reference he is trying to assemble.  He needs to select for those reads assuming the alignment coverage shows an assembly is even worthwhile.

The E. coli genome is closer to 4.6 Mb in size; if there are a large quantity of host reads in there, then they'll take over the reads and there won't be enough for your target plasmid. By default, it will only try to assemble from a maximum of ~40X coverage of the specified genome size (higher values do not add any more information for assembly).
It's possible that miniasm will take all reads into account, rather than just a portion of reads. You may need to filter out the E. coli reads for this assembly to work properly.
For what it's worth, Canu's default parameters do include a correction step.

Yep, but the reads aren't for an E. coli genome. It's for a plasmid within E. coli, which is definitely within the 170-180k region (in fact, I got the number from the ENA...). I'm basically just trying to reconstruct the consensus plasmid that the authors have, using canu instead of miniasm.

source: Stack Exchange

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

sort: votes


text:
What would be a reasonable healthy data set/control group to compare with a cell line from TCGA?
Assume I want to run a simple machine learning of healthy/non healthy classification, how can I use TCGA where it has only data of people with cancer. What do usually people do for example for study of gene expressions?

Depends. What do you want to study? Genome architecture? Gene expression? Cell structure? Metabolic rate? Variants? Usually, you would compare the healthy to non healthy sample from the same individual but we can't really answer unless you specify what sort of analysis you want to do.

@terdon gene expression. But TCGA doesn't contain healthy cells of the same person.

@terdon my final goal is to do metabolic analysis, but I need to do gene expression first.

In cancer research with gene expression profiling, you can find molecular subtypes with machine learning. These subtypes are usually not seen with the standard pathological methods. Here is a famous study where they define a different subtype of breast cancer, that respond differently to therapy. With this gene signature they can identify patients that do not benefit from a certain therapy, but would survive without it (you save them a lot of side effects which come with the therapy).

Consider that is hard to justify to an ethic committee to remove a healthy tissue with an operation for the purpose of acquire data.

Now that I was looking for an analysis in the TCGA data, AFAIK it has data of the same people but from healthy regions (near the tumoral regions), which is usually considered as healthy. I suspect that this samples are the ones used for ML. In the TCGA datasets there is a variable  type used to encode if the sample is a "tumor" or "normal".

@Llopis that what I thought, but I didn't see the healthy tissue sample there. Can you post an example.

This was an analysis for my master. See Figure 3.A. I didn't download the data, which was provided as a summarized experiment object of Bioconductor. But I would expect to be downloadable somehow.

@Llopis thanks. Feel free to add it as an answer. Did you write you thesis in jupyter or just the SI part?

That is only a project for one subject and it was written in rmarkdown. My master thesis was like an article and I used MS Word. I still have to learn  jupyter.

In the TCGA datasets there is a variable (type) containing the information if the sample was from the tumor or from an adjacent region, which is usually considered as healthy.
As this samples come from the same patient and the same tissues they are used in the analysis for the comparisons and classifications.

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

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


text:
Interpreting 0x200 flag in bwa-mem alignments
I am looking at the bwamem.h code in http://github.com/lh3/bwa and found that BWA-MEM will give flag 0x200 to what it calls MEM_F_SOFTCLIP:
#define MEM_F_PE        0x2
#define MEM_F_NOPAIRING 0x4
#define MEM_F_ALL       0x8
#define MEM_F_NO_MULTI  0x10
#define MEM_F_NO_RESCUE 0x20
#define MEM_F_REF_HDR   0x100
#define MEM_F_SOFTCLIP  0x200
#define MEM_F_SMARTPE   0x400
#define MEM_F_PRIMARY5  0x800
#define MEM_F_KEEP_SUPP_MAPQ 0x1000

But looking at the explanation in this website: https://broadinstitute.github.io/picard/explain-flags.html
It labels 0x200 as read fails platform/vendor quality checks.
If we have reason to believe that soft clipping can happen for some of the alignments, should we just ignore the fact that IGV calls this fail QC? Does the failQC tag come from the SAM/BAM specification or is this something IGV came up with?
The SAM/BAM specification labels 0x200 the same way as the Picard flag explaining website:

The question remains as to why BWA-MEM uses the 0x200 flag for softclipping notation rather than Fail QC.

bwa mem is using exactly the same meaning of the 0x200 flag as every other program, including picard. Don't blindly assume that that entry in the header file relates to the output flags in the SAM output. If you search the code, you'll see that bwa is using this internally as a flag for "this alignment is soft-clipped", which it handles differently for primary and secondary alignments (hard-clipping the latter).
For what it's worth, bwa mem (at least the version I just tried) doesn't set the QC failed flag. This behavour is actually pretty common among aligners, since Illumina has changed how they represent QC failure at least once.

I don’t understand this answer, as I think you’re contradicting yourself: you say that BWA-mem handles 0x200 “as every other program, including picard”, implying that it uses it for failed platform/vendor QC. But then you go on to say that it uses it for soft-clipped alignments (exactly as implied by the name in the header!). Well, which is it?

If it sets it then it has the same meaning. I just don't know off-hand if it sets it at all.

The connection between soft clipping and a read level flag is interesting.  Soft clipping occurs at a per base level and is denoted in the CIGAR string.  The failQC flag is applied to a read and usually comes directly from the sequencer. It is likely that BWA just named their internal var weird.

Does the failQC tag come from the SAM/BAM specification or is this something IGV came up with?

The failQC flag is part of the Sam Spec.

If we have reason to believe that soft clipping can happen for some of the alignments, should we just ignore the fact that IGV calls this fail QC?

If you expect soft clipping I would suggest trimming that sequence from the reads prior to alignment.  Also you can look at read in the input reads and see if the failQC flag is set prior to alignment to understand where it is getting set.  IGV is just displaying what is in the BAM file not making any decisions about those reads.

I wouldn't so much say that the variables are oddly named so much as it shows from the bit usage for flags (be they internal or external) that the same author wrote both bwa and SAM.


text:
About CAFA experiment
I have a question about CAFA1 and CAFA2 experiment.
CAFA is an experiment that compare computational methods predicting protein functions. In the paper discussing CAFA2, there is an image comparing 2 methods (CAFA1 versus CAFA2):

How do I interpret the numbers and colors (dark and light)? What do they mean?

Note that in the article published that figure is different from the preprint!
See the figure:

As they aren't substantially different in display other than showing the two set of predictions that are comparable between CAFA1 and CAFA2 I explain them below:
The colors are the margins of $F_{max}$ they refer to the equation found in the section "Comparison between CAFA1 and CAFA2 methods" under methods:
$$\delta(m_2, m_1) = \frac{1}{B}\sum_{b=1}^{B}F_{max}^{b}(m_2)-\frac{1}{B}\sum_{b=1}^{B}F_{max}^{b}(m_1)$$
Where B is the 10000 bootstraps they did and $m_1$ and $m_2$ are the CAFA1 and CAFA2 methods.
Where $F_{max}$ is explained in the article (and in the preprint) as:

A perfect predictor would be characterized with $F_{max} = 1 $

So, the higher the (positive) difference is, the better precision and recall (See wikipedia) are for those methods. Note that the scale is from -0.1 to 0.1 not from -1 to 1 so the differences are quite small (One could argue that the log10 would be better to see those small differences)
One can appreciate that in the bar plots at the right, where the differences between each dataset are small in the B section, so the color is white, and in the C panels they have brighter color.
The numbers in the boxes are the percentage that a CAFA2 method did better than a method in CAFA1 for another dataset in the "intersection of the two GO snapshots" as noted in the methods section. (But it is also not clear to me)


text:
What do 0s in stageR mean?
I'm running stageR as part of a pipeline based on DRIMSeq and wanted to know what it means when the transcript level padj is equal to 0, particularly when the gene is non-significant. I'm pretty sure that in some cases if the padj is 0 it means that the chance of it being different due to chance is 0. There are a lot more 0s than I would expect in sites where the gene padj is non-significant. One example, I have a gene padj that's >0.9 but both of the sites I'm looking at have a padj of 0.
I have a couple of suspicions about why this is: either because it defaults to this if the comparative difference is or it uses 0 as a substitute to NA. Am I thinking along the right lines? Are the 0s real or are they likely to be a quirk of the package?

In stageR, we automatically return stage-wise adjusted p-values of zero for transcripts coming from genes that were found to be significant and only have 1 transcript, when testing for differential transcript expression (DTE).
Similarly, when testing for differential transcript usage (DTU), we return p-values of zero for transcripts derived from a significant gene consisting of no more than two transcripts (based on the tx2gene object you provide).
Transcript-level p-values for non-significant genes should not occur. We provide NA's for genes with only one transcript for DTU, and by default you should get NA p-values for non-significant genes, since these genes do not proceed further than the screening stage.
Could you provide some more details on the analysis? Which multiple testing correction have you selected?
Would it be possible to provide a reproducible example for your analysis?

I'm using DTU and I've set the alpha to 0.5 as we were hoping that some of the genes which were non-significant overall might have some features that were significant. I'll see if I can make a reproducible example in a bit.

I'm really sorry it took so long to get back to you, I ended up getting distracted by something else. It turns out that the issue was probably elsewhere within what I was doing.

No problem. Let us know if you have any remaining questions or issues.


text:
RcppArmadillo installation error
When I am trying to install the RcppArmadillo Bioconductor library in R, it is giving me the following error:
/usr/bin/ld: cannot find -lgfortran
collect2: error: ld returned 1 exit status
make: *** [RcppArmadillo.so] Error 1
ERROR: compilation failed for package ‘RcppArmadillo’
* removing ‘/home/nikita/R/x86_64-pc-linux-gnu-library/3.2/RcppArmadillo’

I am using Ubuntu 12.04, R 3.2.5 and g++-5.4. Any help would be greatly appreciated.

It can't find gfortran. Is it installed?

I'm voting to close this question as off-topic because it is not a question about biology, but rather about software installation.

In a way you are right. The issue is that this question is more appropriate here rather than on stackoverflow because bioinformaticians are here

I tried installing gfortran with sudo apt-get install gfortran but it is already there, so nothing happened.

Do you have r-base-dev installed? It should get the necessary headers.

@NikitaVlasenko There is a bioinformatics.se...

I found a similar question on stackoverflow and the solution was to create a symlink as described by kevin there:
sudo ln -s /usr/lib/x86_64-linux-gnu/libgfortran.so.3 /usr/lib/libgfortran.so

This feels like a dirty hack — better fix your GCC installation.

On Ubuntu 14.04 there is a direct route to installing RcppArmadillo. You could perhaps try the same syntax.
sudo apt-get update
sudo apt-get install r-cran-rcpparmadillo

It is giving me Unable to locate package r-cran-rcpparmadillo


text:
sort a fasta file containing the Oxford Nanopore Technologies (ONT) header by sequencing start_time ascending
I have basecalled ONT reads and converted them to multifasta. The multifasta contains the original ONT headers in this format:
3ebb56cd-3671-4abd-b1ac-0c759bt068d0 runid=eb6214851489c8e00eb0dcbd00d737f7dddxxxx read=1078 ch=422 start_time=2017-11-21T20:35:34Z
By checking the sequencing output as fasta, the sequences are not ordered by start_time.
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATG
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATG
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTTTCCTGTCGCTCTATCTTCAC

I learned how to extract header IDs and time_start from ONT fastq files by using this answer.
I'd like to have a similar approach but with fasta.
Yet, I would like to create a new fasta file where the sequences are written in an ascending order as they were sequenced (start_time) as follows:
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATG
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTTTCCTGTCGCTCTATCTTCAC
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATG

The answer you link to is about fastq, but you talk about fasta in your question and show a line that isn't either. You have received one answer which assumes your input is fasta and one which assumes your input is fastq. Please [edit] your question and include an example of your input file so we can know which is correct. Ideally, show us a few lines of the original file and then show us the output you would like to see from those lines.

Thank you @terdon. I have edited the question to clarify!

Perfect! Thank you, now it's perfectly clear! By the way, this site supports a specific version of markdown syntax for formatting, you can see the details here: formatting tools.

Edit: After doing some reading, if it can be assumed that all the timestamps are ISO8601 and of the same representation of ISO861 (there are a few), then they can be sorted lexicographically, directly, without converting to epoch-time integers. So I have edited this answer to add an answer without conversion steps:

Here's one approach that uses the timestamps directly:
(1) Extract the start times from the headers and sort them to another file, in ascending order:
$ awk '(/^>/){ match($0, /start_time=(.*)/, a); print a[1]; }' my.fa | sort > starts_iso8601.sorted.txt

(2) Use that file of ordered ISO8601-formatted starts to write out FASTA records in your sort order:
$ while IFS='' read -r line || [[ -n "$line" ]]; do awk -v line=$line -v RS=">" '{ m = match($0,line); if (m) { printf(">%s", $0); } }' my.fa ; done < starts_iso8601.sorted.txt > my_ordered.fa

Quite a bit cleaner and faster.
You might separate the timestamp extraction and sort steps into two commands and inspect the intermediate result, just to ensure that you have valid strings. Or use the approach below to convert the timestamps, which will catch bad input.

Here's a second way that may add unnecessary work. As an advantage, however, it adds some error checking by converting the ISO8601-formatted timestamp.
If the string you pull out of the header is not a valid timestamp, the first conversion step (ISO8601-to-epoch) will fail and you will then be able to catch invalid input pretty quickly, which can add confidence that your output is correct.
(1) Extract the start times from the headers:
$ awk '(/^>/){ match($0, /start_time=(.*)/, a); print a[1]; }' my.fa > starts_iso8601.txt

(2) Convert the ISO8601-formatted times to time-since-epoch times:
$ while IFS='' read -r line || [[ -n "$line" ]]; do date -d"$line" +%s; done < starts_iso8601.txt > starts_epoch.txt

(3) Sort the epoch times — which are just integers — to another file, in ascending order:
$ sort -n starts_epoch.txt > starts_epoch.sorted.txt

(4) Convert sorted epoch starts back to ISO8601 format:
$ while IFS='' read -r line || [[ -n "$line" ]]; do date -d"@$line" -u +"%Y-%m-%dT%H:%M:%SZ"; done < starts_epoch.sorted.txt > starts_iso8601.sorted.txt

(5) Use that file of ordered starts to write out FASTA records in your sort order:
$ while IFS='' read -r line || [[ -n "$line" ]]; do awk -v line=$line -v RS=">" '{ m = match($0,line); if (m) { printf(">%s", $0); } }' my.fa ; done < starts_iso8601.sorted.txt > my_ordered.fa

Note: Use gawk, gdate, gsort in place of awk, date and sort if you're using OS X, Homebrew, and GNU core utilities.
I'm showing the steps separately to explain the general procedure.
Many or all of these steps can be glued together in a Unix pipeline to save the time, memory and disk space costs of writing intermediate files.
Also, I'm not doing much error checking here. You might look at cases where headers do not contain start_time, or times are not all in ISO8601, or other cases not anticipated here.

Thank you, @Alex Reynolds.  Everything is working fine, however, the last code (5) seem to be writing some sequences twice. I can notice this because the resulting file my_ordered.fa contains more sequences than the original file my.fa . Could you check on that, please?

Can you provide an example FASTA file which is having this problem?

I've just made a small tweak to my fastx-sort.pl script to get this working. It works on both fasta and fastq sequences, and should even work on a random mixture of the two.
This script can now do a dictionary sort by the bits of an ID string that match an arbitrary search pattern, as well as the previous ID and length-based sorting:
fastx-sort.pl -p '<perl regex pattern>' <fastq file> > sorted.fastq

Here's an example with your fasta input, sorting on the nanopore time string:
$ cat seq.fa | grep start_time
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
$ ~/scripts/fastx-sort.pl -p 'start_time=(.*?Z)' seq.fa | grep start_time
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z

Here's a worked example with one of my fastq files for the nanopore time string:
$ cat called_all.fastq | grep start_time | head -n 3
@3a151247-6720-4071-945e-7a92c1856e45 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=9896 ch=246 start_time=2017-11-27T00:58:42Z
@52e439d6-f14a-4af2-adba-132d240fa967 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=8634 ch=511 start_time=2017-11-27T00:58:41Z
@82e3dd50-c395-4ade-ab6b-0eed530d3c3f runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=8312 ch=299 start_time=2017-11-27T00:58:39Z
$ ~/scripts/fastx-sort.pl -p 'start_time=(.*?Z)' called_all.fastq | grep start_time | head -n 10
@195ad363-1d5a-40bf-abed-820c8acce8fc runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=9 ch=144 start_time=2017-11-26T22:42:32Z
@5df2fc6e-53c8-489c-8002-c97f151104bd runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=11 ch=129 start_time=2017-11-26T22:42:37Z
@23940fe6-0eef-418f-9156-d0a26df9bd80 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=22 ch=330 start_time=2017-11-26T22:42:39Z

This script now does a stable sort, so you can do things like sorting by channel, then by read number:
~/scripts/fastx-sort.pl -n -p 'read=([0-9]+)' called_all.fastq | ~/scripts/fastx-sort.pl -n -p 'ch=([0-9]+)' | grep 'ch='
...
@876327ec-ccb0-44ad-a9fb-e008744d6c15 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=5974 ch=66 start_time=2017-11-27T00:04:52Z
@e3ec5d92-ee8a-43d6-82ed-e359e0539b75 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=6117 ch=66 start_time=2017-11-27T00:06:55Z
@410fe169-c8d1-4686-b92f-cbeb0602b2f0 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=6280 ch=66 start_time=2017-11-27T00:09:26Z
@d27c2827-1746-4827-9fec-3075e08d3d64 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=6375 ch=66 start_time=2017-11-27T00:10:44Z
@43079baa-8591-4498-bbfd-5f5ad677c011 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=32 ch=69 start_time=2017-11-26T22:43:32Z
@2bd9b317-438d-4e45-958b-7f77df1dd51a runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=170 ch=69 start_time=2017-11-26T22:47:32Z
@e3599944-849a-492c-ad6d-8cc319b7c5d7 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=207 ch=69 start_time=2017-11-26T22:49:11Z
@8b9d80d0-b28d-46bb-a4e9-6bba95a79ff4 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=320 ch=71 start_time=2017-11-26T22:50:33Z
@b964f164-c36c-4f8a-8a1a-1ad196c996a3 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=780 ch=71 start_time=2017-11-26T22:59:14Z
@92497c10-a46a-4a31-bc18-cfa334491ca9 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=937 ch=71 start_time=2017-11-26T23:02:32Z
@2a715984-6e14-4f31-bd47-d36ee5cbc3e8 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2271 ch=71 start_time=2017-11-26T23:24:11Z
@ff2a59ae-4909-4f3e-bea2-b695d521bf5e runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2441 ch=71 start_time=2017-11-26T23:27:01Z
@f106792d-31a4-44f2-a4fe-83f35a42a573 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2559 ch=71 start_time=2017-11-26T23:28:54Z
@8de1d1ca-c6f7-43f5-bbd6-3c79187c1bbe runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2930 ch=71 start_time=2017-11-26T23:34:33Z
...

Additional notes:

In order to do the sort, the script currently stores the entirety of the input in memory.

Thank you @gringer.  Can you modify your script to use fasta and write an ascending ordered fasta as output?

It works with both fasta and fastq (and even a mixture of the two). The sort order is ascending by default.

My go-to tools for any type of fasta file manipulation are these two simple awk scripts:

FastaToTbl will convert fasta to "tbl" formato: the ID, a tab and then the sequence, all on the same line.
#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
      if (NR>1)
                     printf "\n%s\t", substr($0,2,length($0)-1)
      else
       printf "%s\t", substr($0,2,length($0)-1)
            else
                    printf "%s", $0
}END{printf "\n"}

TblToFasta will convert tbl back to fasta:
#!/usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld

     if (fld == NF-1)
      {
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls)
  {
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}

Armed with those and a bit of Perl, you can do:
$ FastaToTbl file.fa |
    perl -lne '/start_time=(\S+)/;
               $d=`date -d $1 +%s`;
               $k{$d}=$_;
               END{
                    print "$k{$_}" for sort(keys(%k))
               }' | TblToFasta
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTAT
GTGGTATGTGGTATGTGGTATGTGGTATG
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTT
TCCTGTCGCTCTATCTTCAC
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGG
TATGTTGGTATGTTGGTATGTTGGTATG

Explanation
The perl script will take the tbl format sequences (so everything on one line), and identify the start_time: /start_time=(\S+)/. The parentheses capture the matched pattern so the start time will now be $1. Then, we run the system date command giving it the start time as the date and getting the result in seconds since the epoch (%s). note that this assumes GNU date (the default on Linux systems) which can handle the -d option to specify a date. The next step is using $d, the date in seconds, as a key for a hash whose values are the input lines: $k{$d}=$_. Then, after the whole file has been processed, we print each entry of the hash in sorted order (print "$k{$_}" for sort(keys(%k))). Finally, TblToFasta changes back to fasta format.


text:
How to trim adapter sequences from GSE65360 in order to map the reads?
I am trying to map single-cell chromatin accessibility data from doi:10.1038/nature14590, obtained using scATAC-seq, to the reference genome. Example paired-end reads are here.
The reads contain various adapter sequences, which the authors "trimmed from FASTQs using custom python scripts" [ref]. I am having problems reproducing this trimming. They used Nextera Transposase Adapters

TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
  GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG

and custom barcode adapters which are listed [here]. An excerpt from this list:

Custom Barcodes Adapter 1 (index i5):
  AATGATACGGCGACCACCGAGATCTACACTAGATCGCTCGTCGGCAGCGTCAGATGTGTAT
  AATGATACGGCGACCACCGAGATCTACACTAACCAAGTCGTCGGCAGCGTCAGATGTGTAT
  AATGATACGGCGACCACCGAGATCTACACAGGTTGCCTCGTCGGCAGCGTCAGATGTGTAT
  Custom Barcodes Adapter 2 (index i7):
  CAAGCAGAAGACGGCATACGAGATACCCAGCAGTCTCGTGGGCTCGGAGATGTG
  CAAGCAGAAGACGGCATACGAGATCCCAACCTGTCTCGTGGGCTCGGAGATGTG
  CAAGCAGAAGACGGCATACGAGATCACCACACGTCTCGTGGGCTCGGAGATGTG

Analyzing the FASTA-files with FastQC showed that the reverse complement of those sequences was over-represented in the data and also confirmed Nextera Transposase Adapters.
I tried to trim the sequences using cutadapt:
cutadapt --minimum-length=20 -a CACATCTCCGAGCCCACGAGAC -a CTGTCTCTTATA -A ATACACATCTGACGCTGCCGACGA -A CTGTCTCTTATA -o SRR1780164_1_trimmed.fastq.gz -p SRR1780164_2_trimmed.fastq.gz SRR1780164_1.fastq.gz SRR1780164_2.fastq.gz
with

CACATCTCCGAGCCCACGAGAC being the reverse complement of the 2nd fix part of the i7 adapters
ATACACATCTGACGCTGCCGACGA being the reverse complement of the 2nd fix part of the i5 adapters
CTGTCTCTTATA for the Nextera Transposase Adapters

FastQC does not indicate overrepresented sequences in the resulting FASTA-files, however still indicates Nextera Transposase Adapters in the _1 files. On the other hand, after TrimGalore removes CTGTCTCTTATA, Nextra Adapters are not indicated anymore, for some reason.
Furthermore, mapping the reads using bowtie2
bowtie2 -x GRCh38_primary_assembly_index -1 SRR1780164_1_trimmed.fastq.gz -2 SRR1780164_2_trimmed.fastq.gz > SRR1780164.sam
results in a low alignment rate
1352 (2.79%) aligned concordantly exactly 1 time
3420 (7.07%) aligned concordantly >1 times
11.25% overall alignment rate

which gives me the impression I did something wrong during adapter trimming.
Any help would be much appreciated.
--
Edit: As I am not sure about the sequences to use for trimming, I also tried -a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -A CTGTCTCTTATACACATCTGACGCTGCCGACGA, as found in a suggestion for Nextera Adapters, however it did not remove the sequence over-representation in the -1 file. -a CACATCTCCGAGCCCACGAGAC -a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -A CTGTCTCTTATACACATCTGACGCTGCCGACGA solves adapter-issues according to FastQC, but also does not give good mapping results.

First the Nextera adapters and the custom barcode adapters overlap each other
Nextera 1                                     TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
Barcode 1   GATACGGCGACCACCGAGATCTACACTAGATCGCTCGTCGGCAGCGTCAGATGTGTAT

Nextera 2                                    GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
Barcode 2    CAAGCAGAAGACGGCATACGAGATACCCAGCAGTCTCGTGGGCTCGGAGATGTG

Looking at some read 1s and read 2s you were correct that the adapters appear at the 3' end of the read reverse complimented (The |> points to the adapter in the sequence).  You can see that CTGTCTCTTATACACA would be the ideal sequence to use for the 3' end of the reads.  Since the Nextera adapter and the barcode overlap you only need to use the one that is adjacent to the insert.
SRR1780164_2.fastq: GCCCGAGGGGACGCTGCCCCACGTGCTCCAGG |> CTGTCTCTTATACACATCTGACGCTGCCGACGATGCACGAAGTG
SRR1780164_1.fastq: GGGCTGGACAGCGGCGGGGACGCGGACCTGCAG |> CTGTCTCTTATACACATCTCCGAGCCCACGAGACACAAAGTGA

Of the 64k reads there are ~11k with the Nextera adapter at the 3' end based on a rough sequence grep
On the 5' end of the read there is some Nextera adapter present but it is not reverse complimented and there are only a handful of reads with the adapter present ~250 total between read 1 and read 2.
SRR1780164_2.fastq: TGGGCTCGGAGATGTGTATAAGAGACAG <| CTCCCCCGCCTCCAGCATCCGGGCGAGGTAGTGCATCGACGCGTCCAC

I used this command:
cutadapt --minimum-length=20 -A CTGTCTCTTATACACA -a CTGTCTCTTATACACA -G GATGTGTATAAGAGACAG
-g GATGTGTATAAGAGACAG -o SRR1780164_1_trimmed.fastq.gz -p SRR1780164_2_trimmed.fastq.gz
SRR1780164_1.fastq SRR1780164_2.fastq

Which gives a trimming result of
   === Summary ===

Total read pairs processed:             64,832
  Read 1 with adapter:                  16,721 (25.8%)
  Read 2 with adapter:                  16,996 (26.2%)
Pairs that were too short:                  46 (0.1%)
Pairs written (passing filters):        64,786 (99.9%)

I think this best represents the library construct and the trimming of the adapter sequences


text:
Mini-bulk RNA-seq and scRNA-seq analysis differences?
We have RNA-seq data from libraries prepared using the Smartseq2 single-cell protocol on 500 cells (mini-bulk) / library. The complexity is much better than with single-cells (~14k genes for 1.5M reads).
6 cell types / biological replicates were collected with flow cytometry, and there are 4 control and 5 diseased biological replicates (patients).
We are mostly interested in the differences between control and disease either overall or within each cell type.
Can single-cell tools (e.g. Seurat) be used on this dataset? It is neither bulk nor single-cell data. How would you approach this?

I would process it both ways and check if they are comparable. If not then I would need look into the methods more in detail (To see if they are applicable in this data)

Single-cell tools were developed to address some specific factors unique to scRNA-seq, such as:

a lot of dropouts, so hundreds or thousands of detected genes, not tens of thousands (for mammalian genomes)
hundreds or thousands of replicates
not necessarily well-defined groups

Thus, I would say your experiment is a lot closer to bulk RNA-seq than scRNA-seq.


text:
How can I apply proportional (p) distances (Nucleotide) using bioPython
I am working on phylogenetic tree generation. I used BioPython and ClustalW 2.x.x for this purpose. I have generate the tree using BioPython but when I try to generate tree using "MegaSoftware GUI" tool my tree does not matched with the output of MegaSoftware.
When I use "No. of difference" model then results are matched with biopython tree Like below image shown.

But whenever I choose "p-distance" model my results are not matched.

Output MegaSoftware using p-distance

Output BioPython

Here is my code for generating phylogenetic tree:
def conversion_fasta_dnd(clustalw, filepath, type, clustering, output, pwmatrix):
    """
        ClustalW tool to conversion the file format.
        :param clustalw: path of the executable programm.
        :param filepath: path of `fasta` extension file.
        :param type: dna, protein
        :param clustering: upgma, nj
        :param output: clustal, phylip, fasta
        :param pwmatrix: blosum, pam, id
        :return: dnd, aln (files)
        """
    # esblish the full path to the clustalW program
    clustalw_exe = clustalw

    # check to make sure the program is there
    assert os.path.isfile(clustalw_exe), "Clustal W executable missing"

    # create the appropriate command line
    clustalw_cline = ClustalwCommandline(clustalw_exe,
                                         infile=filepath,
                                         type=type,
                                         clustering=clustering,
                                         output=output,
                                         pwmatrix=pwmatrix)

    # print the command line
    print(clustalw_cline)

    # execute the command
    stdout, stderr = clustalw_cline()

# This function draw the tree on image file
def simple_tree(filepath):
    """
    Generate the simple tree image.
    :param filepath: path of dnd file format.
    """
    def get_label(leaf):
        return leaf.name

    file = filepath
    tree = Phylo.read(file, 'newick')
    Phylo.draw(tree, label_func=get_label, do_show=False, branch_labels=lambda c: c.branch_length)

    # axis management
    pylab.axis('on')
    pylab.savefig(OUTPUT_IMAGE, format='png', bbox_inches='tight', dpi=100, bbox_extra_artists=[])

In above code I have used two things first BioPhython lib and the other is Clustal W for align the sequences.
Could you please clarify what p-distances are?
I only have these reference for p-distance:

p-distance (Nucleotide).
Proportional (p) distances

I need to add p-distance model in my code, someone please help me out How can I do that.

Guide me on how to write this P-distance code in a matrix form.
I want to add P-distance matrix.

I think you should make a new question with a reproducible example.

Based on the sites you mention, it seems that p-distances are pairwise distances between aligned sequences obtained by counting the number of differences and dividing by the length of the alignment. This is the proportion of differences between two sequences.
I'm not sure exactly what you mean by "add p-distance model in my code", but here is how you could compute those distances from an existing alignment (here I use as an example an alignment I have in fasta format):
from Bio import AlignIO

def p_dist(seq1, seq2):
    ali_len = len(seq1)
    assert ali_len == len(seq2)
    return sum(0 if nuc1 == nuc2 else 1 for (nuc1, nuc2) in zip(seq1, seq2)) / ali_len

ali = AlignIO.read("most_targeted_histones/transcript_alignments/H4_top4_loss_spliced_NT.fa", "fasta")
nb_seq = len(ali)
distances = {}
for i in range(nb_seq - 1):
    rec1 = ali[i]
    j = i + 1
    while j < nb_seq:
        rec2 = ali[j]
        distances[(rec1.id, rec2.id)] = p_dist(rec1.seq, rec2.seq)
        j += 1

print(*distances.items(), sep="\n")

Output:
(('ZK131.4_spliced', 'ZK131.8_spliced'), 0.0)
(('ZK131.4_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.4_spliced', 'ZK131.1_spliced'), 0.0)
(('ZK131.8_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.8_spliced', 'ZK131.1_spliced'), 0.0)
(('T10C6.14_spliced', 'ZK131.1_spliced'), 0.10897435897435898)

If you want to build a tree from those distances, I'm not sure you can do this with BioPython. It seems however that there are some possibilities with DendroPy.
You would need to write the distances in a suitable format, and to choose the method between UPGMA and neighbour-joining.


text:
Finding data to validate results
My group has a complex datset and I would like to validate the results and check the methods I use in a similar dataset (better if it has been already studied).
Characteristics of my dataset:

It is related to the inflammatory bowel disease (IBD)
More specifically to Crohn's disease.

It has RNA-seq
I have explored the GEO and there are few projects (10) with RNA-seq on the tissue of my disease, and I don't know how to find the 16S-seq of those projects if it exists (see point below).

It has 16S-seq (and not only from the stools)
The ENA seems to have 16S sequencing but most of them are from stools not from biopsies and still, I don't know how to find the RNA-seq of those projects if it exists (see point above)

It has both RNA-seq and 16S-seq from the same individuals at the same time point.
The NIH Human Microbiome Project for my disease has only one individual with that both 16S-seq and RNa-seq of the same patient at the same time.

How can I find a dataset similar to mine?

The snarky answer is that it's better to find a matching dataset at the experimental design phase (i.e. before these studies are carried out), so that the bioinformatics component on the other side of things ends up being a lot easier and more robust. Although it can be quite time consuming to do a database search at the start, it's even harder to find matching projects if the legwork hasn't been done during the design phase. Researchers should be thinking about how hypotheses will be validated before the experiments are carried out.
Every additional characteristic that is chosen for comparison adds complications with regards to matching data, and reduces the applicability of results when there is a mismatch. You would have to be somewhat lucky to be able to find a IBD-Crohn's dataset with paired RNA-seq data and 16s-seq (presumably short-read; is 454 good enough?) from a specific non-stool tissue.
Apart from the places you've already looked, there's dbGaP and SRA. Given the number of characteristics that have been mentioned, I expect you'll have to make some compromises when comparing other datasets.

Thanks for the tip, but is hard to find projects at the experimental design phase. We have an easier dataset in our lab, but I would prefer to confirm it with other datasets. I had forgotten about these databases, thanks.


text:
NCBI - edirect suite to download all genome sequences associated with a query - troubleshooting
Not sure if this is allowed but I cannot think of a better place to ask. I am attempting to download all genomic sequences from refseq associated with the query Peptostreptococcaceae. I have been been using the following command (based on some reading of the manual):
esearch -db nuccore -query "Peptostreptococcaceae" | efilter -molecule genomic -source refseq | efetch -format fasta

But I keep getting the following error
Unknown option: molecule
Unknown option: source
Must supply -query or -days or -mindate and -maxdate arguments on command
line Db value not found in fetch input No do_post output returned from
'http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=nuccore&term=Peptostreptococcaceae&retmax=0&usehistory=y&edirect=3.80&tool=edirect&email=lgunn@sonic15' Result of do_post http request is .............
(a lot of text which I assume is .html or .xml related)

What operating system are you using? I just installed edirect using the instructions you linked to, then ran the exact command you show and it seems to be working fine. Could you be using an older version of edirect that doesn't support these options?

Quite possible. I tried looking up versions but could not find the info. Ill remove it and reinstall it. Fingers crossed.

That was it @terdon. Thanks. Theres no mention of updates or versions so I assumed the one I had would work fine.

Great! Could you please post a quick answer (and then accept it when the system lets you—you need to wait a bit) explaining this? That way, the question can be marked as solved.

The issue here was I was using an older version of eDirect. There is no command to check version or updates. Your best bet is if you experience problems uninstall and reinstall.

If you want to find version info you can either :
Ask for another edirect tool version (it will be the same version for edirect) :
efetch -version

6.80

Check the info contained in the edirect.pl script :
grep -i 'version = ' edirect.pl

$version = "6.80";


text:
Trying to show the gaps of each seq on bio::Graphics after converting clustalw
I want a box representing each sequence, positioned as they are in the alignment and with gaps shown as breaks in the each box.
I've been having trouble for a while with this and have been trying to get this on my own, but I can't. I'm trying to display the gaps using Bio::Graphics after retrieving my sequences from an alignment file.
I removed the exit portion. I was trying to format the code for this website and accidentally put the exit there. The input fill is a custalw file which is an alignment file with sequences on it. The sequences have "-" that represent gaps.
The command is ./aln.pl test-bioaln.aln > ll.png
My expected output would be a .png file with sequences that are aligned and show gaps as a different color.
The input file is:
CLUSTAL W (1.81) multiple sequence alignment

JD1:1:102:1601:ORFJ00027      ------------------------------atgtataaacaacaatattttatttct--c
94a:1:107:117:orf00001        ------------------------------atgtataaacaacaatattttatttct-ac
118a:1:106:158122218:orf00020 ------------------------------atgtataaacaacaatattttatttct-gc
B31:1:100:4091:ORFB0018       ------------------------------atgtataaacaacaatattttatttctggc
72a:1:105:32:orf00022         ------------------------------atgtataaacaacaatattttatttctggc
64b:1:110:473:orf00001        ------------------------------atgtataaacaacaatattttatttctggc
29805:1:108:171:orf00001      ------------------------------atgtataaacaacaatattttatttctggc
BOL26:1:111:60:orf00001       ------------------------------atgtataaacaacaatattttatttctggc
CA-11.2A:1:109:33:orf00001    ------------------------------atgtataaacaacaatattttatttctggc
WI91-23:1:112:493:orf00001    ------------------------------atgtataaacaacaatattttatttctggc
297:1:103:411:ORFB00012       ttggatagattttatacaaagaaggtaataatgtataaacaacaatattttatttctggc
N40:1:101:1716:ORFK00021      ------------------------------atgtataaacaacaatattttatttctggc
ZS7:1:113:22:orf00001         ------------------------------atgtataaacaacaatattttatttctggc
                                                            ******************************

JD1:1:102:1601:ORFJ00027      aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
94a:1:107:117:orf00001        aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
118a:1:106:158122218:orf00020 aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
B31:1:100:4091:ORFB0018       aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
72a:1:105:32:orf00022         aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
64b:1:110:473:orf00001        aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
29805:1:108:171:orf00001      aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
BOL26:1:111:60:orf00001       aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
CA-11.2A:1:109:33:orf00001    aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
WI91-23:1:112:493:orf00001    aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
297:1:103:411:ORFB00012       aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
N40:1:101:1716:ORFK00021      aaggtgcaaggtgttggttttagatttttcacagagcaaatagcaaataatatgaaacta
ZS7:1:113:22:orf00001         aaggtacaaggtgttggttttagattttttacagagcaaatagcaaataatatgaaacta
                              ***** *********************** ******************************

JD1:1:102:1601:ORFJ00027      aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
94a:1:107:117:orf00001        aaaggatttgtaaaaaatctaaacgatggaagggtagaaattgtagctttctttaatact
118a:1:106:158122218:orf00020 aaaggatttgtaaaaaatctaaacgatggaagggtagaaattgtagctttctttaatact
B31:1:100:4091:ORFB0018       aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
72a:1:105:32:orf00022         aaaggatttgtaaaaaatctaaacgatggaagggtagaaattgtagctttctttaatact
64b:1:110:473:orf00001        aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
29805:1:108:171:orf00001      aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
BOL26:1:111:60:orf00001       aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
CA-11.2A:1:109:33:orf00001    aaaggatttgtaaaaaatctaaacgatggaagggtagaaattgtagctttctttaatact
WI91-23:1:112:493:orf00001    aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
297:1:103:411:ORFB00012       aaaggatttgtaaaaaatctaaacgatggaagggtagaaattgtagctttctttaatact
N40:1:101:1716:ORFK00021      aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
ZS7:1:113:22:orf00001         aaaggatttgtaaaaaatctcaacgatggaagggtagaaattgtagctttctttaatact
                              ******************** ***************************************

JD1:1:102:1601:ORFJ00027      aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattaaa
94a:1:107:117:orf00001        aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
118a:1:106:158122218:orf00020 aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
B31:1:100:4091:ORFB0018       aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattaaa
72a:1:105:32:orf00022         aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
64b:1:110:473:orf00001        aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattaaa
29805:1:108:171:orf00001      aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
BOL26:1:111:60:orf00001       aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattaaa
CA-11.2A:1:109:33:orf00001    aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
WI91-23:1:112:493:orf00001    aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
297:1:103:411:ORFB00012       aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
N40:1:101:1716:ORFK00021      aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattgaa
ZS7:1:113:22:orf00001         aaagaacaaatgaaaaaatttgaaaaattattaaatgggaataagtattcaaacattaaa
                              ********************************************************* **

JD1:1:102:1601:ORFJ00027      aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
94a:1:107:117:orf00001        aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
118a:1:106:158122218:orf00020 aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
B31:1:100:4091:ORFB0018       aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
72a:1:105:32:orf00022         aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
64b:1:110:473:orf00001        aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
29805:1:108:171:orf00001      aacattgaaaaaatagctttagatgaaaattatccttttcaatttaatgattttaaaatt
BOL26:1:111:60:orf00001       aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
CA-11.2A:1:109:33:orf00001    aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
WI91-23:1:112:493:orf00001    aacattgaaaaaatagctttagatgaaaattatccttttcaatttaatgattttaaaatt
297:1:103:411:ORFB00012       aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
N40:1:101:1716:ORFK00021      aacattgaaaaaatagctttagatgaaaattatccttttcaatttaatgattttaaaatt
ZS7:1:113:22:orf00001         aacattgaaaaaatagttttagatgaaaattatccttttcaatttaatgattttaaaatt
                              **************** *******************************************

JD1:1:102:1601:ORFJ00027      tattat
94a:1:107:117:orf00001        tattat
118a:1:106:158122218:orf00020 tattat
B31:1:100:4091:ORFB0018       tattat
72a:1:105:32:orf00022         tattat
64b:1:110:473:orf00001        tattat
29805:1:108:171:orf00001      tattat
BOL26:1:111:60:orf00001       tattat
CA-11.2A:1:109:33:orf00001    tattat
WI91-23:1:112:493:orf00001    tattat
297:1:103:411:ORFB00012       tattat
N40:1:101:1716:ORFK00021      tattat
ZS7:1:113:22:orf00001         tattat
                              ******

Here's the code:
#!/usr/bin/perl
use Bio::AlignIO;
use Bio::Graphics::Panel;

my $line = shift @ARGV;
my $in = Bio::AlignIO->new(-file=>$line,-format=>"clustalw");
while($aln = $in->next_aln()){
    foreach $seqobj($aln->each_seq()){
        my $seq = $seqobj->seq;
        my $id = $seqobj->id;
        my $length = $seqobj->length;
        my $seqobj = Bio::SeqFeature::Generic->new(-start =>1, -end=>$length,-display_name=>$id);
        push (@seq, $seqobj);
    }
    foreach $seq(@seq){
        my @features = $seq->get_SeqFeatures;
        my %sorted_features;
        for my $f (@features) {
            my $tag = "\-";
            push @{$sorted_features{$tag}},$f;
        }

        my $panel = Bio::Graphics::Panel->new(
            -length    => $seq->length,
            -key_style => 'between',
            -width     => 800,
            -pad_left  => 10,
            -pad_right => 10,
        );
        $panel->add_track(generic => Bio::SeqFeature::Generic->new(-start=>1,
                -end=>$seq->length),
            -glyph  => 'generic',
            -bgcolor => 'blue',
            -label  => 1,
        );

        my @colors = qw(cyan orange blue purple green chartreuse magenta yellow aqua);
        my $idx    = 0;
        for my $tag (sort keys my %sorted_features) {
            my $features = $sorted_features{$tag};
            $panel->add_track($features,
                -glyph    =>  'generic',
                -bgcolor  =>  $colors[$idx++ % @colors],
                -fgcolor  => 'black',
                -font2color => 'red',
                -key      => "${tag}s",
                -bump     => +1,
                -height   => 8,
                -label    => 1,
                -description => 1,
            );
        }
        print $panel->png;

    }
}

Unfortunately, $seq->get_SeqFeatures never returns anything. I'm not sure why it doesn't since there should be features on it coming from the Bio::SimpleAlign object. Maybe I should set some new parameters?

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


text:
Link Prediction in Bipartite Networks between biological pathways and drugs
Currently I'm building a big matrix (using microarray, mass spectrometry, RNAseq data) that consist by pathways (rows) and treatment/drugs (columns).
The values of that matrix are scores that describe the effect/relation of that specific treatment/drug on each pathway.
Now, some of the drugs is known to target specific disease (target disease) and thus is expected to have a high score with that disease's pathways. (A)
The rest of the drugs is known to target other diseases.(B)
My thought now is that, if I label the known drugs as "1" (or whatever), could I use a prediction model to find out new relationships (links) between drugs from the B class and the significant disease pathway from A ?
The whole thing, we can represent it as a bipartite network shown below.

For example, here the denoted as green, drugs are the known drugs belong to the class A that mentioned before and they are connected highly with some pathways of the target disease.
From the other hand blue and red ones are drugs taken from other cases (class B) that are used to target other diseases, but as my plot shows, blue seems to has an "affinity" to some of the pathways of the target disease and thus should be an alternative candidate for that disease. The red from the other side, doesn't have any "affinity" to any significant pathway and thus we can not consider it as possible candidate.
To find a solution to that kind of problem, I looked for algorithms and methods developed for social networks (e.g recommendation etc.) and mostly I read about random walks with return.
But before to proceed, I would like to know if someone else tried such an approach here and having more experience.

I am not sure if I understood correctly how you classify  the drugs. But what you attempt to do is similar to what they do here but in that article they use the drugs targeted to a molecule/pathway to find combination of drugs that are better for the disease (not an alternative drug), see the first image:

In that article they mention the pathway cross-talk inhibition (PCI), you might want to calculate it for your pathways to assess the relative reduction of the pathway efficiency.


text:
How to split multifasta based on partial fasta header
I have download a batch of refseq fasta files and want to split them based on strain. this is complicated by accession numbers, project IDs etc. I downloaded the genomes via eDirect from NCBI. I know I can split the output but I want them together for some preliminary parsing (remove ambiguous and plasmid sequences etc.)
See example below.
>NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
>NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
>NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
>NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC

I want two files:
Clostridium botulinum BKT028387.fas and Peptoanaerobacter stomatis OBRC8.fas
My actual dataset is much larger (~70,000 genomes) and not too sure about the headers in each case.

Im thinking the following could work but Im not sure how to execute it:
Using a grep command to extract the lines but split the output.
If I download the summary list from refseq and use the names as the query example:
grep -Fwf refseq.list sequences.fa

Summary file

Peptoanaerobacter stomatis strain OBRC8, whole genome shotgun sequencing project 7,810,226 rc other DNA  NZ_AJXR00000000.1 GI:938474759

Clostridium botulinum BKT028387, whole genome shotgun sequencing project 4,134,593 rc other DNA  NZ_NOJY00000000.1 GI:1231817869

But I could edit it down to
Peptoanaerobacter stomatis strain OBRC8
Clostridium botulinum BKT028387

Am I one the right track. How do I finish this?

What operating system are you using? I asked you in your previous question as well but you didn't reply. The tools available to you and their behavior change depending on your OS.

Unix? Do you mean an actual Unix or macOS or Linux? Please tell us your actual OS, a unix-like system can be assumed if you're doing bioinformatics usually.

Yes Its is a linux system .

Great, then you should have the GNU tools. Good. Have a look at the updated answer. What does the summary list look like?

Ah! That summary file can simplify things enormously! Can you link to the entire file? Does it have all the strains and does it have each strain in exactly the same format as it appears in the fasta file header? I am afraid it won't since I've already found inconsistencies. It might be simpler to start from the beginning and download each strain separately.

For the files you have terdon, use the link and use 'send to: > complete record > file > summary > create file` https://www.ncbi.nlm.nih.gov/nuccore/?term=peptostreptococcaceae Im pretty confident they match any left behind I can curate manually.

OK, no, that summary is useless. It is unstructured and non-standard. For example, it has Peptostreptococcaceae bacterium and Uncultured Peptostreptococcaceae bacterium. I'm afraid this isn't the right way of doing what you need. It would probably be best to use a more sophisticated approach and save each sequence to the right file while downloading by also downloading the meta data for the sequence entry.

The SeqBuddy --pull_records command can grab sequences based on a regular expression against record IDs and/or metadata. Given your example:
$: seqbuddy.py sequences.fa --pull_records "Clostridium botulinum.*BKT028387" full > Clostridium_botulinum_BKT028387.fas
$: seqbuddy.py sequences.fa --pull_records "Peptoanaerobacter stomatis.*OBRC8" full > Peptoanaerobacter_stomatis_OBRC8.fas

This looks great. Is there a way to set this up to loop using an input list of names?

Yep, you just pass in the file name.

I had a similar problem and found a solution using awk. Since this is the first time I use awk, the syntax might not be perfect, but you should get an idea of what is possible.
awk -F' ' '{if(/^>/){sub("strain ", ""); name=$2" "$3" "$4; print > name".fa"}else{print > name".fa"}}' yourfile.fa

(more information below)
This should transform a file like this:
yourfile.fa
>NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
ATGCATGC
ATGCATGC
>NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
ATGCATGC
>NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
ATGCATGC
ATGCATGC
>NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC
ATGCATGC

Into two files like this:
'Clostridium botulinum BKT028387.fa'
>NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
ATGCATGC
ATGCATGC
>NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
ATGCATGC

'Peptoanaerobacter stomatis OBRC8.fa'
>NZ_ALNK01000035.1 Peptoanaerobacter stomatis OBRC8 ctg120007290634 ATGC
ATGCATGC
ATGCATGC
>NZ_ALNK01000031.1 Peptoanaerobacter stomatis OBRC8 ctg120007290635 ATGC
ATGCATGC

How it works

-F' ' tells awk to split the line by space (accessible with $1, $2, ...)
if(/^>/) looks for line starting with the > character, then sub("strain ", "") remove the word strain and name=$2" "$3" "$4; print > name".fa" print the line into a new file with a name based on what you asked for.
else if the line does not start with >, {print > name".fa"} appends it to a file using the current value of name.

This solution will work regardless of the number of lines between the fasta headers.
Note that I won't recommend using spaces in filenames, so I would replace name=$2" "$3" "$4 with something like name=$2"_"$3"_"$4.
I used this solution to split a large fasta file (>5 GB) containing >5k genomes with >25k scaffolds and I found it fast and rather memory efficient.

Thank you so much for the code and explanation. However, when I tried it, it gave me the following error. Any ideas?
awk: syntax error at source line 1  context is  {if(/^>/){ name=$2" "$3" "$4; print > >>>  name".fa" <<<  awk: illegal statement at source line 1 awk: syntax error at source line 1

I'm not sure about your error but you can probably put the code on several lines for easier debugging, and execute it with the awk -f option like this: awk -f myscript.awk mygenomes.fasta.

I would convert the fasta file to a format with one line per record (I have already posted the FastaToTbl and TblToFasta scripts I use for this before), and then look for a specific pattern:
FastaToTbl file.fa | grep 'BKT028387' | TblToFasta > clostridium.fa
FastaToTbl file.fa | grep 'OBRC8' | TblToFasta > peptoanaerobacter.fa

I have downloaded a few sequences using the same command you described in your previous question, and it looks like the headers are not standardized. For example, while you have many that are in the format of $species strain $strain, like these:
Clostridioides difficile strain 7496-D/ST9
Clostridioides difficile strain 5616-DH/ST42,
Peptoanaerobacter stomatis strain OBRC8
Romboutsia maritimum strain CCRI-22766
Tepidibacter mesophilus strain JCM

Others don't have the word strain and just give the strain:
Peptostreptococcus anaerobius 653-L contig00088,
Peptoclostridium acidaminophilum DSM 3953

So, as a first step, you could deal with the ones that do have the word strain as the 4th word of the header line (the 1st is the seqID which I have not included in the examples above) and the actual strain as the 5th and get those out of the way:
FastaToTbl foo.fa |
    awk '($4=="strain"){gsub("/","_"); print $0 >> $2"_"$3"_"$5".tbl"}'

That will look for lines whose 4th field is strain, then it will replace any / in the line with _ (there's a strain called CF/ST37 and that cannot be in a filename) and print into a file whose name is the 2nd, 3rd and 5th fields of the line concatenated. On my test file, this produced the following files (as you can see, some are wrong since there is a comma after the strain, but that's because they're not standardized):
$ ls *tbl
Clostridioides_difficile_5479-NonSp_ST229,.tbl  Clostridioides_difficile_6589-DH_ST42.tbl         Clostridioides_difficile_7457-NonSp_novelST,.tbl
Clostridioides_difficile_5479-NonSp_ST229.tbl   Clostridioides_difficile_6592-G_ST8,.tbl          Clostridioides_difficile_7457-NonSp_novelST.tbl
Clostridioides_difficile_5480-E_ST6,.tbl        Clostridioides_difficile_6592-G_ST8.tbl           Clostridioides_difficile_7458-N_ST10,.tbl
Clostridioides_difficile_5480-E_ST6.tbl         Clostridioides_difficile_6593-Y_ST110,.tbl        Clostridioides_difficile_7458-N_ST10.tbl
Clostridioides_difficile_5482-L_ST55,.tbl       Clostridioides_difficile_6593-Y_ST110.tbl         Clostridioides_difficile_7459-G_ST8,.tbl
Clostridioides_difficile_5482-L_ST55.tbl        Clostridioides_difficile_6595-AL_ST58,.tbl        Clostridioides_difficile_7459-G_ST8.tbl
Clostridioides_difficile_5485-Y_ST2,.tbl        Clostridioides_difficile_6595-AL_ST58.tbl         Clostridioides_difficile_7462-DH_ST42,.tbl
Clostridioides_difficile_5485-Y_ST2.tbl         Clostridioides_difficile_6596-D_ST9,.tbl          Clostridioides_difficile_7462-DH_ST42.tbl
Clostridioides_difficile_5488-DH_ST42,.tbl      Clostridioides_difficile_6596-D_ST9.tbl           Clostridioides_difficile_7463-E_ST6,.tbl
Clostridioides_difficile_5488-DH_ST42.tbl       Clostridioides_difficile_6599-NonSp_ST103,.tbl    Clostridioides_difficile_7463-E_ST6.tbl
Clostridioides_difficile_5489-M_ST15,.tbl       Clostridioides_difficile_6599-NonSp_ST103.tbl     Clostridioides_difficile_7464-Y_ST110,.tbl
Clostridioides_difficile_5489-M_ST15.tbl        Clostridioides_difficile_6600-G_ST8,.tbl          Clostridioides_difficile_7464-Y_ST110.tbl
Clostridioides_difficile_5490-CF_ST37,.tbl      Clostridioides_difficile_6600-G_ST8.tbl           Clostridioides_difficile_7466-DH_ST42,.tbl
Clostridioides_difficile_5490-CF_ST37.tbl       Clostridioides_difficile_6606-BK_ST11,.tbl        Clostridioides_difficile_7466-DH_ST42.tbl
Clostridioides_difficile_5493-DH_ST42,.tbl      Clostridioides_difficile_6606-BK_ST11.tbl         Clostridioides_difficile_7467-DH_ST42,.tbl
Clostridioides_difficile_5493-DH_ST42.tbl       Clostridioides_difficile_6607-DH_ST42,.tbl        Clostridioides_difficile_7467-DH_ST42.tbl
Clostridioides_difficile_5495-AF_ST41,.tbl      Clostridioides_difficile_6607-DH_ST42.tbl         Clostridioides_difficile_7468-NonSp_ST97,.tbl
Clostridioides_difficile_5495-AF_ST41.tbl       Clostridioides_difficile_6610-N_ST10,.tbl         Clostridioides_difficile_7468-NonSp_ST97.tbl
Clostridioides_difficile_5497-Y_ST2,.tbl        Clostridioides_difficile_6610-N_ST10.tbl          Clostridioides_difficile_7471-BI_ST1,.tbl
Clostridioides_difficile_5497-Y_ST2.tbl         Clostridioides_difficile_6612-DH_ST28,.tbl        Clostridioides_difficile_7471-BI_ST1.tbl
Clostridioides_difficile_5500-DH_ST42,.tbl      Clostridioides_difficile_6612-DH_ST28.tbl         Clostridioides_difficile_7472-Y_ST2,.tbl
Clostridioides_difficile_5500-DH_ST42.tbl       Clostridioides_difficile_6613-NonSp_ST3,.tbl      Clostridioides_difficile_7472-Y_ST2.tbl
Clostridioides_difficile_5504-L_ST55,.tbl       Clostridioides_difficile_6613-NonSp_ST3.tbl       Clostridioides_difficile_7473-DH_ST42,.tbl
Clostridioides_difficile_5504-L_ST55.tbl        Clostridioides_difficile_6614-AL_ST34,.tbl        Clostridioides_difficile_7473-DH_ST42.tbl
Clostridioides_difficile_5505-Y_ST2,.tbl        Clostridioides_difficile_6614-AL_ST34.tbl         Clostridioides_difficile_7476-Y_ST2,.tbl
Clostridioides_difficile_5505-Y_ST2.tbl         Clostridioides_difficile_6615-NonSp_ST11,.tbl     Clostridioides_difficile_7476-Y_ST2.tbl
Clostridioides_difficile_5537-D_ST9,.tbl        Clostridioides_difficile_6615-NonSp_ST11.tbl      Clostridioides_difficile_7478-NonSp_novelST,.tbl
Clostridioides_difficile_5537-D_ST9.tbl         Clostridioides_difficile_6616-NonSp_novelST,.tbl  Clostridioides_difficile_7478-NonSp_novelST.tbl
Clostridioides_difficile_5538-BK_ST11,.tbl      Clostridioides_difficile_6616-NonSp_novelST.tbl   Clostridioides_difficile_7479-NonSp_ST3,.tbl
Clostridioides_difficile_5538-BK_ST11.tbl       Clostridioides_difficile_6617-A_ST43,.tbl         Clostridioides_difficile_7479-NonSp_ST3.tbl
Clostridioides_difficile_5539-Y_ST2,.tbl        Clostridioides_difficile_6617-A_ST43.tbl          Clostridioides_difficile_7480-M_ST15,.tbl
Clostridioides_difficile_5539-Y_ST2.tbl         Clostridioides_difficile_6618-NonSp_ST53,.tbl     Clostridioides_difficile_7480-M_ST15.tbl
Clostridioides_difficile_5541-CF_ST37,.tbl      Clostridioides_difficile_6618-NonSp_ST53.tbl      Clostridioides_difficile_7481-NonSp_ST59,.tbl
Clostridioides_difficile_5541-CF_ST37.tbl       Clostridioides_difficile_6619-CF_ST37,.tbl        Clostridioides_difficile_7481-NonSp_ST59.tbl
Clostridioides_difficile_5542-Y_ST2,.tbl        Clostridioides_difficile_6619-CF_ST37.tbl         Clostridioides_difficile_7482-DH_ST42,.tbl
Clostridioides_difficile_5542-Y_ST2.tbl         Clostridioides_difficile_6621-NonSp_ST10,.tbl     Clostridioides_difficile_7482-DH_ST42.tbl
Clostridioides_difficile_5547-DH_ST42,.tbl      Clostridioides_difficile_6621-NonSp_ST10.tbl      Clostridioides_difficile_7486-A_ST43,.tbl
Clostridioides_difficile_5547-DH_ST42.tbl       Clostridioides_difficile_6622-CF_ST37,.tbl        Clostridioides_difficile_7486-A_ST43.tbl
Clostridioides_difficile_5548-Y_ST2,.tbl        Clostridioides_difficile_6622-CF_ST37.tbl         Clostridioides_difficile_7488-NonSp_ST103,.tbl
Clostridioides_difficile_5548-Y_ST2.tbl         Clostridioides_difficile_6623-NonSp_ST54,.tbl     Clostridioides_difficile_7488-NonSp_ST103.tbl
Clostridioides_difficile_5549-NonSp_ST21,.tbl   Clostridioides_difficile_6623-NonSp_ST54.tbl      Clostridioides_difficile_7491-DH_ST42,.tbl
Clostridioides_difficile_5549-NonSp_ST21.tbl    Clostridioides_difficile_6625-Y_ST2,.tbl          Clostridioides_difficile_7491-DH_ST42.tbl
Clostridioides_difficile_5550-AH_ST67,.tbl      Clostridioides_difficile_6625-Y_ST2.tbl           Clostridioides_difficile_7492-CF_ST37,.tbl
Clostridioides_difficile_5550-AH_ST67.tbl       Clostridioides_difficile_6626-Y_ST110,.tbl        Clostridioides_difficile_7492-CF_ST37.tbl
Clostridioides_difficile_5551-CF_ST37,.tbl      Clostridioides_difficile_6626-Y_ST110.tbl         Clostridioides_difficile_7494-E_ST6,.tbl
Clostridioides_difficile_5551-CF_ST37.tbl       Clostridioides_difficile_6628-AL_ST58,.tbl        Clostridioides_difficile_7494-E_ST6.tbl
Clostridioides_difficile_5552-DH_ST42,.tbl      Clostridioides_difficile_6628-AL_ST58.tbl         Clostridioides_difficile_7496-D_ST9,.tbl
Clostridioides_difficile_5552-DH_ST42.tbl       Clostridioides_difficile_6629-NonSp_novelST,.tbl  Clostridioides_difficile_7496-D_ST9.tbl
Clostridioides_difficile_5553-L_ST55,.tbl       Clostridioides_difficile_6629-NonSp_novelST.tbl   Clostridioides_difficile_7497-NonSp_novelST,.tbl
Clostridioides_difficile_5553-L_ST55.tbl        Clostridioides_difficile_6631-DH_ST42,.tbl        Clostridioides_difficile_7497-NonSp_novelST.tbl
Clostridioides_difficile_5555-DH_ST42,.tbl      Clostridioides_difficile_6631-DH_ST42.tbl         Clostridioides_difficile_7498-M_ST15,.tbl
Clostridioides_difficile_5555-DH_ST42.tbl       Clostridioides_difficile_6632-NonSp_ST3,.tbl      Clostridioides_difficile_7498-M_ST15.tbl
Clostridioides_difficile_5557-DH_ST42,.tbl      Clostridioides_difficile_6632-NonSp_ST3.tbl       Clostridioides_difficile_7499-CF_ST37,.tbl
Clostridioides_difficile_5557-DH_ST42.tbl       Clostridioides_difficile_6634-DH_ST42,.tbl        Clostridioides_difficile_7499-CF_ST37.tbl
Clostridioides_difficile_5558-AH_ST67,.tbl      Clostridioides_difficile_6634-DH_ST42.tbl         Clostridioides_difficile_7501-BI_ST1,.tbl
Clostridioides_difficile_5558-AH_ST67.tbl       Clostridioides_difficile_6635-BI_ST1,.tbl         Clostridioides_difficile_7501-BI_ST1.tbl
Clostridioides_difficile_5559-DH_ST42,.tbl      Clostridioides_difficile_6635-BI_ST1.tbl          Clostridioides_difficile_7503-Y_ST2,.tbl
Clostridioides_difficile_5559-DH_ST42.tbl       Clostridioides_difficile_6637-AL_ST58,.tbl        Clostridioides_difficile_7503-Y_ST2.tbl
Clostridioides_difficile_5573-CF_ST37,.tbl      Clostridioides_difficile_6637-AL_ST58.tbl         Clostridioides_difficile_7505-Y_ST2,.tbl
Clostridioides_difficile_5573-CF_ST37.tbl       Clostridioides_difficile_7450-DH_ST42,.tbl        Clostridioides_difficile_7505-Y_ST2.tbl
Clostridioides_difficile_5616-DH_ST42,.tbl      Clostridioides_difficile_7450-DH_ST42.tbl         Clostridioides_difficile_NCKUH-21,.tbl
Clostridioides_difficile_5616-DH_ST42.tbl       Clostridioides_difficile_7451-NonSp_ST2,.tbl      Clostridioides_difficile_RC10,.tbl
Clostridioides_difficile_5754,.tbl              Clostridioides_difficile_7451-NonSp_ST2.tbl       Clostridioides_difficile_RC10.tbl
Clostridioides_difficile_5754.tbl               Clostridioides_difficile_7452-Y_ST110,.tbl        Clostridioides_difficile_RF17,.tbl
Clostridioides_difficile_6583-AH_ST67,.tbl      Clostridioides_difficile_7452-Y_ST110.tbl         Clostridioides_difficile_RF17.tbl
Clostridioides_difficile_6583-AH_ST67.tbl       Clostridioides_difficile_7453-Y_ST2,.tbl          Peptoanaerobacter_stomatis_OBRC8,.tbl
Clostridioides_difficile_6584-AQ_ST26,.tbl      Clostridioides_difficile_7453-Y_ST2.tbl           Peptoanaerobacter_stomatis_OBRC8.tbl
Clostridioides_difficile_6584-AQ_ST26.tbl       Clostridioides_difficile_7455-A_ST43,.tbl         Romboutsia_maritimum_CCRI-22766,.tbl
Clostridioides_difficile_6588-G_ST8,.tbl        Clostridioides_difficile_7455-A_ST43.tbl          Romboutsia_maritimum_CCRI-22766.tbl
Clostridioides_difficile_6588-G_ST8.tbl         Clostridioides_difficile_7456-DH_ST42,.tbl        Tepidibacter_mesophilus_JCM.tbl
Clostridioides_difficile_6589-DH_ST42,.tbl      Clostridioides_difficile_7456-DH_ST42.tbl

Finally, you can convert these back to fasta with:
for f in *tbl; do TblToFasta $f > ${f/.tbl/.fa} && rm $f; done

This will only get those sequences whose naming scheme is as described above. But that's a start. You can then see what else you have and extract the rest manually or by finding a similar approach that matches those.

@AudileF we can't know that unless we know the details of your sequences. it all depends on how the strains are defined. I assume this is about the sequences you downloaded using edirect from your other question but that is the type of thing you should mention in your question. It's only luck that I happen to have read the other one. I am trying to give you a more general solution but the names of the few sequences I downloaded are not standard so it will not be easy.


text:
State of the art mutation simulation software
There are many features affecting mutation probabilities, e.g. CpG mutations are 10-fold more likely than other types of mutations.
Is there a model (preferably with software) which can take two aligned genomic regions, estimate parameters of the neutral mutation procces and then simulate mutations for another region? Obviously, there are models like HKY and GTR, but they are somewhat simplistic and do not take the context into account.
One model I found is this, but it seems that there is no software implementing the model. Is there something more recent?
I am also looking for software which can account for indels.

How good the simulation must be? Should it include that in transcription factors and translated genes are less mutated aside from CpG and indels transverstions...? (I am not aware of any software, but I think that with a single reference any tool will have problems estimating good parameters)

@Llopis the idea is to use it on all the regions with very weak selection (e.g. intergenic non-functional), only to account for mutations process, not selection, which of course affects genes and regulatory elements.

I'm not sure these tools allow parameter estimation from the data, but I've heard good things about Indelible https://academic.oup.com/mbe/article/26/8/1879/980884 and Evolver https://www.drive5.com/evolver/EvolverUserGuide.pdf

Does this help:
https://www.biorxiv.org/content/early/2017/11/22/223297 New synthetic-diploid benchmark for accurate variant calling evaluation

FoldX is the state of the art method for mutation modelling (cited by countless papers).  It does not do exactly what you want without extra work though.

The MMOS company developing Citizen science solutions for these kind of tasks.

Question I suspect this question is about the phylogenetics of methylation and the approach the investigator is proposing would be the last approach to use.
Summary
The approaches to assess the phylogenetics of methylation in order of preference are:

dN/dS between CpG sites and non-CpG sites,
An explicit molecular clock between CpG sites and non-CpG sites
A likelihood ration test (LRT) based on a null distribution of mutations generated by a Monte Carlo algorithm

Background/Rationale
The approaches you are using have been around for a long while and generate the random distribution of mutations for a given phylogeny. Its phylogenetics in reverse and the technique is known as Monte Carlo, you start with a randomised probability and send it through a parameterised model to predict the amino acid/nucleotide. Thus it is ML(max likelihood), Bayesian, "phylogenetic-HMM" in reverse. It is used within likelihood ratio test calculations. The model is determined by a standard ML, Bayesian phylogeny algorithm, i.e. its circular mathematics because there is no independent means of calculating the mutation behaviour, so the precise context you are using this for needs to be carefully considered.
Packages There are lots of packages that do this type of Monte Carlo simulation, SOWHat being a good one (LRT) and you generate a large number of replicate data sets (100 or 1000). The one of the base algorithms is "seq-gen", although PAML may have implemented this.
Considerations To use this style of approach you need to carefully consider your question. The "mutation simulations" when analysed via a phylogenetic program will produce the same parameters that you initially set and trees of basically the same length. If you are using this to generate a null distribution for a phylogenetics test these approaches are useful, you then compare the observed likelihood against the null distribution. If you are using it to work out whether the mutation rate at CpG sites is higher than at other sites, its one approach amongst a number of alternatives.
Drawbacks Calculating a de novo null distribution of the mutation rate is computationally very expensive and therefore tends to the last calculation performed. It will address a single hypothesis, usually regarding topology.
Crux of the problem To really do what (I think) you want to do you need an independent measure of mutation and mutational behaviour and its not trivial to achieve. You would have to consider it like a machine learning calculation, with a formal training set, I don't know whether such an approach has been implemented.
Strengths/Summary In summary, Monte Carlo phylogenetics simulation has a limited application because of circularity of the calculation, BUT, BUT, BUT what I have omitted is when it is appropriate it is a very powerful test indeed.

Points of interest I do like the population dynamic simulations within Beast, this certainly "a priori" and its potential IMO needs further exploration. However, I don't think you looking at how molecular epidemiology might impact the mutational behaviour methylation. I have not looked at Indelible, but looks interesting.
I was going to talk about "processor war simulations" but its probably off-topic.

Thanks for your response. My question is not so much about phylogenetics of methylation. What I need essentially is a null distribution for neutral mutations. MCMC or maximum likelihood helps a lot given that there is a model :) But I am not aware of a non-time revesible model accounting for sequence context and indels.

I don't know what "sequence context" is but non-reversibility of mutations (and by inference time) is Beast, all other ML or Bayes (MCMCMC) algorithms have a symmetrical mutation matrices. MCMC in Bayes is different from the Monte Carlo used for nucleotide simulation. To honest the lack of clarity in the question is of concern, so I can't help any further.

In classical phylogenetic models a probability of mutation at position N depends only on the nucleotide at that position and some measure of time. What I mean by sequence context here is the fact that probability of mutation will also depend on nucleotides and mutations at positions N-M..N+M. We have a wealth of empirical evidence that there is such a dependency. But not many models accounting for it (such as this one).
I will try to clarify the question and to give some background information there.

So do I, but I bet you can't name the mechanism.  Anyway Nick Goldman did this years ago, its not new


text:
R package equivalent to RSeQC infer_experiment to get strandedness of RNA-Seq
I am currently writing an R package that includes a module to run featureCounts (gene quantification tool) from Rsubread. I wanted to be able to specify the correct strandedness option to featureCounts without the user needing to specify whether their protocol is unstranded, forwardly stranded, or reversely stranded. Typically when I don't know the protocol for a given BAM file, I use RSeQC's infer_experiment module, which uses the BAM plus an annotation (GTF) file to give an output like so:
Unstranded:
This is PairEnd Data
Fraction of reads failed to determine: 0.0172
Fraction of reads explained by "1++,1--,2+-,2-+": 0.4903
Fraction of reads explained by "1+-,1-+,2++,2--": 0.4925

Reversely stranded (R2 same strand as parent gene):
This is PairEnd Data
Fraction of reads failed to determine: 0.0072
Fraction of reads explained by "1++,1--,2+-,2-+": 0.0487
Fraction of reads explained by "1+-,1-+,2++,2--": 0.9441

Does anyone know of an R package that does something similar?
If I really needed to I could always require that users have RSeQC installed and just run infer_experiment as a system call. But it would really be better to run it from an existing R package if available.
I know I could also write something up myself (use rsamtools to downsample the BAM, check ratio of SAM flags 99/147 vs. 83/163 in regions matching positive and negative strand genes). But would rather not reinvent the wheel here if there is something already available to do what I need.


text:
Folded Protein Chunk Dimensional Classification?
Are there known dimensional measurements for the classification of folded proteins given a starting chunk/domain as defined by something like the clustering functionality of MSM Builder?
Examples of what I would be looking for would be dimensions such as:

pH

Salinity

Dissolved Sugars

Reducing Agent Concentration

Ambient Electronegativity

With corresponding scales and precisions such as:

logorithmic vs linear

step size

distribution (e.g. dissolved/flat vs fractal/inclusions)

The above parameters would be included in folding software after a known starting conformation provided by MSM Builder or similar chunking software to derive a GUID for proteins without (or with extremely minimal and tuneable) collisions?
For instance, if hundreds of GUIDs on average end up equating to the same folded protein that's fine, but if one GUID could map to two or more thermodynamically stable folded proteins at greater than some incredibly small (e.g. 0.0000001%) then it wouldn't fit what I'm looking for.  I know this is inherently NP-hard, but the complete ID would be something like:
{
    sequence,
    chunkId,
    pH,
    salinity,
    dissolvedSugar,
    reducingAgentConcentration
}

So I would expect the chunkId and dimensional parameters to play a significant role in reducing collisions for different starting parameters.
It's worth noting I'm not seeking some magic bullet where you type an ID and get a fully folded protein, just a reliable system of GUIDs with enough starting information embedded in the ID to start the folding process.

CATH
The CATH database classifies protein by fold: https://www.cathdb.info/
So the value from that is probably the most useful for you.
Crystallographic conditions

pH
Salinity
Dissolved Sugars
Reducing Agent Concentration

Your crystallising conditions do not mean much. They are a solution that is close to precipitating your protein but slowly enough for it to crystallise. A protein can crystallise under multiple different conditions and a crystal may be chosen simply because it is a lower index well in the Hampton screening conditions...
These mean very very little. And unless you are planning on doing ML to determine the best screening condition these should be avoided.
Data from the PDB
You have probably noticed that the only stats PDB/PDBe give are crystallographic, like space group, unit cell dimensions, resolution, completeness in shells, Rfree, Rmerge etc. So nothing there is of much use... except for bound ligand. Often this is the cofactor, but often it's a cofactor analog or a product (as opposed to the substrate). But that information can come from Uniprot etc.
pH
A solution has a pH, a protein has a pI. At a pH equal to the pI it will be neutral —and less insoluble. ProtParam can calculate the pI from a linear sequence. However, surface charge totally shows different spots. APBS can map the solution to the Poisson-Boltzmann equation to your protein. There is a PyMOL wrapper, a Pyrosetta wrapper and more. But you still would need to some stats and 3D operations to convert it into a small vector of values —rototranslations and projects are not for the faint of heart or in a rush.
Melting temperature
BRENDA has lots of enzymatic data on enzymes including melting temperature or kinetics at different temperatures (which follow the Eyring equation (TST), until MMRT kicks in) which could be useful.
For a given protein you can calculate the relative folding ∆∆G (eg. with Rosetta Score), which is a proxy to melting temperature in a buffer it is happy in.
Despite paining me to say... some fanciful Poisson-Boltzmann map derived vector and ∆∆G are probably all you could get from analysing structures. Whereas the rest are from CATH classification and other databases that are not structure based.


text:
Correlation network representation
I have troubles understanding the figure 10 of this paper:

  Similarity network of participating methods for BPO. Similarities are computed as Pearson’s correlation coefficient between methods, with a 0.75 cutoff for illustration purposes. A unique color is assigned to all methods submitted under the same principal investigator. Not evaluated (organizers’) methods are shown in triangles, while benchmark methods (Naïve and BLAST) are shown in squares. The top-ten methods are highlighted with enlarged nodes and circled in red. The edge width indicates the strength of similarity. Nodes are labeled with the name of the methods followed by “-team(model)” if multiple teams/models were submitted

This study calculated the pairwise Pearson correlation between methods on a common set of gene-concept pairs and then visualized these similarities as networks. How is the place of the methods determined in this visualization?
For example, the relationship between Pannzer (1) and Pannzer (2) is greater than the relationship between Pannzer (1) and Pannzer (3), but the Panzer (1) and Panzer (3) are closer to each other than Panzer (1) and Panzer (2). Is this proximity positioning random?  or is there a certain rule?

Where is this image from? If it's from a published article, please include the reference both so we can look it up and so it is properly cited, and so we can get a better quality image.

This image belongs to article that an expanded evaluation of protein function prediction methods shows an improvement in accuracy

Thanks, but we needed something like the edit Llopis did. Whenever you show something from a paper you must always include the source.

I did not know that I needed to add the source. Thank you for saying that to me

When plotting networks where only the weight of an edge is known like this one. The relative position between elements of the networks is irrelevant. The only informative part of the graph is (usually) the size of the edges, where is the thicker means more weight. This isn't a PCA, so the distance between different nodes is not relevant.
For those reasons there are several methods to build these networks, some of them give more emphasis to represent better the degree of the network or other characteristics of the networks.

I consider this representations bad for informing about the networks, because the cut off is arbitrary and the distance between nodes is not informative. One could first do a PCA or MDS and then add the edges between the nodes, this way the relative distance between the nodes would be meaningful as well as the edges. Or one could also use a heatmap to represent all the correlations between the nodes.

The description you quote only states that edge thickness has relevance, not edge length. Usually with this sort of image, the program generating the graph (most probably Cytoscape) will simply arrange the nodes to maximize legibility.
While the edge length can, sometimes, be relevant it usually isn't and instead it is the thickness of the edges and the size of the nodes that carry the information. The positioning of the nodes is most likely arranged to make the image as clear as possible.


text:
What are the differences between GWAS and QTL mapping?
What are the differences between the two methods?
What advantages does one have over the other, and what are their limitations relative to one another?

This question is too broad, and is prone to subjective answers that aren't very helpful. Please try to narrow your question, preferably thinking about a specific task that you would like to carry out.

Most papers that I've looked at distinguish between looking for genotype/trait associations in mapping populations created by experimental crosses (QTL mapping) and looking for genotype/trait associations in populations of unrelated individuals (GWAS). No distinction is made between susceptibility to disease (often referred to as a case-control GWAS) and continuous traits (all of the big studies on genetic associations with height use the term GWAS, for example).
Many of the papers in my first link use both methods, so you can take a look at them to get an idea of the advantages and limitations of each.


text:
Identification of Microorganism
I have the 16S DNA Sequence of a Microorganism. Can I confirm its identity by only running a BLAST or do I have to follow other methods. If so please mention the methods.

Can you please clarify this question? Do you have a single 16s DNA sequence, or multiple sequences? What proportion of the 16s gene is represented by this sequence?

As always with bioinformatics, any reasonable answer will be dependendant on what you want to do with this sequence, and you haven't provided enough information in that regard.
For most exploratory analysis BLAST is good enough. You don't have to follow any existing methods. The benefit of using NCBI's web BLAST server is that it provides you with a very quick way to find a match from the sequence information in GenBank.
BLAST does a local alignment. If you have the full 16s sequence, and it's a novel organism, you might get different match results from a program that does global alignment, but then again you might get different results from a different local alignment tool.


text:
scater: SingleCellExperiment, Error in seq_len(ncol(assay))
I am trying to go through the following scRNA-seq tutorial.
But the line sce <- newSCESet(countData=all.counts) is not working anymore with the most up-to-date version of scater. Now we should use SingleCellExperiment function instead, so when trying to do the following:
sce <- SingleCellExperiment(assays = list(counts = all.counts))

I am getting the error:
Error in seq_len(ncol(assay)) :
  argument must be coercible to non-negative integer
In addition: Warning message:
In seq_len(ncol(assay)) : first element used of 'length.out' argument

Any help would be greatly appreciated.

You can better ask questions about bioconductor packages on https://support.bioconductor.org/, you'll probably get an answer from Aaron Lun himself within a day.

Which version of R and of the packages are you using? That workflow is for version 3.4 of Bioconductor, and the current  version is 3.6

I think you might want to use the most recent workflow, which should be for the latest version on Bioconductor.

You simply need to transform your data.frame into a matrix. Just use these  two lines:
library(SingleCellExperiment)
sce <- SingleCellExperiment(assays = list(counts = as.matrix(all.counts)))