text:
How to loop multiple function in shell script?
I need to extract sequences one after another consecutively from a large fasta files (multiple fasta files) and each extracted files to be saved in new fasta file (I mean the first sequence extracted from multiple fasta file to be saved as gene1 and second sequence extracted from multiple fasta file to be saved in another file and so on). For that, I tried a script like this,
for file in *.fasta do cat "$file" \
  | awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),$0);N++;next;} \
          {printf("%s",$0);} \
     END {printf("\n");}' \
  | awk '{if ( NR==1 ) print $0;}' \
  | awk '{print $1"\n"$2}' | uniq > "new_${file}"

But, its not serving my purpose. Firstly, I could not loop my sequence extraction procedure from multiple fasta file and save as gene1 and gene2 and so on.
For example,

I have fasta files as follows,
org1.fasta
>locus1
ATGCGTAGAG
>LG12
TATGCATGAT
>bpl
TAGATGGAT

org2.fasta
>ysr
ATGTAGCGA
>inc
TGATCTATG
>pym
ATGCATAGA

org3.fasta
>siv
TAGTAGTAT
>tdp
ATGAGTCGA
>tem
ATGCTHATG

I need to extract first positioned sequences from all the org1.fasta org2.fasta and org3.fasta files and save in a new file gene1.fasta and second positioned sequences as gene2.fastaand third positioned sequences to be saved as gene3.fasta.
Expected output,

gene1.fasta
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT

gene2.fasta
>LG12
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA

gene3.fasta
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

Please help me to do the same.
Thanks in advance.

You don't need a loop. You can do the whole thing with a simple awk one-liner:
awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org*fasta

I ran this on your files and got:
$ for file in gene*; do echo "=== File: $file ==="; cat $file;  done
=== File: gene1.fasta ===
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA
>siv
TAGTAGTAT
=== File: gene2.fasta ===
>LG12
TATGCATGAT
>inc
TGATCTATG
>tdp
ATGAGTCGA
=== File: gene3.fasta ===
>bpl
TAGATGGAT
>pym
ATGCATAGA
>tem
ATGCTHATG

Explanation

-vRS='>' : set the record separator to >. This will make awk consider every fasta sequence as a single "line".
FNR>1{foo} : FNR is the line number of the current file. Because > is the first character of the entry, the first "line" will always be empty, so it needs t obe skipped.
printf ">%s",$0 > "gene"FNR-1".fasta" : this does the work. It prints a > followed by the current "line" ($0, this is the entire sequence + the header) and redirects it into a file called gene, then the current line number -1 (-1 is needed because of the issue mentioned above) and a .fasta extension.

=====
Note that org*fasta will be expanded with the default collation order. This means that org10.fasta will be processed before org1.fasta:
$ printf "org1.fasta\norg10.fasta\n" | sort
org10.fasta
org1.fasta

To avoid this and have the files in numerical order, you can use brace expansion instead of the *. So, if you know your files are all orgN.fasta where N can take values from 1 to 20, do this:
awk -vRS='>' 'FNR>1{ printf ">%s",$0 > "gene"FNR-1".fasta"}' org{1..20}.fasta

@K.Dineshkumar you're welcome! If one of the answers here solved your issue, please take a moment to accept it by clicking on the checkmark on the left. That will mark the question as answered and is the way that thanks are conveyed on the Stack Exchange sites.

Thank you so much @terdon. However, I am facing a problem, this script perfectly fine for up to 9 number of fasta files. If it exceed 9, it takes 10th fasta file and printing first instead of 1st fasta file sequence. Same way, it prints 20th fasta file second instead of 2nd fasta file. I do not know how to solve this.

@K.Dineshkumar ah yes, that is because of the order of collation. When you do echo org*fasta, it will sort alphabetically, so org10.fasta comes before org1.fasta (try with printf "org1.fasta\norg10.fasta\n" | sort). Try the updated answer.

Many thanks @terdon for your detailed explanation. With your explanation, I learned how to play around the script, in order to serve our purpose.

The title of your question indicates a shell solution, however you have added R and Python as tags. Here is an R solution making use of Biostrings pacakge (using two input files for simplicity):
library(Biostrings)

# create Biostring objects from each multifasta file
org1 <- readDNAStringSet("org1.fasta")
org2 <- readDNAStringSet("org2.fasta")

# loop over the Biostring objects
for(i in 1:length(org1)){

  output_file <- paste0("gene_", i, ".fasta")

  # create a Biostring object for each index
  genes_at_index_i <- c(org1[i], org2[i])
  assign(gene_name, genes_at_index_i)

  # save Biostring object as a multifasta file
  writeXStringSet(genes_at_index_i, output_file)
}

Output, gene 1:
>locus1
ATGCGTAGAG
>ysr
ATGTAGCGA

Output, gene 2:
>LG12
TATGCATGAT
>inc
TGATCTATG

Warning: The code above won't work if you don't have the same number of fasta entries in each multifasta file.

Thank you @haci for your help. However, I could not use this approach for 100 fasta files or above. Creating variable for each fasta file is tedious job.

You can do so in a loop, i.e. creating a list of Biostring objects with for or lapply and then looping over that list.

Thank you @haci, I would try to do the same.


text:
Deep learning RNA sequences
Currently I'm working on a project, which combines deep learning with RNA sequences. I'll try to predict pseudotorsion angles [1] from raw rna sequence. The ideas is to train a neural network with raw rna sequences and for each nucleotide their corresponding pseudotorsion angles, and than to predict the angles of the remaining sequences in the test set.
How my data is structured:
example:
seq1:   A  C  G    G   U   A   C
Eta:   169 87 110 87  45  187 78
Theta: 123 10 45  168 132 34  100

[1] These are angles describing the backbone conformation of a rna molecule.
I'm pretty new to the field of deep learning, and so far I build a simple feedforward neural network, but it's prediction accuracy is pretty low with only one percent.
Has anyone some tips for me how to improve this?
How do I preprocess this kind of data correctly for deep learning?
I appreciate any help.

Note that RNA has over 110 known modifications, most of which can't be detected using existing automated sequencing tools. These would probably influence the structure of the molecule, and reduce the accuracy of your prediction.

The question of folding of RNA sequences looks slightly similar to protein folding - perhaps searching in this domain might bring more suggestions.
An example of (current) state of the art of deep learning approach to protein folding is AlphaFold developed by DeepMind for CASP competition . They basically decompose the angles to an image representation and look for structures in the proximity of each atom. While it requires a significant processing power, RNA structures might be easier due to the size of the sequence.

I come from a protein background and this problem is analogous to protein torsion angle prediction, which in turn is a variant of protein secondary structure prediction.
Conventional ways to go here would be to use a sliding window of x bases either side to make a prediction for each point, or use a CNN/LSTM-type model to run over the whole sequence and output angles.
There are other technical concerns such as how to split the training and test sets, how to output the angles from the model (usually transformed by sin or cos), and whether to predict a single value or the probability of the value being in certain ranges.
The reference to protein tertiary prediction and CASP in another answer isn't directly relevant here unless you intend to build 3D models from your torsion angles, which you might very well do.
Some references on torsion angle prediction for proteins:

A recent review.
A recent state of the art method.
A classic paper describing some of the approaches.

Thanks for your post and highlighting your work. Could you clarify that CNN is convolution neural network and when do you use it as opposed to LHTM?

Correct, CNN is convolutional neural network and LSTM is long short-term memory, a variant of recurrent neural networks (RNNs). LSTMs would typically be used for sequential data, where you aim to predict a property for each point in the sequence or predict the next point in the sequence.
CNNs apply the same filters to every point in the sequence and rely more on picking up local features than treating the input as a sequence. They are usually applied to images (2D) but can be 1D, 2D, or 3D (etc.) as required.

LSTM, even I know this acroynm. Many thanks that finally makes  sense in context to how CNNs operate in this field and more generally. We're using your (DJ's) stuff, it always nice to know how it works :) . For the OP getting deep learning to work isn't always trivial, this team is expert.

I think the solution is more in the DL.
I would look at simpler ML, using a random forest as a starting point and then maybe SVC. If you get some reasonable predictive power on this then move towards neural network/Tensorflow type calculation.

source: Stack Exchange

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

sort: votes


text:
Any way to quantify the variation of genes that expressed in Affymetrix expression data?
am experimenting preprocessed Affymetrix microarrays expression data matrix (Affymetrix probe-sets in rows (32830 probesets), and RNA samples in columns (735 samples)) for my downstream analysis. Here is how the data looks like:
> dim(eset_HTA20)
[1] 32830   735
## load assayData
HTA20_rma <- load("HTA20_RMA.RData")

> eset_HTA20[1:10,1:3]
             Tarca_001_P1A01 Tarca_003_P1A03 Tarca_004_P1A04
1_at                6.062215        6.125023        5.875502
10_at               3.796484        3.805305        3.450245
100_at              5.849338        6.191562        6.550525
1000_at             3.567779        3.452524        3.316134
10000_at            6.166815        5.678373        6.185059
100009613_at        4.443027        4.773199        4.393488
100009676_at        5.836522        6.143398        5.898364
10001_at            6.330018        5.601745        6.137984
10002_at            4.922339        4.711765        4.628124
10003_at            2.689344        2.771010        2.556756

objective:
I am interested in to reduce the dimension of the dataset by tossing off low expressed in genes in the experiment. To do so, I need to quantify the variation of the gene that expressed. I checked the possible approach and using the coefficient of variation (cv) could be one of the approaches I could try. It is not intuitive for me how to quantify the variation of genes that expressed.
Can anyone point me out how to lay out this procedure on the above expression data? Is there any Bioconductor package to do this kind of task? any idea to make this happen? thanks

You can use the following code to calculate the coefficient of variation:
# expr is your expression matrix.
SD <- apply(expr, 1, sd)
CV <- sqrt(exp(SD^2) - 1)

It might be implemented in some package but it is so brief that you can write it again yourself. Then you can filter out those that are below certain percentage of the distribution of CV (like 0.1)
expr[CV > quantile(CV, probs = 0.1), ]

Another way to reduce dimensionality of your dataset is via principal component analysis or PCA. It is implemented as pcaMethods R package, for instance.

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


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

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

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

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

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

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

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

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

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

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

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

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


text:
How to use extractfeat?
I am working on a script on bash linux to get the CDS (coding sequence) of a gene using extractfeat by EMBOSS. It gives me the error: Warning: No sequences written to output file
I am unsure if my syntax is right, but this is what I have:
extractfeat inputFile outputFile -type CDS

Hi, KA97, what is the error? If you can include it, this will help establish what the problem is.

. The error is :   Warning: No sequences written to output file
It produces a blank output file with no sequence in.

Perhaps,removing introns (if it is a higher eukaryote)
extractfeat tembl:x65921 -type CDS -join

It is difficult to know without knowing the gene.
BTW NCBI's Entrez "efetch" will do essentially the same thing, but its more complicated

I tried that syntax with the -join at the end which didn't work either
If it matters, I am currently testing it on a Histone H1 gene of C. elegans
I'd prefer If I could get the CDS by processing the sequence file in my bash script. Since I'm using efetch to download the sequence anyway, I don't mind trying that out if it's not too complicated.

Okay, the syntax you are using is correct, but you don't need to specify the output file

I managed to fix it. I was accidentally using FASTA format files instead of Genbank flat file format which seemed to be the problem.

That's good: fasta format has no features for extractfeat to extract, so it would have no chance. The error trapping in the code is clearly absent and it might be worth reporting the bug to the Emboss team.


text:
What are the mathematical assumptions and algorithmic procedure behind MACS2?
I wonder if there's some mathematical details that outline what the algorithms does.

Hi and welcome, it would be helpful if you could put your question in context. What is the objective, how does the algorithm help and why it's published information is insufficient for you.

Have you read the papers?

See here and here.

@TheSimpliFire Thanks.


text:
Biological replicates on Chip-seq Transcription factor data
I have 2 biological replicates for chip-seq transcription factor data and I have applied macs2 peak calling for each replicate separately.
How can i make use of the biological replicates to extract the most out the data?
Thanks!

Hello Dimtiris, that's quite a broad question, is there a way to narrow it? What authors of macs2 say? How it's treated in other publications? How do you think you should treat your two replicates?

This question is somewhat generic, so a generic answer is that ENCODE has a Transcription Factor ChIP-seq Data Standards and Processing page that can give you a useful starting point.
For TF ChIP-seq data with replicates, the Irreproducible Discovery Rate (IDR) method helps leverage replicates to produce higher confidence peak calls, producing both "optimal" and "conservative" peaks, defined by some IDR threshold. Here's ENCODE's description of the method:

A statistical procedure that operates on the replicated peak set and compares consistency of ranks of these peaks in individual replicate/pseudoreplicate peak sets. Peaks with high rank consistency are retained. IDR can operate on peaks across a pair of true replicates resulting in a “conservative” output peak set, or across a pair of pseudoreplicates resulting in an “optimal“ output peak set. Peaks in the conservative peak set can be interpreted as high confidence peaks, representing reproducible events across true biological replicates and accounting for true biological and technical noise. Peaks in the optimal set can be interpreted as high-confidence peaks, representing reproducible events and accounting for read sampling noise. The optimal set is more sensitive, especially when one of the replicates has lower data quality than the other.

The Python code to run IDR on peak sets can be found on ENCODE's GitHub Chip-seq pipeline repo.

You could first look at the degree of correlation between the two replicates - what proportion of peaks are shared between the two, versus peaks found in one sample only?  This will give an idea of how repeatable the analysis is, and how many peaks are variable between samples or due to artefacts of the method.
When it comes to interpreting the biology behind the binding patterns, you are probably going to want to only interpret peaks which you find in both replicates.


text:
How to convert BAM file to bigWig in Python?
I am manipulating some BAM files using pysam package which seems to be very fast and handy. However, I ran into problem when I am trying to generate bigWig files from BAM files in Python. I used bamCoverage from deeptools to do this previously, but it is a standalone tool that cannot be used from within my Python scripts.
Therefore, I am wondering what Python package would allow me to read in BAM files and convert them to bigWig files.
Thanks

deepTools has a (somewhat poorly documented) API, since it's a python package too. The basic code framework is:
from deeptools import writeBedGraph
from deeptools.getScaleFactor import get_scale_factor

wr = writeBedGraph.WriteBedGraph(options and input files)
func_args = {'scaleFactor': get_scale_factor(...options...)}
wr.run(writeBedGraph.scaleCoverage, func_args, output_file_name, ...)

Thanks, I am wondering what this scale factor means. Is it possible to calculate a RPKM or CPM like bamCoverage?

The scale factor dictates how the coverage track is transformed before writing the results to a bigWig file. Among other things, this then handles RPKM and other normalizations. Please notes that you'll have to pass an argparse object into it, which I understand is a bit annoying (I had't considered external use of that function).

pyranges:
import pyranges as pr
gr = pr.read_bam("your_file.bam")
gr.to_bigwig("out.bw", chromosome_sizes=pr.data.chromsizes()) # for hg19


text:
How can extract the list of genes name from the raw data in GEO?
The format of the raw data that upload in GEO is different, like TXT, CEL or other form. How can I extract the list of genes name from these formats? And then how can I normalize them with R?
I would be really glad if anyone write step by step a complete set of codes for R with an example.

Answer from @benn, converted from comments:
You can find a step by step description in the limma user's guide. Just download the CEL files, and follow the user guide.
[please edit to improve this answer, if possible]


text:
Accuracy of large peptide fragment in LC-MS/MS
Peptide digestion leads to generation of small peptides (ideal) as well as some large peptides (3-4kDa). The large peptides often don't fly well and therefore do not detect well. Does it affect the mass accuracy (delta ppm error) of the detected peptide? I am taking following rationale here: lower peptides are well detected and therefore a high confidence data can be obtained. While in larger peptides that confidence cannot be obtained due to low detection).


text:
Fragment mass tolerance in proteomics analysis
Proteomics newbie here. I am using proteome discoverer software to analyze the fragmentation of a chemically modified single protein. In my understanding, when we set the Fragment mass tolerance to larger value, we should observe more peptide coverage/modification compared with stricter mass tolerance (for example 1Da vs 0.02 Da).I played with the fragment mass tolerance a bit keeping all other parameters constant. To my surprise, I observed better peptide coverage and newer modifications with lower mass tolerance compared to larger value. How this is possible?

I agree with @DanielGeiszler's explanation. Maybe some illustration may help.
Software typically use a scoring function to separate the true and false identifications. In the illustration below, all peptides to the right of the cut off line are the one to keep (meaning within 1% FDR)

Now using a large fragment ion threshold will artificially increase the scores of the false peptides (it may affect the true peptides as well, but to a much lesser degree). The distribution will now look like below, where the false peptides' curve shifted to the right. To maintain the same FDR cut off, the line needs to be drawn to the right as well, thus less confident peptide identification.

If your data is high mass accuracy, there isn’t a benefit to using low mass accuracy in your search. You should be hitting nearly all of the true peptide fragment ions with your tight tolerances.
When you use larger mass tolerances, you’re going to match more fragment ions for peptides that are incorrect matches. The result is that incorrect peptide hits score higher while correct peptides don’t increase their scores, so it becomes harder to separate your false peptide hits from your true peptide hits. Then, when you go to filter to a given FDR, this manifests as fewer true peptides passing FDR filtering.
MSFragger (within FragPipe) implements a procedure to optimize your fragment tolerances described briefly here: https://www.nature.com/articles/s41467-020-17921-y .
Full disclosure: I’m an author on the above paper.

Thanks for the answer. "Then, when you go to filter to a given FDR, this manifests as fewer true peptides passing FDR filtering" I didn't fully understand this text. Could you elaborate a bit?


text:
Pre-filtering genes for Principal Component Analysis
I have a raw counts data-set of 20,502 genes and 137 samples. I want to find out Principal Components which best explain variation between samples in different stages of tumor.
I am new to Machine learning and would like some help in selecting a criteria to reduce the number of genes before doing a PCA.
I have tried doing a DESeq2 analysis on the data-set but I am confused as to what to do next.

You can filter our genes that do not show variation across samples, they would not be helpful in differentiating samples anyhow.

So you mean genes which do not fit the dispersion curve in DESeq2?

More like they did here: http://pklab.med.harvard.edu/scw2014/subpop_tutorial.html

This is really helpful. Thanks a lot :)

As you found the approach useful, I will re-write this as a formal answer for future reference.

You can filter out genes that do not show variation across samples, they would not differentiate samples anyway. For the specifics on how to do so, please see the Subpopulation Analysis section of a nice single cell RNA-seq workflow from the Kharchenko Lab. I believe there won't be fundamental differences when applying this approach to bulk RNA-seq.

I think "the best" way to do this is going to depend on what your question is, unfortunately.
For example, let's say you're trying to find which genes vary by tumour stage.
Variance of a given gene might be high across all samples, in which case it has no relation to the stage and might be useless for you.
Alternatively, you may have a lowly-to-moderately expressed gene that does not vary much between samples of the same stage, but differs between stages, which is useful for you.
However, because it's not highly-expressed, its variance may be small compared to the variance of highly-expressed genes and skipped over.
But if you're trying to do sequencing batch correction to identify batch effects, looking at raw sequencing counts and doing some type of quantile normalization might be all you need.
So I don't think there's a one-size-fits-all answer to your question, unless you're able to refine it a little bit


text:
Downloading dataset from SRA (SOLiD Platform)
Trying to download colorspace data from SRA, but getting an error
abi-dump -A SRR1657115.sra
abi-dump.2.9.6 err: item not found while constructing within virtual database module - the path 'SRR1657115.sra' cannot be opened as database or table.


text:
Getting Unique Identifier List for GEO Datasets NCBI
AIM: Download "Unique Identifier List" for the following query from GEO DataSets.
Query: ("Expression profiling by high throughput sequencing"[DataSet Type] AND ("Homo sapiens"[Organism] OR "Mus musculus"[Organism] OR "rattus norvegicus"[Organism])) AND ("2020/01/01"[PDAT] : "3000"[PDAT])
which means, all RNASeq studies deposited on GEO in the year 2020 for humans, mice or rat.
Problem: I need the GSE ID list for ~9k datasets, but while trying to download the list of ids, it loads to a blank page and nothing happens. Also, clicking on "Next Page" gives error.
I have been trying for the last 3-4 days but it doesn't work.
Steps to generate file:
"Send To" -> "File" -> "Format" (Unique Identifier List) -> "Sort By" (Default Order) -> "Create File"

You can use Entrez Direct for this. The following returns Unique Identifiers which are just bare integers.
$ geo_query='"Expression profiling by high throughput sequencing"[DataSet Type] AND ("Homo sapiens"[Organism] OR "Mus musculus"[Organism] OR "rattus norvegicus"[Organism]) AND ("2020/01/01"[PDAT] : "3000"[PDAT])'
$ esearch -db gds -query "$geo_query" | efetch -format uid > gds_results.txt
$ wc -l gds_results.txt
9981 gds_results.txt
$ head -n2 gds_results.txt
200134092
200120931

Instead, if you are looking for a way to get the GSE accessions, you can use the built-in xtract command to parse the XML returned by esummary as follows:
$ esearch -db gds -query "$geo_query" | esummary | xtract -pattern DocumentSummary -first Accession > gse_accs.txt > gse_accs.txt
$ wc -l gse_accs.txt
9981 gse_accs.txt
$ head -n2 gse_accs.txt
GSE165829
GSE165824

I think I would just try to do this with GEOquery.
https://bioconductor.org/packages/release/bioc/html/GEOquery.html

GEOquery would download datasets for a given accession ID. R provides a library https://www.bioconductor.org/packages/release/bioc/html/GEOmetadb.html which contains metadata for microarray studies. I couldn't find anything similar to this for RNASeq studies so had to resort to the GUI.

Oh you're right, I did something similar before, i though it was GEOquery but i guess i don't remember... Ok, can you download the text file like:
"Send To" -> "File" -> "Format" Summary (text)
And then do something like:
cat gds_result.txt | grep 'Accession: GSE' | cut -f 3 | cut -f 2 -d ' ' | sort > gse.txt


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

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

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

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

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

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

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

1/2^bitscore

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

(query length * database length * p-value)

which simplifies a bit to:

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

I hope this helps!

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

Use a Poison or binomial distribution.
Use randomization

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


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

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

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

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

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

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

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

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

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


text:
Creating msk file for Illumina Array Analysis CLI Copy number train
recently I'm researching on how to detect CNV from Illumina microarray data. I found that Illumina has Array Analysis CLI software for this task, and I need to train the CN model myself, since the official .dat file is no where to be found.
However, this requires another file: the manifest file (.bpm) must be accompanied by a "mask file" (.msk), which I cannot find online.
Could you please help me on this issue please? And could you suggest some CNV detector for Illumina microarray? Thank you.


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

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

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

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

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

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

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

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


text:
bedtools: get name column of alignment file A
I'm using the coverage function of bedtools to check whether a certain set of regions (file B) has any overlap with known regions of epigenetic markers (file A). When that is the case, I want to access the name of that marker, which is given as the third column (name) in file A.
However, when I run bedtools coverage -a fileA.bed -b fileB.bed, I only get the columns of file B as an output, together with the overlap fraction. Is there a built-in way of accessing also the columns of the file A to which I aligned, or an easy way of cross-referencing to it?

I would recommend using bedtools intersect for that type of comparison.  This will allow you to see the specific segments that overlap.  bedtools coverage is used more for determining how many things cover a region but not which ones.
The documentation linked above will provide a much better desciption of intersect than I would be able to in this post but I would recommend using the -loj or -u depending on what you want.
-loj    Perform a “left outer join”. That is, for each feature in A report each overlap with B. If no overlaps are found, report a NULL feature for B.

-u  Write original A entry once if any overlaps found in B. In other words, just report the fact at least one overlap was found in B. Restricted by -f and -r.

If you're working with fileB.bed (a sorted BED file) and fileA.txt (a three-column text file, where the name is in the third column), then you could pipe single-base pseudo-BED elements to bedmap.
$ bedmap --echo --echo-map-id --delim '\t' fileB.bed <(awk -vFS="\t" -vOFS="\t" '{ print $1, $2, ($2+1), $3 }' fileA.txt) > answer.bed

The --echo-map-id option reports the list of marker IDs from fileA.txt, where the single-base marker interval from fileA.txt overlaps a region in fileB.bed.
If all the marker IDs are not unique (if there are duplicate marker IDs) it may make sense to instead use --echo-map-id-uniq.
This assumes that marker positions in the second column and that the intervals in fileB.bed are using a 0-based index scheme. If marker positions are 1-based and fileB.bed is 0-based, then a small adjustment is needed to make both inputs use the same indexing, to calculate correct overlaps:
$ bedmap --echo --echo-map-id --delim '\t' fileB.bed <(awk -vFS="\t" -vOFS="\t" '{ print $1, ($2-1), $2, $3 }' fileA.txt) > answer.bed

I mention this, as SNVs (such as those from VCF) can often be reported with a 1-based index scheme, which adds some nuance to doing set operations against BED files.


text:
Convert or extract Seurat object as the input of FateID
I am currently using Seurat to do the scRNA-seq clustering analysis. After this, I am planning to use FateID to do the downstream analysis, could someone teach me how to convert or extract a Seurat object as the input for FateID?
Thanks in advance.

Seurat objects have various slots where different types and forms of data are stored. For FateID you would need an expression matrix (with cells as columns and genes as rows) and this can be extracted via:
seurat_object@assays$RNA@data

Here RNA is the assay (you can have multiple assays per sample, RNA, protein, ...) and data is normalized data. You can extract counts or scaled data instead.
If you are using Seurat v2, drop the @assays$RNA part.

Dear haci, thx, I will try.

Thank you @haci, I have figured out it by myself.


text:
Where to Download Cancer Raw Reads (fastq)?
Does anyone know where I can download cancer raw reads (fastq files), tumor and Germline for non-humans?
I wanted to make a study with human data but I don't have the control access to download raw reads.

So do you want human or non-human? Your question is confusing. For human, have a look at Is there a publicly available tumor-normal sample?. Is that what you need?

Please tell us when you cross-post your question. This question is also asked on biostars.

You can try to browse SRA (sequence read archive) at https://www.ncbi.nlm.nih.gov/sra.
To download fastq files from SRA, you'll need their sratoolkit: https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software.
The easy command to download fastq files by their accession number is :
fastq-dump --split-files accession_number


text:
How do you go about grabbing proper genotyping infomation from dbSNP using entrez API?
I made an R function to grab rsID genotypes from dbSNP using the entrez API. For example, here is the xml output that I am trying to parse:
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=snp&id=rs2656176&retmode=xml
the 'SPDI' tag contains the genotype info. so I can easily use xpathSApply to get this tag for many IDs. But for this dbSNP there are two different values, T>C and T>A.
The T>C change is the one that has more support in all pop freq databases, and this is easily grabbed from the dbSNP webpage:
https://www.ncbi.nlm.nih.gov/snp/rs2656176
The "Alt Allele" frequencies show A being nearly 0, while C is the major alt allele.
However, I don't understand how I can parse the data from the XML from above using this... the 'global maf' tag in the XML shows some MAF values, but the 1000G value is only for the reference allele, not the alt.. There are other dbs that show the alt "C" allele but that might not always be the case.
Does anyone have an idea of what tags I should look at? Or is there a better way to hit a different API to get this data instead?
Thanks!


text:
Is this the correct method for determining Power at a given alpha, VAF, and coverage threshold?
I'm trying to calculate power at a given VAF, target coverage, and alpha. For example, what would my power be for a 100x coverage site, detecting VAF of 0.10 with an alpha of 0.05. This is what I came up with, in R, but I would like confirmation that this is the correct way to proceed.
coverage <- 100
VAF <- 0.10
alpha <- 0.05
alternative_VAF <- 0.5

mean <- coverage * VAF
std_dev <- sqrt(coverage * VAF * (1 - VAF))
critical_value <- qnorm(1 - alpha, mean, std_dev)

alt_mean <- coverage * alternative_VAF
alt_std_dev <- sqrt(coverage * alternative_VAF * (1 - alternative_VAF))
power <- 1 - pnorm(critical_value, alt_mean, alt_std_dev)

In theory this should allow me to calculate the power of an assay.

Please [edit] your question and start from the beginning. What are you doing? Are you analyzing genetic variants? What is "power" in this context? Why do you want to calculate it? What will it tell you? We can't know if it's the right way to proceed unless we know where it is you want to proceed towards.


text:
bedtools unable to open file despite being tab-delimited
I got the following error when attempting to intersect two files. The file below is supposed to be the -a file.
Error: Unable to open file compare_me_to.txt. Exiting.

When I attempted to make it tab-delimited through Ubuntu (see code below)
perl -p -i -e 's/ /\t/g' compare_me_to.txt

I got the following error:
Can't open compare_me_to: Permission denied.

Additionally, the ownership of the folder the files are in is of my username.
I can't tell if this is a bedtools problem or a Ubuntu problem or something else.

Just a quick check since you said you have ownership of the file... The file isn't open anywhere else (like in another shell or on desktop) when you try to convert right?

I suspect the reason is the "-i". The following will work,
perl -pe 's/  */\t/g' compare_me_to.txt > compare_me_to.txt2

-i will want to write onto the file and create a backup. I suspect you'll have spotted chmod +w compare_me_to.txt.
Just make sure that a double space is inserted before the * character, otherwise it will go wrong. Copying and pasting will avoid this.

This is a clearer solution because it avoids missing a space character,
perl -e 's/[ ][ ]*/\t/g' compare_me_to.txt > compare_me_to.txt2

and the solution you used .. (I did on OS X), what I'd recommend is,
perl -p -e 's/[ ][ ]*/\t/g' compare_me_to.txt # to see if it works by printing on screen without altering the file

perl -i -e 's/[ ][ ]*/\t/g' compare_me_to.txt

The  -i flag will alter the original file and store a backup without requiring renaming the file.
What's it doing?
The regex [ ][ ]* structure is saying look a space and if one or more spaces occur adjacent to each other replace everything with a single tab. [ ]* means a "space" that is present or absent and if present it can be present an infinite number of times. If the [ ]* is by itself the command risks saying replace '' (zip) or 'space' with a tab, so tabs can appear at every '' (zilch), so tabs can appear everywhere. The "*" symbol in regex work differently to the "*" wildcard character in shell.

The following error appears: Can't open compare_me_to.txt: Permission denied.

I was able to give myself permission (https://askubuntu.com/questions/466605/cannot-open-output-file-permission-denied), but the file still can't be opened when intersecting. I ran the line of code you wrote; when I opened the resulting file, every character was separated by a tab.

Oh I know, you put a single space into the term ... its a double space followed by a *. so its ... 's/  */\t/g' ... if you only put a single space into the first / */ then it will result in tabs going everywhere. I can explain later.

chmod 700 file_name is another way in, I'm surprised chmod +w file_name didn't work if you were the file owner. I assumed you knew about chown and chmod


text:
Truble to run a multiprocessing kmer count script
Hi there I have this code:
from functools import partial
import multiprocessing

def kmer_count(sequence, alphabet, k):
    """Returns a dictionary with kmers and it counts."""
    seq = sequence.upper()
    seq_len = len(seq)
    kmers = [seq[i:i+k] for i in range(0, seq_len - k + 1)]
    filterd_kmers = [kmer for kmer in kmers if
                     all(base in set(alphabet) for base in kmer)]
    return Counter(filterd_kmers)

pool = multiprocessing.Pool(4)

args = [''.join(alphabet.iupac_dna), 4] # in the real code it came from argparse
f = partial(kmer_count,seq)
pool.starmap(f, *args)

Result:
---------------------------------------------------------------------------
RemoteTraceback                           Traceback (most recent call last)
RemoteTraceback:
"""
Traceback (most recent call last):
  File "/opt/anaconda3/lib/python3.7/multiprocessing/pool.py", line 121, in worker
    result = (True, func(*args, **kwds))
  File "/opt/anaconda3/lib/python3.7/multiprocessing/pool.py", line 47, in starmapstar
    return list(itertools.starmap(args[0], args[1]))
TypeError: kmer_count() missing 1 required positional argument: 'k'
"""

The above exception was the direct cause of the following exception:

TypeError                                 Traceback (most recent call last)
<ipython-input-79-0f5ba8381a51> in <module>
      1 f = partial(kmer_count,seq)
----> 2 pool.starmap(f, *args)

/opt/anaconda3/lib/python3.7/multiprocessing/pool.py in starmap(self, func, iterable, chunksize)
    274         `func` and (a, b) becomes func(a, b).
    275         '''
--> 276         return self._map_async(func, iterable, starmapstar, chunksize).get()
    277
    278     def starmap_async(self, func, iterable, chunksize=None, callback=None,

/opt/anaconda3/lib/python3.7/multiprocessing/pool.py in get(self, timeout)
    655             return self._value
    656         else:
--> 657             raise self._value
    658
    659     def _set(self, i, obj):

TypeError: kmer_count() missing 1 required positional argument: 'k'

If I run like this:
f(*args)

I got this:
Counter({'GTTC': 9780,
         'TTCG': 11291,
         'TCGC': 14803,
         'CGCC': 14428,
         'GCCA': 19254,
         'CCAG': 11512,
         'CAGA': 11641,
         'AGAG': 10387,
         'GAGC': 13164,
         'AGCG': 16025,
         'GCGG': 13067,
         'CGGT': 12329,
         'GGTT': 14217,
         'GTTT': 18733,
         'TTTT': 25488,
         'TTTG': 22006,
         'TTGA': 20004,
         'TGAC': 11351,
         'GACT': 7262,
         'ACTA': 5663,
         'CTAG': 4190,
         'TAGC': 8543,
         'AGCT': 12898,
         'GCTT': 19114,
         'CTTG': 16557,
         'TGAA': 17703,
         'GAAC': 9334,
         'AACA': 14050,
         'ACAC': 8812,
         'CACA': 11050,
         'ACAT': 9705,
         'CATC': 17357,
         'ATCC': 11157,
         'TCCC': 5754,
         'CCCG': 5926,
         'CCGT': 8601,
         'CGTC': 7401,
         'GTCC': 3811,
         'TCCT': 5890,
         'CCTG': 7985,
         'CTGC': 15324,
         'TGCG': 14856,
         'TTTA': 15607,
         'TTAG': 8366,
         'TAGG': 5393,
         'AGGC': 9580,
         'GGCA': 13939,
         'GCAC': 11686,
         'CACC': 14496,
         'ACCA': 16950,
         'CCAC': 13028,
         'CACT': 13194,
         'ACTT': 12585,
         'TTGC': 18877,
         'TGCA': 13770,
         'GCAT': 13307,
         'CATT': 14102,
         'ATTA': 9776,
         'TTAA': 12585,
         'TAAG': 9048,
         'AAGC': 18834,
         'AGCC': 12734,
         'CCAT': 13799,
         'CATA': 8188,
         'ATAA': 12278,
         'GCGT': 12546,
         'CGTT': 12599,
         'GTTG': 15136,
         'TTGG': 18903,
         'TGGC': 19391,
         'GGCC': 7685,
         'GCCC': 7404,
         'CCCA': 9752,
         'CCAA': 18755,
         'CAAC': 14852,
         'AACT': 11886,
         'ACTC': 10028,
         'CTCA': 13518,
         'TCAC': 15385,
         'ACCG': 12027,
         'CCGC': 13105,
         'CGCA': 14853,
         'ATAG': 7066,
         'AGGG': 5273,
         'GGGG': 4927,
         'GGGT': 7972,
         'GAAA': 16511,
         'AAAC': 18215,
         'ACTG': 12639,
         'TGCT': 16663,
         'CTTT': 18228,
         'TTTC': 17497,
         'TCGG': 10116,
         'TGAT': 20593,
         'GATG': 17301,
         'ATGA': 14218,
         'GATC': 14479,
         'ATCG': 16604,
         'CGCG': 13405,
         'CGTA': 8086,
         'GTAG': 6829,
         'TAGA': 6656,
         'AGAT': 11833,
         'ATCA': 20344,
         'TCAG': 13868,
         'ATCT': 11892,
         'TCTT': 14706,
         'CTTC': 13897,
         'TTCC': 9522,
         'TCCG': 7099,
         'CCGG': 5935,
         'CGGG': 5786,
         'GGGA': 5615,
         'GGAT': 11054,
         'GATA': 12085,
         'ATAT': 9205,
         'TATC': 11941,
         'TCGT': 8802,
         'GTAT': 8661,
         'TATT': 12620,
         'ATTT': 18149,
         'TGAG': 13150,
         'GAGA': 8658,
         'GAGG': 5644,
         'GCCT': 9673,
         'CCTT': 11377,
         'CTTA': 9100,
         'TTAC': 11418,
         'TACG': 8250,
         'ACGC': 12774,
         'CACG': 10171,
         'ACGA': 8951,
         'CGAT': 16617,
         'ATTG': 17372,
         'GGCG': 14257,
         'TTGT': 13287,
         'TGTA': 7808,
         'ACGG': 8529,
         'CGGA': 7109,
         'ATGT': 9975,
         'TGTC': 8286,
         'GTCG': 9088,
         'TCGA': 11986,
         'GATT': 16232,
         'GACC': 7608,
         'CATG': 11552,
         'ATGG': 13703,
         'TGGG': 9728,
         'CGAC': 9020,
         'GACG': 7394,
         'ACGT': 8085,
         'GTAA': 11360,
         'TAAT': 9895,
         'AATA': 12342,
         'GGCT': 12713,
         'AAAT': 17834,
         'AATC': 16220,
         'TTCA': 18526,
         'TCAA': 20304,
         'CAAA': 21752,
         'AATG': 14200,
         'TGGT': 16952,
         'GTTA': 10179,
         'GCGA': 14624,
         'AGGT': 9338,
         'GGTA': 10326,
         'CGGC': 12423,
         'TAAC': 9956,
         'AACG': 12524,
         'GCGC': 18306,
         'CTGG': 11344,
         'TGGA': 10805,
         'GCTG': 17455,
         'CTGA': 13532,
         'CGAG': 8649,
         'GAGT': 10014,
         'AGTA': 8668,
         'GTAC': 8038,
         'TACC': 10217,
         'ACCT': 9204,
         'ATAC': 8527,
         'TACA': 7453,
         'ACAA': 12591,
         'AAAA': 24650,
         'AAAG': 17848,
         'AAGG': 11240,
         'AGTT': 12431,
         'CAAT': 17456,
         'GTCT': 6235,
         'AGGA': 5937,
         'GGAG': 5294,
         'AGAA': 12516,
         'AAGA': 13980,
         'CGTG': 9948,
         'GTGC': 11762,
         'TCAT': 14388,
         'AGTC': 7451,
         'GTCA': 11557,
         'CAGT': 13149,
         'GGTC': 7553,
         'GTGA': 15043,
         'TAAA': 15634,
         'GGGC': 7200,
         'GGTG': 14495,
         'TTCT': 13230,
         'TCTA': 6621,
         'CTAC': 6510,
         'TACT': 8573,
         'CTAT': 6856,
         'TATG': 7993,
         'ATTC': 11395,
         'GAAG': 13449,
         'AAGT': 12578,
         'GCTA': 8346,
         'CTAA': 8310,
         'TAGT': 5859,
         'GACA': 8060,
         'TTAT': 12293,
         'TCCA': 11138,
         'CAGC': 17389,
         'GCTC': 13661,
         'CTCG': 8724,
         'CCCC': 5036,
         'TGCC': 13824,
         'GCCG': 12340,
         'GGAA': 9364,
         'CTGT': 9651,
         'TGTT': 14580,
         'ATGC': 13150,
         'AGTG': 13050,
         'GTGT': 9052,
         'TCTG': 11772,
         'AATT': 13738,
         'CGAA': 11026,
         'GTGG': 12927,
         'TCTC': 8866,
         'ACAG': 9505,
         'CCTC': 5686,
         'CTCC': 5391,
         'CCTA': 5236,
         'CGCT': 16302,
         'AACC': 13897,
         'ACCC': 8037,
         'GCAA': 18697,
         'CAAG': 16287,
         'AGCA': 16273,
         'GCAG': 15145,
         'CCGA': 9751,
         'CAGG': 7851,
         'TATA': 4461,
         'GGAC': 3754,
         'CTCT': 10608,
         'TGTG': 11291,
         'AGAC': 6199,
         'CCCT': 5517,
         'GAAT': 11315}

What I am doing wrong with the multiprocessing step. I want to speed up my analysis, because I have a lot of genomes in a directory to process.
I just need to figure out how can I implement the pool, because probably all other functions a similar.
Thanks for your time and attention

I would guess that the kmer counting will be much faster with tools like KMC. Is there any particular reason you want your own counter in python?

I know what you mean, but I am learning to write code and I just want to use python to do that! There are many tools out there like khmer, jellyfish, kat etc... Thank you

I guess that's fair. The question is related to this one about kmer counting in python. Although it does not specificly talk about multiprocessing, I thought it's still relevant.

counting the mers it is not the problem. My problem is  just the multiprocesing step thats giving me a hard time

This code help me to get the answer I was looking for.
r = pool.apply_async(f, a)
print(r.get())


text:
Using python to call Emboss wordcount
Hi guys I would like to run wordcount from Emboss in a bunch of fasta files. I tried this:
#!/usr/bin/env python
import sys
import subprocess

filenames = ['Downloads/Data/GCF_000016725.1_chr.fna']
for filename in filenames:
    name = filename.split('/')[2]
    print(filename)
    subprocess.call (['wordcount', 'Input sequence(s): ' + filename, 'Word size [4]: ' + sys.argv[1]])

But I got error messages:
Count and extract unique words in molecular sequence(s)
Error: Failed to open filename 'Input sequence(s): Downloads/Data/GCF_000016725.1_chr.fna'
Error: Unable to read sequence 'Input sequence(s): Downloads/Data/GCF_000016725.1_chr.fna'
Died: wordcount terminated: Bad value for '-sequence' and no prompt

But if a use shell with Downloads/Data/GCF_000016725.1_chr.fna  it works fine!
Any tip would be helpful!
I am using python, but I would have a method to compare my results!
Thank you by your time!
Paulo

You have to use command line options, not interactive prompts, like so:
subprocess.call(['wordcount', '-sequence', filename, '-wordsize', sys.argv[1]], '-outfile', filename + '.wordcount.txt')

Or, using the f-strings or formatted string literals:
subprocess.call(f'wordcount -sequence {filename} -wordsize {sys.argv[1]} -outfile {filename}.wordcount.txt')

Note that I had to add -outfile option to your system call, otherwise wordcount prompts for the output file. See wordcount manual.

I tested your first solution and worked very nicelly. Thank you by your time and patience. Paulo


text:
How do protein-protein interaction databases work?
I am a student and just trying to gain a base level understanding of PPI databases and how to interpret their data.
For example, when inputing multiple genes into String, does it only give protein-protein interactions between those identified in my gene list? Or does it compare the input gene list to a background of your selected species? If it's the former, is it normal to only want protein-protein interactions between those in your input list?
Is there a resource/text where I can find more information about PPI databases? Any help in the right direction or where I can learn more would be appreciated.


text:
ipfam or pfam domain-domain interactions
I'd like to download domain-domain interactions that were in ipfam (which doesn't seem to be up-and-running anymore).
for an example of what I would like to be able to download via ftp or another source, see this list of interactions with Death domain:
https://pfam.xfam.org/family/PF00531#tabview=tab8
Is this available in one of these files in the current release??
ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release

I think you're looking for ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/database_files/pfamA_interactions.txt.gz. Based on the file's name, that's a list of domain interactions. The file has two columns, so presumably each line represents an interaction.
Here's what I get for the Death domain (PF00531):
$ zgrep PF00531 pfamA_interactions.txt.gz
PF00036 PF00531
PF00531 PF14786
PF00531 PF00531
PF00791 PF00531
PF01335 PF00531
PF13499 PF00531
PF13833 PF00531
PF14786 PF00531

Those are 8 interactions, which matches what the PFam link you gave in your question shows. If you want to convert the accessions to something more informative, you could try Pfam's REST API and do something like:
$ zgrep PF00531 pfamA_interactions.txt.gz |
    while read a b; do
    name1=$(wget -O - "https://pfam.xfam.org/family/${a}?output=xml" 2>/dev/null |
     perl -00lne 'm|<description>.*?CDATA\[\s*(.*?)\s\].*?</description>|s && print $1');
    name2=$(wget -O - "https://pfam.xfam.org/family/${b}?output=xml" 2>/dev/null |
     perl -00lne 'm|<description>.*?CDATA\[\s*(.*?)\s\].*?</description>|s && print $1');
    printf "%s\t%s\t%-20s\t%-25s\n" "$a" "$b" "$name1" "$name2";
done
PF00036 PF00531 EF hand                 Death domain
PF00531 PF14786 Death domain            Tube Death domain
PF00531 PF00531 Death domain            Death domain
PF00791 PF00531 ZU5 domain              Death domain
PF01335 PF00531 Death effector domain   Death domain
PF13499 PF00531 EF-hand domain pair     Death domain
PF13833 PF00531 EF-hand domain pair     Death domain
PF14786 PF00531 Tube Death domain       Death domain

Obviously, if you need anything more fancy than that, it's time to write a proper parser using Perl or Python or any other language that can understand XML.

This is it! I had found the pfamA_interactions.txt.gz file but didn't know how to get the IDs from accession. This is great thanks.

@Beeba you're welcome :) But this is going to be horribly slow and inefficient if you have more than a dozen or so IDs. If you need to deal with many, do consider writing a proper parser. If you need help with that, you can always ask another question here!

Thanks will do, I'm pretty proficient with using R's xml package so hopefully I've got it.


text:
I get the following error when using the GSEA software 'Enrichment Map Visualization' tool using the Cytoscape 3.7.1 application
I am trying to visualize a GSEA analysis I conducted using the Broad Institute's GSEA Enrichment Map Visualization, however I get the following error and am not sure what the issue is: any help would be greatly appreciated.
INFO - successfully connected to cytoscape rest
INFO - successfully connected to cytoscape rest
INFO - Found enrichment map command
INFO - Get URL:http://localhost:1234/v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt
DEBUG - CookieSpec selected: default
DEBUG - Auth cache not set in the context
DEBUG - Connection request: [route: {}->http://localhost:1234][total kept alive: 0; route allocated: 0 of 2; total allocated: 0 of 20]
DEBUG - Connection leased: [id: 4][route: {}->http://localhost:1234][total kept alive: 0; route allocated: 1 of 2; total allocated: 1 of 20]
DEBUG - Opening connection {}->http://localhost:1234
DEBUG - Connecting to localhost/127.0.0.1:1234
DEBUG - Connection established 127.0.0.1:51968<->127.0.0.1:1234
DEBUG - Executing request GET /v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt HTTP/1.1
DEBUG - Target auth state: UNCHALLENGED
DEBUG - Proxy auth state: UNCHALLENGED
DEBUG - http-outgoing-4 >> GET /v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt HTTP/1.1
DEBUG - http-outgoing-4 >> Host: localhost:1234
DEBUG - http-outgoing-4 >> Connection: Keep-Alive
DEBUG - http-outgoing-4 >> User-Agent: Apache-HttpClient/4.5.2 (Java/1.8.0_201)
DEBUG - http-outgoing-4 >> Accept-Encoding: gzip,deflate
DEBUG - http-outgoing-4 >> "GET /v1/commands/enrichmentmap/gseabuild?edbdir=C%3A%5CUsers%5CHatim%5Cgsea_home%5Coutput%5Cmar15%5Cmy_analysis.Gsea.1552686798457%5Cedb&pvalue=0.005&qvalue=0.1&overlap=0.5&similaritymetric=OVERLAP&combinedconstant=0.5&expressionfile=C%3A%5CUsers%5CHatim%5CDesktop%5CJunior%5Csem+2%5CKarras+Lab%5CGSEA%5CKarrasEXP.gct.txt HTTP/1.1[\r][\n]"
DEBUG - http-outgoing-4 >> "Host: localhost:1234[\r][\n]"
DEBUG - http-outgoing-4 >> "Connection: Keep-Alive[\r][\n]"
DEBUG - http-outgoing-4 >> "User-Agent: Apache-HttpClient/4.5.2 (Java/1.8.0_201)[\r][\n]"
DEBUG - http-outgoing-4 >> "Accept-Encoding: gzip,deflate[\r][\n]"
DEBUG - http-outgoing-4 >> "[\r][\n]"
DEBUG - http-outgoing-4 << "HTTP/1.1 500 Internal Server Error[\r][\n]"
DEBUG - http-outgoing-4 << "Date: Sun, 17 Mar 2019 06:55:04 GMT[\r][\n]"
DEBUG - http-outgoing-4 << "Content-Type: text/plain[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Origin: *[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Headers: origin, content-type, accept, authorization[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Credentials: true[\r][\n]"
DEBUG - http-outgoing-4 << "Access-Control-Allow-Methods: GET, POST, PUT, DELETE, OPTIONS, HEAD[\r][\n]"
DEBUG - http-outgoing-4 << "Content-Length: 20[\r][\n]"
DEBUG - http-outgoing-4 << "Server: Jetty(9.4.10.v20180503)[\r][\n]"
DEBUG - http-outgoing-4 << "[\r][\n]"
DEBUG - http-outgoing-4 << "Failed: empty String"
DEBUG - http-outgoing-4 << HTTP/1.1 500 Internal Server Error
DEBUG - http-outgoing-4 << Date: Sun, 17 Mar 2019 06:55:04 GMT
DEBUG - http-outgoing-4 << Content-Type: text/plain
DEBUG - http-outgoing-4 << Access-Control-Allow-Origin: *
DEBUG - http-outgoing-4 << Access-Control-Allow-Headers: origin, content-type, accept, authorization
DEBUG - http-outgoing-4 << Access-Control-Allow-Credentials: true
DEBUG - http-outgoing-4 << Access-Control-Allow-Methods: GET, POST, PUT, DELETE, OPTIONS, HEAD
DEBUG - http-outgoing-4 << Content-Length: 20
DEBUG - http-outgoing-4 << Server: Jetty(9.4.10.v20180503)
DEBUG - Connection can be kept alive indefinitely
INFO - status:Internal Server Error
INFO - Unable to create Enrichment Map:Internal Server ErrorResponseEntityProxy{[Content-Type: text/plain,Content-Length: 20,Chunked: false]}
DEBUG - Connection manager is shutting down
DEBUG - http-outgoing-4: Close connection
DEBUG - Connection manager shut down

Could you please explain the steps to reproduce this error? What data did you load, the settings and steps to connect to Cytoscape... and what do you expect from us. To me it seems like we should debug the GSEA software in order to solve this problem...

@ham123 Could you kindly review the answer below? I believe this should be definitely be upvoted and "accepted". The site member is knowledgeable in this specific problem.

This is a tough one, because the only relevant lines from your GSEA debug log were these two:
DEBUG - http-outgoing-4 << "Failed: empty String"
DEBUG - http-outgoing-4 << HTTP/1.1 500 Internal Server Error

My recommendations to you at the time would've been

try closing all the programs, restarting them, possibly rebooting and restarting them, and even analyzing a completely different set of data to see if you get the same results, —or—
update to the latest Java 8, reinstall Cytoscape 3.7.whatever and the EnrichmentMap plugin from scratch, and re-download GSEA, —or—
lastly, try reinstalling everything on a completely different machine/OS and see if you get different results (e.g., try on a Mac or Linux if you had problems with Windows originally, or vice-versa)

I know that's lame, but sometimes that's how you solve weird computer errors, when the error message the application itself gives you is non-specific and unhelpful.
If you had seen the error message
None of the gene sets have passed the filter. Try relaxing the gene set filter parameters.

in the GSEA application logs, then you would've had the problem described in EnrichmentMap's GitHub issue #358, and the solution would be to use less stringent p- and FDR Q-value cutoffs, as prescribed in the "Tips on Parameter Choice" section of the EnrichmentMap documentation.
My problem was that GSEA's default cutoffs (p-value < 0.005 and FDR Q-value < 0.1) were entirely too stringent for my results set, and using the "Very permissive" settings of p-value < 0.05 and FDR < 0.25 from the documentation solved it for me.
Hopefully this benefits someone who ends up here after googling for "Unable to create Enrichment Map:Internal Server ErrorResponseEntityProxy," which is just a stupendously unhelpful error message when "None of the gene sets have passed the filter" is actually what went wrong. Maybe go thumbs-up this comment in the EnrichmentMap issue tracker, if you think the maintainers need to work on that.


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

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


text:
Pipeline for extracting gene from multiple genomes for use in HyPhy selection analyses?
I have been trying to obtain some preliminary data from HyPhy selection analyses to inform a larger project. I have obtained a number of assembled mammalian genomes from NCBI with the initial goal of extracting from each of them their corresponding sequence for a specific gene. I wanted to use HyPhy to determine if there are signatures of positive selection in any of the mammalian lineages for that specific gene. However, I have been struggling to make headway with this for some time and am realizing that I may be out of my depth.
I was hoping that people could suggest, or at least point me in the direction of, a pipeline/guide/best practices for how one would go about achieving what I am trying to do. I would like to start from scratch and make sure I am performing my analyses in a proper way, using advice from other users.
As I mentioned, I currently have a collection of assembled genomes from NCBI which I would like to pull specific genes from. My understanding is that I will then need to extract the exons (I imagine there is a way to do this without having to extract the entire gene first), align the sequences from each species while preserving reading frames, and then use that alignment as input in HyPhy's codon models. It is very possible that I am missing out some fundamental steps in the process.
If anybody could offer direction/advice, I would be very appreciative.

Very interesting, personally I'd use PAML. @NatWH is good advice. Keep us posted.

That is interesting. I believe that others in my lab have more experience with HyPhy, hence why that is the goal at the moment. I will look into PAML as well. Thank you!

Okay well, because you're using HYPHY, you're going to need a gene tree of every gene you intend to analyze. This requires you to understand the homology relationships between the genomes at hand. Depending on your project, you might already have a list of homologs you're trying to analyze, or you can start from scratch.
If starting from scratch, you will need the proteomes of the species you're interested in (that is, the translation of the coding sequences of each gene, or to put it another way, the concatenated exons). Genome releases should typically have a proteome available for download, along with its corresponding coding sequence.
Next, you need some way to get the homolog groups for the genes in your genomes. You can use a pipeline like OrthoFinder, or iterative tree building and trimming, to infer homolog clusters (OrthoFinder calls them orthogroups but they're the same thing - homologous genes for the species of interest including both single copy orthologs and paralogs).

For this step, if you already have a set of genes and you know the copy numbers, or know that they're single copy orthologs, then it is more easy - you can just pull the sequence from each proteome and go from there.

Once you have your homologs, you need to align the amino acids using a tool like MAFFT, and then align the coding sequences to the amino acid alignment to attain a coding sequence alignment. For this you can use a tool like pal2nal, or check out pxaa2cdn from the phyx package written by some of my collaborators.

If the sequences aren't too divergent, and intron-exon structure is conserved, you could possibly also include alignments of the introns to provide more phylogenetic information. If you do use this approach however, you'd need to align the exons and introns separately. Also, the introns would not be included in the selection analysis.

Next, I would typically reinfer the gene tree of each gene using nucleotide models, which depending on the time, level of selective constraint, etc. can give more signal for the phylogenetic relationships. For this you can use a tree inference program like RAxML. With this tree and your coding sequence alignment, you'll be ready to run analyses in HYPHY (aBS-REL, MEME, FUBAR, etc., depending on your requirements).

Note that if you want, you can also do this in a Bayesian framework to integrate over uncertainty in the model parameters, tree topology, etc. Several Bayesian phylogenetic inference tools (e.g. MrBayes, PhyloBayes), allow you to infer putatively positively selected sites at the same time as tree inference.

Thank you, NatWH! Your answer seems perfect - comprehensive and easy to follow! I will try to follow this pipeline.

@PollardMD no problem, feel free to comment here if you run into any problems.


text:
Is this RNA seq data good based on the MDS plot
I am analyzing some RNA sequencing data from a collaborator where the effect of some ligands on the transcriptome is being looked at. I am using DESeq2 for my analysis. I am looking to compare each of the ligand treatments vs the DMSO. I have attached the MDS plot with this post. There seems to be a lot of variation between the 3 replicates across all the treatments. Based on this plot, is it worth it to use this data to run a differential expression?
It is also at this link below.
https://i.sstatic.net/ZXhk5.jpg


text:
Samtools index fails within nextflow workflow but not outside of it
I am a beginner to Nextflow  and using it to setup Chipseq workflow and in the workflow below the samtools_index process throws the following error:
[E::hts_idx_push] Chromosome blocks not continuous
  [E::sam_index] Read 'FCC1L20ACXX:4:1101:1645:2230#0/1' with ref_name='chr17', ref_length=83257441, flags=16, pos=58851543 cannot be indexed
  samtools index: failed to create index for "TCDD_plus_E2_ERa_Rep_1_ChIPSeq_first_1000.Aligned.sorted.bam": Inappropriate ioctl for device

The command executed is:
/Users/siddhaduio.no/Desktop/All_omics_tools/samtools-1.14/samtools index TCDD_plus_E2_ERa_Rep_1_ChIPSeq_first_1000.Aligned.sorted.bam

However, when I index the BAM file outside of nextflow, it runs without any problems. I was wondering what might be the issue here?
#!/usr/bin/env nextflow

nextflow.enable.dsl = 2

params.genome = '/Users/siddhaduio.no/Desktop/All_omics_tools/bowtie2-2.4.4-macos-x86_64/hg38_index/human_hg38_index.*.bt2'
params.reads = '/Volumes/matthews/AHR_ER_project_files/FASTQ/TCDD_plus_E2_ERa_Rep_{1,2,4}_ChIPSeq_first_1000.fastq'
params.results = '/Volumes/matthews/AHR_ER_project_files/FASTQ/results'

/*
* Process to align chip seq reads to genome with bowtie2
*/

process mapping {

publishDir '/Volumes/matthews/AHR_ER_project_files/FASTQ/results', mode: 'copy'

input:

path genomeDir

tuple val(replicateId), path(reads)

output:

tuple val(replicateId), path("${replicateId}.Aligned.sam")

script:

def btw2_index = genomeDir[0].simpleName

"""
/Users/siddhaduio.no/Desktop/All_omics_tools/bowtie2-2.4.4-macos-x86_64/bowtie2 -p 64 -q -x ${btw2_index} -U ${reads.join(',')} -S ${replicateId}.Aligned.sam 2> ${replicateId}.bowtie.log

"""

}

process convert_sam_to_bam {

publishDir '/Volumes/matthews/AHR_ER_project_files/FASTQ/results', mode: 'copy'

input:

tuple val(replicateId), path("${replicateId}.Aligned.sam")

output:

tuple val(replicateId), path("${replicateId}.Aligned.bam"), path("${replicateId}.Aligned.sorted.bam")

script:

"""

/Users/siddhaduio.no/Desktop/All_omics_tools/samtools-1.14/samtools view -S -b -o ${replicateId}.Aligned.bam ${replicateId}.Aligned.sam
/Users/siddhaduio.no/Desktop/All_omics_tools/samtools-1.14/samtools sort --output-fmt BAM -o ${replicateId}.Aligned.sorted.bam ${replicateId}.Aligned.bam
rm ${replicateId}.Aligned.sam

"""

}

process samtools_index {

publishDir '/Volumes/matthews/AHR_ER_project_files/FASTQ/results', mode: 'copy'

input:

tuple val(replicateId), path("${replicateId}.Aligned.sorted.bam")

output:

tuple val(replicateId), path("${replicateId}.Aligned.sorted.bam.bai")

script:

"""

/Users/siddhaduio.no/Desktop/All_omics_tools/samtools-1.14/samtools index ${replicateId}.Aligned.sorted.bam

"""

}

workflow{

reads_ch= Channel.fromPath(params.reads).map { fastq -> tuple(fastq.simpleName, fastq)}
index_ch= files(params.genome)
reads_ch
index_ch
mapping(index_ch,reads_ch)
convert_sam_to_bam(mapping.out)
samtools_index(convert_sam_to_bam.out)

}

```

You should have also received the following Nextflow warning:
WARN: Input tuple does not match tuple declaration in process `samtools_index`

Specifically, the problem is that the output of your convert_sam_to_bam process is:
tuple val(replicateId), path("${replicateId}.Aligned.bam"), path("${replicateId}.Aligned.sorted.bam")

But the samtools_index process is expecting:
tuple val(replicateId), path("${replicateId}.Aligned.sorted.bam")

If you would like to declare multiple output channels, you can use the emit keyword (this is documented under additional options) to define names for the output channel(s). For example:
process convert_sam_to_bam {

    input:
    tuple val(replicateId), path("${replicateId}.Aligned.sam")

    output:
    tuple val(replicateId), path("${replicateId}.Aligned.bam"), emit: aligned_bam
    tuple val(replicateId), path("${replicateId}.Aligned.sorted.bam"), emit: sorted_bam

    script:
    """
    samtools view \\
        -b \\
        -o "${replicateId}.Aligned.bam" \\
        "${replicateId}.Aligned.sam"

    samtools sort \\
        -o "${replicateId}.Aligned.sorted.bam" \\
        "${replicateId}.Aligned.bam"
    """
}

workflow {

    ...

    samtools_index( convert_sam_to_bam.out.sorted_bam )
}

Otherwise, either change the convert_sam_to_bam output declaration to:
tuple val(replicateId), path("${replicateId}.Aligned.sorted.bam")

Or use the map operator to get what you want, for example:
workflow {

    ...

    samtools_index(
        convert_sam_to_bam.out.map { rep, aligned, sorted -> tuple( rep, sorted ) }
    )
}

Note that a separate indexing process can be avoided using newer versions of samtools with the --write-index option:

To convert a SAM to a compressed BAM using BAI indexing: samtools view --write-index in.sam -o out.bam##idx##out.bam.bai

Thanks a lot! This was helpful.


text:
How do I start a nextflow run from last executed command
I am trying to use nextflow to process some ATAC seq data. The data is paired end and the files are huge. So, I first trimmed the adapters and now I want to align them to the genome.
My script is shown below. Now, I want to align them using bowtie2. However, when I run the script, the trimming is being rerun again instead of starting at the alignment part.
Below is my call to nextflow at the command line
nextflow run run_nextflow_ATACseq.nf -w /Volumes/FrailtData/ATAC_Seq_downloaded/results -resume

I was wondering what I am doing wrong such that the resume flag is not working?
Thanks for any suggestions!
#!/usr/bin/env nextflow

nextflow.enable.dsl = 2

params.reads = '/Volumes/FrailtData/ATAC_Seq_downloaded/PB_CD4plusCD45RAmin_D*_T72_*ADA_ATAC_{1,2}.fastq.gz'
params.index = '/Users/siddhaduio.no/Desktop/All_omics_tools/bowtie2-2.4.4-macos-x86_64/hg38_index/human_hg38_index.*.bt2'

/*
* Trim atac seq reads to align to genome with cutadapt
*/

process cut_adapt_trim {

publishDir '/Volumes/FrailtData/ATAC_Seq_downloaded/results', mode: 'copy'

input:

tuple val(replicateId), path(reads)

output:

tuple val(replicateId), path("${replicateId}.1.trimmed.fastq"), emit: trimmed_pair_1
tuple val(replicateId), path("${replicateId}.2.trimmed.fastq"), emit: trimmed_pair_2

script:
"""

cutadapt -a TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG -A GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG -o ${replicateId}.1.trimmed.fastq -p ${replicateId}.2.trimmed.fastq ${reads}

"""

}

/*
* Align trimmed atac seq reads to genome
*/

process align_to_genome {

publishDir '/Volumes/FrailtData/ATAC_Seq_downloaded/results', mode: 'copy'

input:

path index

tuple val(replicateId), path("/Volumes/FrailtData/ATAC_Seq_downloaded/results/${replicateId}.1.trimmed.fastq")
tuple val(replicateId), path("/Volumes/FrailtData/ATAC_Seq_downloaded/results/${replicateId}.2.trimmed.fastq")

output:

tuple val(replicateId), path("${replicateId}.sam"), emit:aligned_sam
tuple val(replicateId), path("${replicateId}.bam"), emit:aligned_bam

script:

"""
def bwt_index = index[0].simpleName

/Users/siddhaduio.no/Desktop/All_omics_tools/bowtie2-2.4.4-macos-x86_64/bowtie2 -p 64 -q -x ${bwt_index} -1 ${replicateId}.1.trimmed.fastq -2 ${replicateId}.2.trimmed.fastq -S ${replicateId}.Aligned.sam 2> ${replicateId}.bowtie2.log | /Users/siddhaduio.no/Desktop/All_omics_tools/samtools-1.14/samtools view -S -b - > ${replicateId}.bam

"""

}

workflow{

reads_ch= Channel.fromFilePairs(params.reads)
index_ch= Channel.fromPath(params.index)
Channel.fromFilePairs(params.reads).view()
index_ch.view()
cut_adapt_trim(reads_ch)
align_to_genome(index_ch,cut_adapt_trim.out)

}
```

You can achieve your desired checkpointing with some if statements and a user defined flag pretty easily. Something like:
#!/usr/bin/env nextflow
nextflow.enable.dsl = 2

params.skip_trim = null
params.reads = '/Volumes/FrailtData/ATAC_Seq_downloaded/PB_CD4plusCD45RAmin_D*_T72_*ADA_ATAC_{1,2}.fastq.gz'
params.index = '/Users/siddhaduio.no/Desktop/All_omics_tools/bowtie2-2.4.4-macos-x86_64/hg38_index/human_hg38_index.*.bt2'
params.trimmed_dir='01_trimmed'

process cut_adapt_trim {
   publishDir ${params.trimmed_dir}, mode: 'copy'

   input:
   tuple val(replicateId), path(reads)

   output:
   tuple val(replicateId), path("${replicateId}.1.trimmed.fastq"), emit: trimmed_pair_1
tuple val(replicateId), path("${replicateId}.2.trimmed.fastq"), emit: trimmed_pair_2

   script:
   """
   cutadapt -a TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG -A GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG -o ${replicateId}.1.trimmed.fastq -p ${replicateId}.2.trimmed.fastq ${reads}
   """
}

process align_to_genome {
   publishDir '/Volumes/FrailtData/ATAC_Seq_downloaded/results', mode: 'copy'

   input:
   path index
   tuple val(replicateId), path("/Volumes/FrailtData/ATAC_Seq_downloaded/results/${replicateId}.1.trimmed.fastq")
tuple val(replicateId), path("/Volumes/FrailtData/ATAC_Seq_downloaded/results/${replicateId}.2.trimmed.fastq")

   output:
   tuple val(replicateId), path("${replicateId}.sam"), emit:aligned_sam
tuple val(replicateId), path("${replicateId}.bam"), emit:aligned_bam

   script:
   """
   def bwt_index = index[0].simpleName

   /Users/siddhaduio.no/Desktop/All_omics_tools/bowtie2-2.4.4-macos-x86_64/bowtie2 -p 64 -q -x ${bwt_index} -1 ${replicateId}.1.trimmed.fastq -2 ${replicateId}.2.trimmed.fastq -S ${replicateId}.Aligned.sam 2> ${replicateId}.bowtie2.log | /Users/siddhaduio.no/Desktop/All_omics_tools/samtools-1.14/samtools view -S -b - > ${replicateId}.bam

   """
}

workflow{
   index_ch= Channel.fromPath(params.index)
   if( params.skip_trim == false ) {
       reads_ch= Channel.fromFilePairs(params.reads)
       cut_adapt_trim(reads_ch)
       align_to_genome(index_ch,cut_adapt_trim.out)
   } else if( params.skip_trim ) {
       trim_ch= Channel.fromFilePairs("${params.trimmed_dir}/*_{1,2}.trimmed.fastq")
       align_to_genome(index_ch,trim_ch)
   }
}

Then to execute the skip you would use:
nextflow run run_nextflow_ATACseq.nf -w /Volumes/FrailtData/ATAC_Seq_downloaded/results --skip_trim

Play around with the trimmed_dir variable, I don't know what you are mounting to. But using variables for a lot of the paths and options into the params section makes it easier to adjust in the future for other files.
There is another point I would like to address. Nextflow emits and inputs files in random order. So by splitting the R1 and R2 reads into separate channels, you are going to have a lot of mixed samples as your output. Make the output a single tuple: tuple val(id), path(r1), path(r2), emit: trimmed_reads

Thank you again for the feedback. Very helpful!

Your task cache is being invalidated because your process inputs are non-deterministic. Like in the example, the solution is to join the two output channels before process invocation. Another way is to simply output the trimmed reads in a single channel, for example:
process cut_adapt_trim {

    tag "$replicateId"

    publishDir '/Volumes/FrailtData/ATAC_Seq_downloaded/results', mode: 'copy'

    conda 'cutadapt'

    input:
    tuple val(replicateId), path(reads)

    output:
    tuple val(replicateId), path("${replicateId}.1.trimmed.fastq"), path("${replicateId}.2.trimmed.fastq")

    script:
    """
    cutadapt \\
        -a TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG \\
        -A GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG \\
        -o "${replicateId}.1.trimmed.fastq" \\
        -p "${replicateId}.2.trimmed.fastq" \\
        ${reads}
    """
}

process align_to_genome {

    tag "$replicateId"

    publishDir '/Volumes/FrailtData/ATAC_Seq_downloaded/results', mode: 'copy'

    conda 'bowtie2'

    input:
    tuple val(replicateId), path(fastq_1), path(fastq_2)
    path index

    output:
    tuple val(replicateId), path("${replicateId}.bam")

    script:
    def bwt_index = index[0].simpleName

    """
    bowtie2 \\
        -p 64 \\
        -q \\
        -x "${bwt_index}" \\
        -1 "${fastq_1}" \\
        -2 "${fastq_2}" |
    samtools view \\
        -b \\
        -o "${replicateId}.bam" \\
        -
    """
}

workflow{

    reads_ch = Channel.fromFilePairs( params.reads )
    index_ch = Channel.fromPath( params.index )

    cut_adapt_trim( reads_ch )

    align_to_genome( cut_adapt_trim.out, index_ch )

    align_to_genome.out.view()

}

Run using:
nextflow run run_nextflow_ATACseq.nf -with-conda -resume


text:
Struggling to install scvi-tools on Mac M1
I was wondering if somebody here has experience with installing scvi-tools on a Mac with M1 chip. I tried the following command:
conda install scvi-tools -c conda-forge

But got the following error:
Collecting package metadata (current_repodata.json): failed
CondaSSLError: OpenSSL appears to be unavailable on this machine. OpenSSL is required to
download and install packages.
Exception: HTTPSConnectionPool(host='conda.anaconda.org', port=443): Max retries exceeded with url: /conda-forge/osx-64/current_repodata.json (Caused by SSLError("Can't connect to HTTPS URL because the SSL module is not available."))

Then I tried to install OpenSSL with following command:
conda install -c anaconda openssl

and I got the exact same error as the conda install scvi-tools -c conda-forge command.
I do not understand why collection of metadata is failing for openssl. Any inputs would be highly helpful.
Thanks!

There might be system incompatibilities, thats what my error log suggests. Thus I was unable to easily install as well. However this might be a Python version  issue (could be Python 3.6, I used 3.11 [not the best idea]). However the error log should have flagged that. scvi-tools is a bioinformatics tool, but please try and edit the question with the biological problem. I'd guess this is RNA-seq PCA type analysis

Wait, how did you end up with a non-working conda install on that M1 mac in the first place?  Looks like it's not complaining about a missing package dependency; it's complaining about not having openssl to do HTTPS downloads within conda.

Did you try:
pip install scvi-tools

?

Would need venv (virtual env) with pip, but I agree I think this is the best approach in this particular scenario

https://scvi-tools.org/


text:
Error in FindIntegrationAnchors- Seurat package
I am working on integrating a labelled single cell RNA seq cell atlas with an unlabelled one. I am wondering how do I determine the value of max.features when integrating two single cell RNA seq datasets?
Below is the code I am using:
epithelial = Read10X(data.dir="/projects/b1025/sdi0596/Covid19_Single_Cell_RNASeq_project/filtered_matrices_count/epithelial/filtered_feature_bc_matrix/",gene.column=1,cell.column=1)
habermann = Read10X(data.dir="/projects/b1025/sdi0596/Covid19_Single_Cell_RNASeq_project/filtered_matrices_count/Habermann_Single_Cell_Data/",gene.column=1,cell.column=1)

data.list=c(epithelial,habermann)

normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = data.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

Perform integration
immune.achors= FindIntegrationAnchors(object.list = ifnb.list,max.features = 200,verbose = TRUE)

Error in FindIntegrationAnchors(object.list = ifnb.list, max.features = 200, :
Max dimension too large: objects 1, 2 contain fewer than 30 cells.
Please specify a maximum dimensions that is less than the number of cells in any object (4).

Regarding the dimension of these objects, the epithelial one is 33538 x 49109 and the habermann is 33694 x 220213.
Any insights would be highly helpful.
Thank you.

There's a few problems with your code, first, when you do Read10X() it returns you a sparse matrix, and you need to put this into a Seurat object with meta data, before doing the integration.
So for example, i create some example data that is similar to your output from Read10X() :
library(Matrix)
epithelial = Matrix(matrix(rnbinom(100000,mu=200,size=1),ncol=50),sparse=TRUE)
rownames(epithelial) = paste0("gene",1:2000)
colnames(epithelial) = paste0("e",1:50)

habermann = epithelial + rnbinom(100000,mu=20,size=1)
colnames(habermann) = paste0("h",1:50)

Then convert to a seurat object:
data.list=c(epithelial = CreateSeuratObject(epithelial),
            habermann = CreateSeuratObject(habermann))

Then find variable features:
data.list <- lapply(data.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

Should look like this:
$epithelial
An object of class Seurat
2000 features across 50 samples within 1 assay
Active assay: RNA (2000 features, 2000 variable features)

$habermann
An object of class Seurat
2000 features across 50 samples within 1 assay
Active assay: RNA (2000 features, 2000 variable features)

Now perform integration, below I have to reduce k.filter because I have very little cells in this example. In your case, you can simply use the default settings.
features <- SelectIntegrationFeatures(object.list = data.list)
anchors= FindIntegrationAnchors( data.list,max.features = 200,
k.filter=50,k.anchor = 3,verbose = TRUE)

It is possible that some of your samples have too few cells, so it is recommended to check the data and remove samples with low cell counts.

As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center.


text:
Clustering by individual tissue population post integration of two datasets
Sorry if this question has been asked. Post batch correction and integration of two datasets, where cells of multiple different tissue types are present, I want to run clustering on cells of corresponding tissue types between the two datasets. Should I simply subset the integrated object to get each tissue type, scale the data and run PCA and then run clustering on it?
Thanks so much for the clarifications.


text:
how do I convert a vcf file for a GWAS study (using R package vcfR)
I have a vcf file with individual level genetic data, that I read using the R package vcfR
gen <- read.vcfR('test.vcf')

Scanning file to determine attributes.
File attributes:
  meta lines: 9
  header_line: 10
  variant count: 152
  column count: 3798
Meta line 9 read in.
All meta lines processed.
gt matrix initialized.
Character matrix gt created.
  Character matrix gt rows: 152
  Character matrix gt cols: 3798
  skip: 0
  nrows: 153
  row_num: 0
Processed variant: 153
All variants processed

How do I convert this file to a data frame that I can use for a GWAS or Mendelian Randomization study?
I was told that this should be a data frame where each column is a genetic variant and it take values 0, 1, or 2, depending on the number of "mutations".
I never worked with genetic data before and do not know much about genetics.

Welcome to Bioinformatics.StackExchange! GWAS is a broad area, and there are numerous tools you could use. The format of your data will depend on which tool you are using, and hopefully this will be clearly documented in the tool's manual. If on the other hand you're trying to implement your own tool, you'll have to be much more specific about what you want (or are required) to accomplish.

Welcome to the site! Could you please explain what have you tried? For instance, what happens if you use as(gen, "data.frame")?

I found that I can use extract.gt, then transform the matrix. Entries in this matrix are "0/0", "0/1", "1/0", and "1/1". I now have to recode these to 0, 1, 1, 2.

If you want to do GWAS or mendelian randomisation, you can do it with Plink (v2) which should be faster than R.


text:
What does this genome file from NCBI means?
I obtained a genome from NCBI. For instance, Haloquadratum walsbyi strain AM180088.1. But the FASTA file contain two parts. The first part marked as
>AM180088.1 Haloquadratum walsbyi DSM 16790 complete genome
The second part marked as
>AM180089.1 Haloquadratum walsbyi DSM 16790 plasmid PL47 complete genome
What does that mean? In particular, the second one "plasmid"? If it is called complete genome should it all in one? Why divide into parts?

Archea can have plasmids last time I checked...

what does that mean? So, plasmid have different gene from other part of archaea?

Yes, bacteria and sometimes both archea and eukaryotes can have plasmids (google the term if you don't know it) in addition to their normal genome.

I know they have plasmid. But will the plasmid have different gene from the bacteria/archaea itself?

Typically, yes.

I'm voting to close this question as off-topic and not worth migrating to the biology stackexchange.


text:
The best tool to count the number of insertion sequences in a set of bacterial genome
I have a set of 414 GenBank genomes and I need to count a number of insertion sequences (ISs, transposases) in each of them. Is there any traditional way to perform such an analysis?
I came up with three possible ways:
1) Count a number of "transposase" as a product of CDS for all CDS in annotation file of each genome. -- it depends on annotation
2) Use a tool named OASIS. -- it depends on annotation
3) Use a web-server named ISsaga. -- it can analyse only 10 replicons at once
4) Perform blast(n or p) against ISfinder database and count a number of hits with low e-value. -- it seems not reliable
Any help is appreciated!


text:
Fasterq-dump: --split-spot or -concatenate-reads?
After using files that I downloaded from the SRA with fasterq-dump, I realize I am not 100% sure that I have all the data.
I noticed in my downstream analysis that I seem to be missing the .1 and .2 numbers in the code associated with individual reads with the same spot ID. Then I noticed that my initial fastq file doesn't seem to have those either.
So I'm wondering if the -split-spot option that I used is doing what I want it do do, which is 'get all the sequence data', and how does it differ from --concatenate-reads. I am afraid I did not understand this from the program description. Please help?

Hi Laura, do you have paired-end reads? Can you provide the SRA accession too?

@StupidWolf It looks like they're paired. Eg. SRR6462984.63.1 and SRR6462984.63.2 .

Fasterq comes from the latest version of sratools. So if you check the manual , it says the equivalence is:

fastq-dump SRRXXXXXX --split-3 --skip-technical
fasterq-dump SRRXXXXXX

In older versions of sratoolkit, if you use fastq-dump without specifying --split-3 for paired-end reads, you get the format mentioned, spotID.1 for forward, spotID.2 for reverse:
fastq-dump SRR6462984.sra
more SRR6462984.fastq

You get:
@SRR6462984.1 1 length=301
NATCTGCCCGTTCCACATAAACACCGTTAAACATCGATAGAGCGAAATAAAGCTGCTGAGTGATGTGATGATAAAGCTTTTCCGTTCTGACACGATCTTCCATCTCGTCTATCATTGCATCGAGAGCGTCAGAGATCGCAAGCAGTGCGGATATTGCCGTCCGTCTTGAAGCTTTGGCAGAGCCAGGCGGAATTTGTTTATCGGATAGTGTTTATGCGCAAATTTAGCATATGATTTCCGACAACTTTTGCTTCATCGGGCCGCAAAGCGTGAAAAACAGTGTCGTCACCGTTGATGTGTG
+SRR6462984.1 1 length=301
#/AAAEE/EEEA/EAAEEAEE/EE/EEA/EEEEE/AAE/EE/EAEAE//EEAEE/EEEEE<EEEEE/EE<EE/AEEEE/EAEEEEA/E<EEAEAEEEEEEAEEEE<EEAEEEEEAAAEEE/EEEEE/E//AEEEE6EA66EE<<<EAAA/AAAAAE/E6E/EA/AE/AEEEEEEAAEEEEA6EEEEEEEEE/EEE/AEEAAAEEEAE/E/EEEEEA//E//EEE/EEE/A<//AE/EE//6A///AA/E//E/E/EE//AE//AAEEE////E/A<AAE/EEAE//A//6/A</<6/</E/
@SRR6462984.2 2 length=302
NTCAGCAGGCGCCGTAACTTCAGCAAAACTGTAGCTATAACTTTTCCGAACTTTATACCTAACTTCGCTTTCCTAGCACTGTATAAAATGCCTAAAAGAAATCAGCAGGCGCCGTCACTTCAAAAAATCAGTATGGATCAATCAGACGCGGAGCTTATTAAAACAGAGCTTATTAAAAATTATGTCTTTCTTCTTCTTTCAGCAATTCGCAACGTAAAATTCAAAATTTTGGACGCCGTAATTGCGCAGGAACCAATTCAAAGAGGCGGCGGCCAATTAATTCCAACGTCTCGAAGAGTGGC
+SRR6462984.2 2 length=302

With the latest version, it by defaults splits the paired end sra file into forward and reverse fastq (i.e --split-3):
fasterq-dump SRR6462984.sra
head -2 SRR6462984_*.fastq

And you get
==> SRR6462984_1.fastq <==
@SRR6462984.1 1 length=150
NATCTGCCCGTTCCACATAAACACCGTTAAACATCGATAGAGCGAAATAAAGCTGCTGAGTGATGTGATGATAAAGCTTTTCCGTTCTGACACGATCTTCCATCTCGTCTATCATTGCATCGAGAGCGTCAGAGATCGCAAGCAGTGCGG

==> SRR6462984_2.fastq <==
@SRR6462984.1 1 length=151
ATATTGCCGTCCGTCTTGAAGCTTTGGCAGAGCCAGGCGGAATTTGTTTATCGGATAGTGTTTATGCGCAAATTTAGCATATGATTTCCGACAACTTTTGCTTCATCGGGCCGCAAAGCGTGAAAAACAGTGTCGTCACCGTTGATGTGTG

So if you need to get back the .1 and .2 format, use fastq-dump from the latest installation.
I would suggest sticking to the newest format, with 2 files, because fastq-dump will be deprecated in the future.

Can you have a look at this https://bioinformatics.stackexchange.com/questions/21407/multiple-srr-id-downloaded-for-single-gsm-ids-input-issue


text:
Can I get longer contigs by changing MEGAHIT settings?
My current settings are megahit -r input.fq --num-cpu-threads 32 --min-contig-len 300 --presets meta-large -o output. I picked the 'presets meta-large' because I am dealing with a quite large and diverse metagenome out of which I only need some genes.
I've got a contig of interest ~19kb, with a depth of 16. Could I make this contig larger with more relaxed settings? Help would be appreciated!

You can use the option --k-step 10 and --kmin-1pass for low-coverage datasets.