text:
Finding annotated counterpart after BLASTn with efetch (Biopython)
I am creating a pipeline for the identification of unknown transcripts. After a local BLASTn search of the transcripts, I have a large list of the respective hits with different genomes. I have the accession code of the respective genome, start/end coordinates, and the strand where the hit is found. Now, I would like to find if the genome where the hit was found is annotated within the coordinates returned by BLAST.
I have found the possibility of using efetch to search for each genome via Biopython, but this only returns the whole annotation, even though I limit the search to the specific coordinates with:
Entrez.efetch(db="nuccore", id="CP054847.1", strand="2", seq_start="326053", seq_end="326786", retmode="xml")

Of course, I could parse the XML and look for possible annotations within the given range, but not only would this be quite computationally intense since I expect around 10.000+ hits to be searched, but it also seems to return the same results for the plus/minus strand, even though the annotations should be strand-specific.
I have also intended to use the standalone tool of efetch, but it would be easier to integrate the search to the pipeline using Biopython.
Do you have any ideas why Entrez.efetch is returning the whole xml or is there any fast way of filtering the output? Do you have maybe another alternative/idea for this step? I would appreciate any suggestion!

can you find any difference between https://www.ncbi.nlm.nih.gov/nuccore/CP054847.1/ and https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP054847.1 database entries ?

It seems like the seq_start and seq_end are sensitive to the direction. Try swapping the two.

Lets assume you can get the blastn output in a BED format:
genomeA_contig001    hit_start    hit_end    query_transcript    score   strand
genomeB_contig021    hit_start    hit_end    query_transcript    score   strand

which will contain all the unknown hits to all genomes.
All you need is a set of GTF files specific for genomeA, genomeB and so on where you add the suitable prefix (genomeA_ for GTF describing genes/exons in that genome and so on). For transcript hits you may start with just the annotated exons. Next run:
bedtools intersect -wao -a all.transctipts.hits.bed -b all_annotated_exons.gtf my_result.tsv
to get all the overlaps with number of overlapping bases.
Filter and sort the my_result.tsv possibly importing it into a data frame using pandas or polars.


text:
Annotating gene names or gene IDs to a dataframe containing SNPs?
I have a large data-frame (excel file) of SNPs with genotyping data. I need to filter the data to get SNP information of a specific gene alone. The list is too long to do it manually. I was wondering of annotating the SNPs with the gene names so that I can filter using the gene name. Does anyone know of a good R package to perform the task?

I assume your SNP file contains genomic positions - please clarify if not. For example, you could have an rsID column, chromsome, and position column that contain "rs3088379", "chr1", "54225281".
Approach #1: use R, as you suggest.
ensembldb is an R package that may be what you're looking for: https://bioconductor.org/packages/release/bioc/html/ensembldb.html.
You could map the genomic positions of your SNPs to gene names by adapting their workflow: https://bioconductor.org/packages/release/bioc/vignettes/ensembldb/inst/doc/coordinate-mapping.html.
After you've mapped the genomic position, you could filter for your gene name.
Approach #2: use bedtools.
You can convert your Excel file of genomic positions to a tab-delimited BED format. https://bedtools.readthedocs.io/en/latest/content/overview.html has an overview. So, your first SNP would turn into
chr1 54225280 54225281 rs3088379

--> make sure you turn 1-based to 0-based coordinates!
You could then download the coordinates of your gene in BED format from your favorite source, like UCSC Table Browser.
Then, you run bedtools intersect with your SNP BED and your gene BED. What remains is just that gene's SNPs.

source: Stack Exchange

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

sort: votes


text:
Filter rows of VCF file for Match=EXACT?
How do I remove rows in a VCF file on 1kGenome column, where Match=EXACT using bash ?

Please edit your question to include a couple of relevant lines from your .vcf.

You can filter records in Python for example using pyvcf.VcfFrame.filter_flagany() method I wrote:
>>> from fuc import pyvcf
>>> data = {
...     'CHROM': ['chr1', 'chr1', 'chr1', 'chr1'],
...     'POS': [100, 101, 102, 103],
...     'ID': ['.', '.', '.', '.'],
...     'REF': ['G', 'T', 'A', 'C'],
...     'ALT': ['A', 'C', 'T', 'A'],
...     'QUAL': ['.', '.', '.', '.'],
...     'FILTER': ['.', '.', '.', '.'],
...     'INFO': ['DB', 'DB;H2;Match=EXACT', 'DB;H2;Match=EXACT', '.'],
...     'FORMAT': ['GT', 'GT', 'GT', 'GT'],
...     'Steven': ['0/0', '0/1', '0/1', '0/0'],
... }
>>> vf = pyvcf.VcfFrame.from_dict([], data)
>>> # vf = pyvcf.VcfFrame.from_file('your_vcf.vcf')
>>> vf.df
  CHROM  POS ID REF ALT QUAL FILTER               INFO FORMAT Steven
0  chr1  100  .   G   A    .      .                 DB     GT    0/0
1  chr1  101  .   T   C    .      .  DB;H2;Match=EXACT     GT    0/1
2  chr1  102  .   A   T    .      .  DB;H2;Match=EXACT     GT    0/1
3  chr1  103  .   C   A    .      .                  .     GT    0/0
>>> vf.filter_flagany(['Match=EXACT'], opposite=True).df
  CHROM  POS ID REF ALT QUAL FILTER INFO FORMAT Steven
0  chr1  100  .   G   A    .      .   DB     GT    0/0
1  chr1  103  .   C   A    .      .    .     GT    0/0
>>> # vf.to_file('filtered.vcf')

Hard to say without seeing VCF rows, but generally speaking with text files you can use grep, if I understand the question:
grep -vw "Match=EXACT" file.vcf > filtered.vcf

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


text:
How to find genes whose expression vary the most in time-series?
Hi I am reading this paper and trying to analyse the same datasets in my work:
dynGENIE3: dynamical GENIE3 for the inference of gene networks from time series expression data
In this paper, the authors are working on several time course gene expression datasets, including "Drosophila melanogaster embryogenesis" dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6186) . This dataset has the whole genome ( 14064 genes) of Drosophila melanogaster measured over 31 time-points.
In the dynGenie3 paper, we analyse a subset of these  14064 genes, as they write:
"We focused our analysis on the 1000 genes whose expression vary the most across the time series."
I want to know how do we measure the variance in expression across time series? What statistical tests might have been performed? I couldn't find any more information in the paper so I am assuming it might be something simple and frequently done.  Need some hints.

https://hbctraining.github.io/DGE_workshop_salmon_online/lessons/08b_time_course_analyses.html

Maybe my thoughts are too simple, but wouldn't the computation of the variance per gene suffice for this? You really do not need a statistical test, since they are not talking about DEGs, just high-variant genes.

I think so too, so I just calculate the variance of each gene and pick the top 1000? The expression values are normalised. Variance will just tell how much the expression of genes differ from the mean?

That would be my way to go: Get the genes that spread the most and label them as "most variant". A non-variant gene would just be a flat profile, hence the distance to the mean for each time point would be quite low.

Thank you. I did that now and I hope it is consistent with that they might have done.

Hope it helps! Keeps me posted if you find that they used another method.

Actually i contacted the authors and this is what they did. Thanks a lot.


text:
How to know the number of SNP in my samples after running the population module in stacks
I am new to stacks for analyzing RAD-Seq data, I have run all the modules. I will like to know to get the total SNP from my sample after running the population module. Is there a command line to use and inside which file?
command:
populations -P -M (denovo output) (popmap) -t 5 --vcf --fasta-samples --fasta-loci -O (pop-output) -r 0.8 -p 5

Thanks

Update
Based on the newly provided command (which has --vcf), there should be a VCF file with suffix .vcf in the output directory.
Original answer
It's pretty hard to answer without more information. For example, it's helpful to link to the software and the module that you are working from (I had never heard of Stacks). What do you mean by "total SNP"? Also, what is the actual command that you ran? There are a lot of possible options!
But (based only on reading the docs for the command), it looks like there are a variety of formats in which the SNP data can be output:
File output options:
--ordered-export — if data is reference aligned, exports will be ordered; only a single representative of each overlapping site.
--fasta-loci — output locus consensus sequences in FASTA format..
--fasta-samples — output the sequences of the two haplotypes of each (diploid) sample, for each locus, in FASTA format.
--vcf — output SNPs and haplotypes in Variant Call Format (VCF).
--genepop — output results in GenePop format.
--structure — output results in Structure format.
--radpainter — output results in fineRADstructure/RADpainter format.
--plink — output genotypes in PLINK format.
--hzar — output genotypes in Hybrid Zone Analysis using R (HZAR) format.
--phylip — output nucleotides that are fixed-within, and variant among populations in Phylip format for phylogenetic tree construction.
--phylip-var — include variable sites in the phylip output encoded using IUPAC notation.
--phylip-var-all — include all sequence as well as variable sites in the phylip output encoded using IUPAC notation.
--treemix — output SNPs in a format useable for the TreeMix program (Pickrell and Pritchard).

I would suggest examining the command that you used, specifically what you passed as the output format. Here are the example commands from the docs:
populations -P dir -b batch_id [-O dir] [-M popmap] (filters) [--fstats] [-k [--window_size=150000] [--bootstrap [-N 100]]] (output formats)
populations -V vcf -O dir [-M popmap] (filters) [--fstats] [-k [--window_size=150000] [--bootstrap [-N 100]]] (output formats)

Can you track down what you specified as output? That will contain the SNP data / genotypes.
Possibly you will need to rerun the command while specifying a different, preferred format.


text:
Demultiplexing FASTQ file without index information
I am trying to understand how data that was uploaded to SRA (https://www.ncbi.nlm.nih.gov/sra?LinkName=biosample_sra&from_uid=4510743) can be analyzed with the assumption that the FASTQ file should contain replicate information (3 replicates within the sample). I downloaded the file and ran it with Nextflow, and it was mapped to the human genome and determined to contain unstranded reads of 50bp, both directions.
Which steps can I take to infer if the FASTQ file indeed is not demultiplexed?
Many thanks,
Geo
The header of the FASTQ looks like this:
@SRR3467208.1 1/1
NTCTCTATGTCCACTCTGGAGCCTTTAAGTGCCACTTGAGGGCCTTTAAC
+
#4:D;DDEHHFHHIJJIJJGHGHJJI@HFIIJJHJJIG?9CFHIIIHHHG
@SRR3467208.2 2/1
NGTGCTGTGTGTGCATGTGTGTGCGTGTGTGTGCTGTGCGTTTGTGTGTG
+
#1:B:BDDFFFFFFFBGFHIFIIFGFFGFGFGIIIIIIEFGGIIFFIFII
@SRR3467208.3 3/1
NGTCACTTGGCATCTGGCCATCGGGCTGGATGCCGTGTTCCAGGCAGTAG
+
#1:A=DDDB<<CFAFGI@FGG9A@68CEG;DB?<DGGGE49?BBF2B=FF
@SRR3467208.4 4/1
NTCCGGGGCTCCAGCAACCAGAAGAGGGAAAAGCTGTCTTCGGTGATGTC
+

I would suggest looking at sample metadata:

Library:

Name: TruSeq RNA

Instrument: Illumina HiSeq 2000

Strategy: RNA-Seq

Source: TRANSCRIPTOMIC

Selection: PolyA

Layout: PAIRED

You can see that it's paired. You can look at the spot descriptor at the page you link to see more information about the design.
There is no multiplexing information available there. I would instead look at the overall project and its metadata to find the different replicates. It looks like there are 2 samples x 3 replicates each, deposited independently. So it is likely to already be demultiplexed, and you can just download them individually.
"Unstranded" isn't really a meaningful term for this, if I understand correctly. Forward and reverse read pairs don't correspond to forward and reverse on the chromosome, but to the way that sequencing-by-synthesis happens on the machine. So the mapping orientation is not really relevant.
I infer that you are working from a single FASTQ file. If you are working from a single FASTQ file for paired-end out of SRA, then the FASTQ file is not yet resolved into F+R reads. I would suggest rerunning the FASTQ conversion with --split-3 or similar option.

Thanks! The term 'unstranded' was determined by RSeQC infer_experiment.py. The reads were 47% forward, and 47% reversely mapped. So the resulting counts from the single FASTQ file worked even w/o the FASTQ conversion. The experiments from the entire project are all 6 independent (6 x 1 sample), not 2 x 3 samples. I worry that the reported replicates are actually missing in the uploads and was hoping to confirm this from the FASTQ files somehow.

@GeoVogler That is helpful background for trying to answer the question. Where do you see the 6 x 1 info? I was not able to track down the actual study that originated these samples, linking that would be helpful. All I am seeing is the 3 samples x 2 treatments on SRA.

The three samples are from three individuals (mother/father/proband), not a single one (see Biosamples). The layout of SRA with Project/BioSamples/Reads could really be improved. If there is no indication of replicates I will need to contact the lab directly, I just thought maybe there is an easy way to gauge the presence of Index sequences inside a FASTQ file.

@GeoVogler the index reads were probably split out in initial processing, but they are often preserved in FASTQ headers by e.g. bcl2fastq. SRA unfortunately has a tendency to strip out useful information from the FASTQ header, which does sometimes include the index sequence for each read. The authors should have the original files in theory, yes.


text:
Challenging benchmarks for supervised learning on sparse scRNA-seq data
One challenging aspect of modeling scRNA-seq data is data sparsity, that is, scRNA-seq measurements typically suffer from large fractions of observed zeros (i.e. dropouts), where a given gene in a given cell has no unique molecular identifiers or reads mapping to it [1].
What would be an appropriate dataset to benchmark supervised learning algorithms that operate on single-cell RNA-seq data? In particular, I am looking for challenging datasets (preferably with labels, i.e. cell type classification) where existing methods struggle due to the sparsity of the data. Additionally, does it ever happen that the dropouts are confounders of the labels? (i.e. for a given cell type A the dropout probability of a gene X is higher than for another cell type B).
Any pointers are welcome.

Can you say what is the purpose of your supervised learning method? to annotate cell types?

Yes, that's right, cell type annotation is the first thing that comes to mind, but I am also willing to consider other tasks (e.g. classification of cancer cells). Even regression tasks would be fine too. Essentially, I have a method in mind that might be useful for modeling non-linear+sparse data and I am looking for an appropriate benchmark dataset.

You could try the PBMC 3k data from the Satija lab:
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
For well-annotated data, there's the Single-cell proteo-genomic reference map:
https://www.biorxiv.org/content/10.1101/2021.03.18.435922v1
Data for that are available here.
If these datasets are too large (i.e. too good), you can use the existing counts as weighting for a random removal of reads. Genes with high expression are likely to be present regardless of sequencing depth, hence the count-adjusted removal.
FWIW, I don't think that supervised learning is a good fit for dealing with data sparsity in cDNA read counts. There are plenty of existing attempts to solve the issues by adjusting the assumed models, including research by the authors of commonly-used bulk differential expression tools (e.g. here).

Thanks! I like the second reference, it is fresh from the oven. If the data is "too good" (i.e. traditional models can easily distinguish the different cell types), you would then introduce random dropouts to increase the sparsity?

Random dropouts at the read level; I wouldn't randomly drop out genes.


text:
Nextflow: Runs without error, produce few but not all output files into output channel
I was trying to develop a pipeline using nextflow and plink. My code runs without any error, but got only two output files instead of 22 (one for each chromosome). Here is part my code block, any help?
process step5 {
publishDir params.publish, mode: 'copy'
input:
val(hrc) from HRC_1000G_check_bim
val(hrc_site) from hrc_sites_file
tuple val(prefix), val(chr_name), path(bfiles) from step4_frq
output:
tuple val("${prefix}"), val(chr_name), path("LOG-${prefix}${chr_name}.step4-HRC.txt") into step4_log
tuple val("${prefix}"), val(chr_name), path("Exclude-${prefix}${chr_name}.step4-HRC.txt") into step4_exclude
tuple val("${prefix}"), val(chr_name), path("Chromosome-${prefix}${chr_name}.step4-HRC.txt") into step4_chromosome
tuple val("${prefix}"), val(chr_name), path("Position-${prefix}${chr_name}.step4-HRC.txt") into step4_position
tuple val("${prefix}"), val(chr_name), path("Strand-Flip-${prefix}${chr_name}.step4-HRC.txt") into step4_strand
tuple val("${prefix}"), val(chr_name), path("Force-Allele1-${prefix}${chr_name}.step4-HRC.txt") into step4_allele
tuple val("${prefix}"), val(chr_name), path("ID-${prefix}${chr_name}.step4-HRC.txt") into step4_id
script:
"""
perl "${hrc}" -b "${prefix}${chr_name}.step4.bim" -f "${prefix}${chr_name}.step4.frq" -r "${hrc_site}" -h
"""
}

The output files are:
gsa-out/Chromosome-chr22.step4-HRC.txt     gsa-out/Force-Allele1-chr3.step4-HRC.txt  gsa-out/Position-chr22.step4-HRC.txt
gsa-out/Chromosome-chr3.step4-HRC.txt      gsa-out/ID-chr22.step4-HRC.txt            gsa-out/Position-chr3.step4-HRC.txt
gsa-out/Exclude-chr22.step4-HRC.txt        gsa-out/ID-chr3.step4-HRC.txt             gsa-out/Strand-Flip-chr22.step4-HRC.txt
gsa-out/Exclude-chr3.step4-HRC.txt         gsa-out/LOG-chr22.step4-HRC.txt           gsa-out/Strand-Flip-chr3.step4-HRC.txt
gsa-out/Force-Allele1-chr22.step4-HRC.txt  gsa-out/LOG-chr3.step4-HRC.txt

I am not sure why I am getting only 2, instead of 22 output files. With similar output channel previously generated 22 output files in the same script. For eaxmple:
process step1b {
publishDir params.publish, mode: 'copy'
input:
tuple val(prefix), val(chr_name), path(bfiles) from step1a_out
output:
tuple val("${prefix}"), val(chr_name), path("${prefix}${chr_name}.step1c.{bed,bim,fam}") into step1b_results
path("${prefix}${chr_name}.step1c.log") into step1b_logs
path("${prefix}${chr_name}.step1b.snplist.txt") into step1b_snplist
script:

"""
awk '{ if ((\\\$5=="T" && \\\$6=="A")||(\\\$5=="A" && \\\$6=="T")||(\\\$5=="C" && \\\$6=="G")||(\\\$5=="G" && \\\$6=="C")) print \\\$2, "ambig" ; else print \\\$2 ;}' "${bfiles[1]
}" | grep ambig > "${prefix}${chr_name}.step1b.snplist.txt"
plink1.9 --bfile "${prefix}${chr_name}.step1b" --exclude "${prefix}${chr_name}.step1b.snplist.txt" --make-bed --out "${prefix}${chr_name}.step1c"
"""
}

I suspect your HRC_1000G_check_bim and hrc_sites_file channels are queue channels when they should be value channels. Please again see how multiple channels work. Also, third-party scripts, like your Perl script, would be better placed in a directory called bin/ at the root of your project repository. Nextflow adds this directory to your $PATH, so scripts in this directory will be automatically accessible in your pipeline without needing to localize them using a separate channel. This way also avoids having to specify an absolute path to invoke them. Just make sure your script is executable: chmod +x your_script.pl
hrc_sites = file( params.hrc_sites_file )

process step5 {

    input:
    tuple val(prefix), val(chr_name), path(bfiles) from step4_frq
    path hrc_sites

    """
    your_script.pl \\
        -b "${prefix}${chr_name}.step4.bim" \\
        -f "${prefix}${chr_name}.step4.frq" \\
        -r "${hrc_sites}" \\
        -h
    """
}

Thank you @Steve. I just avoided using channel for those two files. I need to use them from absolute path, because other users also use this pipeline. However, it works now. One thing, I am using slurm executor with a config file, so far no problem with this. But I want to run the script in background (similar to sbatch submit, after submitting the job I will logout and the job will continue to run). How could I do this?

One thing, After finishing the jobs, transferring files to publish directory taking too much time. There are only 264 files but it is taking hours to transfer. Any solution for this?

No worries at all @zillurrahman. It's generally better to avoid specifying absolute paths in your process scripts. Avoiding absolute paths will mean that your process tasks can be run in isolation from each other and that your users could run your workflow on platforms like AWS Batch where there is no shared filesystem.

I don't have experience with Slurm, but you should be able to submit a job non-interactively using sbatch like you said. If it's submitted non-interactively, you can logout and the job running your workflow will continue on.

I'm not too sure about your last question; 264 files doesn't sound like a lot, but how much data are we talking? Could you be maxing out your network bandwidth? Are you writing to archival storage? Your sysadmin would be the right person to answer these sorts of questions and they might have some useful suggestions on how best run the workflow. For example, it might be better to write to a high performance disk array and then have a separate, but dependant, job copy these files to their final location.

You may also want to consider other publishing modes and parameterizing these, for example: mode: params.publish_mode. Consider also parameterizing the 'enabled' option so to let your users decide what outputs they want: enabled: params.publish_everything || params.publish_step5. HTH.


text:
What kind of BLAST do I need to do to accomplish this task?
(I am a complete newbie to bioinformatics, so please bear with me!).
I recently used BLAST to compare a 500 bp nucleotide sequence in S. cerevisiae to a bunch of other 500 bp nucleotide sequences in various fungi (these other sequences are homologs of the S. cerevisiae gene). Ideally, when I click "BLAST", I want to get nothing back. What I'm trying to double check is the assumption that my S. cerevisiae sequence is NOT similar to my other input sequences.
Here's what I ran on the computer, using the BLAST website:

After I run this query, I only get one sequence back:

which tells me that the sequence I got as an output needs to be excluded from my dataset, as it's too similar to the S. cerevisiae sequence.
The problem is, I need to run this type of query and analysis about 6000 more times for 6000 different S. cerevisiae genes and their homologs. I don't know how to do this programmatically. The Biopython packages seem to only compare one sequence against the whole NCBI database, and I was really struggling to understand how to do this via the command line (was struggling with the package installations a bit, and it also seemed to only be for one sequence input).
Is what I'm doing here a multiple sequence alignment? When I looked into programmatically doing an MSA, the output seemed to only be the alignment, but I want to see an output like the website gives me. I'm not really sure what I should be looking into.
Any guidance would help!

Double checking first, you are using DNA sequences: this is for a DNA related experiment as opposed to gene homology, right?
You most like would run a local blast (standalone) for this, not the web.
Majorly, albeit less sophisticated, Bio.pairwise2 (from Biopython) might actually be a simpler solution for you given your skillset as it does a local alignment and gives a score without the headache.
MSA, not what you need.

@MatteoFerla Sorry, I didn't see these messages for some reason! Yes, I'm using DNA sequences, and this is for a DNA experiment.

Clarification
Just making sure - you have ~6000 datasets - one for each S.cerevisiae gene of interest - which are made up of homologues of each respective gene? And you want to filter out genes that are too similar?
In that case, it would be helpful to know how you're defining a sequence as 'similar enough' for it to be filtered. Given they're all homologues, it's likely there's going to be a high baseline of sequence similarity. In your example, you consider 97.40% identity across the whole query as different enough to be filtered out, so I'm getting the feeling you specifically want to filter out identical sequences?
Also, when you talk about genes here, do you want to filter out gene products (i.e. proteins) that are too similar to the actual coding sequences (CDS)? The CDS  of two homologues might be X% different, but that difference might not be as high once translated into proteins due to how codons work. Proteins can go down to levels of similarity as low as ~30% identity and still have the same function, so depending on the context of your work (taxonomic, structural, functional) you might want to translate these CDS into protein, so further clarification would be helpful!
The actual bioinformatics question
To answer your question of programming, if the context of this comparison is always X S.cerevisiae gene vs. Y dataset of homologues, a MSA is definitely not the answer (you would do that if you wanted to see how all of the homologues compare to each other, rather than to specifically X gene). Also, using the web-client for NCBI BLAST will take a lot of time, so that's also not recommended for the kind of throughput you need.
Pairwise sequence alignments that you do on your local machine/server between X genes and all of its respective Y dataset sequences is probably the best way to go. I would suggest looking into needleall - it allows you to compare X sequences to Y sequences in a pairwise way, though it performs a global alignment rather than local like BLAST, so it might be slower.
I don't know what your computational resources are, so in case that might be too computationally heavy, the faster alternative would be to install the BLAST+ suite, which allows you to run BLAST queries like you did on the web browser, except on your own machine. Specifically, you could make your own custom individual BLAST databases out of each Y dataset using makeblastdb from the suite. Then you use blastn (or blastp if you decide to translate your sequences) on each dataset with its respective S.cerevisae sequence.
BLAST+ also allows you to customise the actual output as well, making it easier to parse your results programatically for when you want to remove sequences that are too similar. There are many examples online of how to do this, but have a go at it on your own  first with the BLAST+ documentation I linked :)


text:
How to use SBOL (Synthetic Biology Open Language)?
One of my students is working on the development of a SynBio tool to design and represent genetic circuits in a standardized format. I know there are some libraries out there.
What I am interested in finding out is to have a hands-on tutorial how to use an SBOL library (say, pySBOL) to represent genetic circuits in SBOL data and visual representations.

For pySBOL there is a tutorial you can follow in the documentation. Additionally, there is a youtube video walking through a similar tutorial, note the tutorial starts at about 1:09:00.
There are lots of additional tools and libraries out there such as SBOL Canvas (for visualisation), libsbolj (a java library for SBOL), and other libraries and tools listed on the SBOL website.

Thanks @PlantPromoter.


text:
How low does an E-Value need to be for a gene to code for the same protein as query?
I'm looking at a bacterial genome and trying to find specific PEGs (Protein Encoding Genes) by BLAST searching. I'm taking the genes found with the lowest E-value in order to psiBLAST them to check it codes for the same protein. I'm not exactly sure of how high the E-Value needs to be for me to conclude that the protein is not coded for by the subject gene.
Does anyone have any advice on how to make this process a little more efficient/scientific?
Thank you.
(By the way, I'm using BLAST with a BLOUSUM62 matrix.)

Just for the sake of clarity, could you please give more detail about what exactly your inputs are for BLAST?
For the first BLAST, I understand it as you having a list of CDS sequences that you're using as individual BLAST queries to a specific bacterial genome, keeping the lowest e-val hit for each query CDS, creating a list of potential PEGs.
For the PSI-BLAST, I understand it as you using the potential PEGs from the first BLAST as query to a list of proteins that you want to match with PEGs.
Let me know if I'm wrong!

I assume that you mean "is homologous to" by "encodes the same protein as". Identical ("same") proteins should be fairly trivial to identify using a tool such as tblastx.
You should probably not use E-value for this. E-value is similar to a p-value, in that it represents some aspect of the data (in this case a sequence and a database) compared to some null hypothesis (meaning in this case statistics on the database). Therefore, it depends both on the specifics of the data and the hypothesis, and there is no universal threshold that will work across all sequences and all databases. (See discussion of the meaning of e-value.)
E-value does not in any way represent the content of the sequence or its relationship to a sequence in the database. You could have E-values close to zero or E-values > 10 for a perfect sequence match or a very loose sequence match, depending on the input sequence and the database.
However, BLAST output does include other information that might be helpful. You could (for example) look at the length of the alignment and the number of mismatches/identities in the alignment. This will give you a pretty good guess at how similar the two sequences are; unless they are very divergent you would expect them to align across most of the length with somewhat high protein identity or similarity (>25% probably?). This process will still not work 100% of the time, but I think it is closer to your intent than any given approach involving E-values.


text:
Nextflow: run sbatch (slurm job script) inside a process
I was trying to run a slurm job script inside nextflow process. Is it possible? For instance, I have a job script 'PRE_PHASING_QC.JOB' and outside nextflow it runs perfectly, I want to integrate it in a nextflow script. May be my example code could be an explanation of what I want.
prephase='PRE_PHASING_QC.JOB'
process prephasing {
input:
file(prephasing) from prephase
"""
sbatch "${prephasing}
"""
}

The error:
sbatch: error: This does not look like a batch script
...

I tried to make my job script executable, put it in bin and then include it in the nextflow process according to @Steve. So, the script runs without any error but got no output. The original job script takes 3 files (bed,bim,fam) as input and produce multiple outputs. I guess that I have to declare outputs inside nextflow process. So far I tried as follows:
process prephasing {
publishDir params.out, mode: 'copy'
output:
tuple val(chr),val(num),val(steps),path("${chr}${num}.${steps}.{bed,bim,fam,log}") into out1
script:
"""
PRE_PHASING_QC.job
"""

It gives:
Error executing process > 'prephasing'
Caused by:
Missing value declared as output parameter: chr

My expected output would be chr{1-22}.step{1-10}.{bed,bim,fam,log}
#!/bin/bash
#SBATCH --job-name=prephase
#SBATCH --account=p697_norment
#SBATCH --time=1:00:00
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=7600M
#SBATCH --cpus-per-task=7
#SBATCH --array=1-22

source /cluster/bin/jobsetup

set -o errexit

module load plink

export prefix=chr${SLURM_ARRAY_TASK_ID}
export hrc_sites_file=/cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/HRC.r1-1.GRCh37.wgs.mac5.sites.tab.gz
export HRC_1000G_check_bim=/cluster/projects/p697/github/norment/moba_qc_imputation/software/HRC-1000G-check-bim.pl

# Step 1a (exclude duplicated), 1b (exclude indels and non-ACTG alleles), 1c (exclude ambiguous SNPs)

plink --bfile hce-ec-eur-fin-pass-qc --chr ${SLURM_ARRAY_TASK_ID} --out $prefix --make-bed || true
plink --bfile $prefix --list-duplicate-vars ids-only suppress first --out $prefix.stepla || true
plink --bfile $prefix --exclude $prefix.stepla.dupvar --make-bed --out $prefix.stepla || true
plink --bfile $prefix.stepla --snps-only just-acgt --make-bed --out $prefix.steplb || true

awk '{ if (($5=="T" && $6=="A")||($5=="A" && $6=="T")ii($5=="C" && $6=="G")||($5=="G" && $6=="C")) print $2, "ambig" ; else print $2 ;}' $prefix. steplb.bim | grep ambig > $prefix.steplb.snplist.txt || true
plink --bfile $prefix.steplb --exclude $prefix.steplb.snplist.txt --make-bed --out $prefix.steplc || true

# Step 2. Zero out Mendelian errors
plink --bfile $prefix.steplc --me 1 1 --set-me-missing --make-bed --out $prefix.step2 || true

# Step 3. Missingness and HWE + Step 4. Compute allele frequencies
# NB! --mind 0.02 filter is removed from the next command, to avoid sample QC per chunk, which may leave inconsistent number of individuals avail able per chunk.
plink --bfile $prefix.step2 --hwe le-6 --geno 0.02 --make-bed --out $prefix.step3 || true
plink --bfile $prefix.step3 --freq --make-bed --out $prefix.step4 || true

# Step 4-10. Run Oxford perl script and apply first 5 lines of its Run-plink.sh command, followed by '--update-name' command (update marker name s)
perl ${HRC-IRC_1000G_check_bim} -b $prefix.step4.bim -f $prefix.step4.frq -r ${hrc_sites_file} -h || true
plink --bfile $prefix.step4 --exclude Exclude-$prefix.step4-HRC.txt --make-bed --out $prefix.step5 || true
plink --bfile $prefix.step5 --update-chr Chromosome-$prefix.step4-HRC.txt --make-bed --out $prefix.step6 || true
plink --bfile $prefix.step6 --update-map Position-$prefix.step4-HRC.txt --make-bed --out $prefix.step7 || true
plink --bfile $prefix.step7 --flip Strand-Flip-$prefix.step4-HRC.txt --make-bed --out $prefix.step8 || true
plink --bfile $prefix.step8 --a2-allele Force-Allele1-$prefix.step4-HRC.txt --make-bed --out $prefix.step9 || true
plink --bfile $prefix.step9 --update-name ID-$prefix.step4-HRC.txt --real-ref-alleles --make-bed --out $prefix.step10 || true

When you say your script ran without error but you got no output files, are you looking in the workDir or publishDir? When I'm developing, I usually only write the output block declaration after I've tested the tool/script a couple of times. This is because your task process can still run successfully without an output block declaration. I can then look in the workdir to see what files have been created, etc, then finally declare them in the output block. Then I start on the next (downstream) process(es) and repeat.

The publishDir is empty and in the workDir there are only log files (chr.log, chr.step1.log, chr.step2.log etc.). It is strange!

So only files declared in your output block will get published into your publishDir. Ignore output in the publishDir while testing and developing. Or better, comment out the publishDir directive until you get your script block working in some capacity. In any case, your job script will show you what's going on. Happy to take a look if you need help.

Thank you @Steve. I added the job file as a png (not sure how to add text file). In line 21 hce-ec-eur-fin-pass-qc is the only input which create chr(1-22).(bed,bim,fam) files (for 22 chromosomes) which are used in the subsequent commands. I tried to replace hce-ec-eur-fin-pass-qc with $1 but got no output. Maybe you could give a better idea.

Please correct me if I'm wrong, but it looks like this script kicks off a job array. This isn't a good fit inclusion "as-is" into a Nextflow workflow, since it won't be portable at all. I think breaking it up into two scripts is the best option here. I'll see if I can give an example of this. Also, you can paste text in a code-block. I've updated your question with the text. Click edit to see the changes I've made.

See also: https://bioinformatics.stackexchange.com/editing-help

Just realized my edit hasn't been approved yet, but you should be able to see it.

Thanks @Steve . Yes, it starts a job array 1-22 for each chromosome. The initial input is one single bfile (bed,bim,fam), in step1a it creates 22 bfiles (bed,bim,fam) for each chromosome. The rest of steps operate on those 22 files to create final output. Any relatable example to break the job into two would be great.

No worries at all @zillurrahman. I've updated my answer with an example on how to do this. I also refactored some of your code, hopefully it still works...

Using Nextflow to submit jobs this way makes little sense, since:

sbatch exits immediately after the script is successfully transferred
to the Slurm controller and assigned a Slurm job ID. The batch script
is not necessarily granted resources immediately, it may sit in the
queue of pending jobs for some time before its required resources
become available.

Assuming your job script is basically a shell script, make it executable (e.g. chmod +x PRE_PHASING_QC.JOB) and move it into a folder called bin in the root directory of your project repository. Nextflow automatically adds this folder to your $PATH, making scripts in this folder accessible in your pipeline without needing localize them and without needing to specify an absolute path to invoke them. Your code then becomes:
process prephasing {

    """
    PRE_PHASING_QC.JOB [options] <args>
    """
}

Finally, add process.executor = 'slurm' to your nextflow.config file. You may need to configure some other options, like resource requests, but that should get your script integrated into your pipeline without much refactoring, if any.

Contents of PRE_PHASING_QC.sh:
#!/bin/bash

set -euo pipefail

hrc_sites_file="${1}"
input_prefix="${2}"
chrom="chr${3}"

# Step 1a (exclude duplicated), 1b (exclude indels and non-ACTG alleles), 1c (exclude ambiguous SNPs)

plink --bfile "${input_prefix}" --chr "${3}" --out "${chrom}" --make-bed
plink --bfile "${chrom}" --list-duplicate-vars ids-only suppress first --out "${chrom}.step1a"
plink --bfile "${chrom}" --exclude "${chrom}.step1a.dupvar" --make-bed --out "${chrom}.step1a"
plink --bfile "${chrom}.step1a" --snps-only just-acgt --make-bed --out "${chrom}.step1b"

awkcmd='
($5 == "T" && $6 == "A") || ($5 == "A" && $6 == "T") ||
($5 == "C" && $6 == "G") || ($5 == "G" && $6 == "C")
{
    print $2, "ambig"
    next
}
{
    print $2
}'

awk "${awkcmd}" "${chrom}.step1b" | grep "ambig" > "${chrom}.step1b.snplist.txt"
plink --bfile "${chrom}.step1b" --exclude "${chrom}.step1b.snplist.txt" --make-bed --out "${chrom}.step1c"

# Step 2. Zero out Mendelian errors
plink --bfile "${chrom}.step1c" --me 1 1 --set-me-missing --make-bed --out "${chrom}.step2"

# Step 3. Missingness and HWE + Step 4. Compute allele frequencies
# NB! --mind 0.02 filter is removed from the next command, to avoid sample QC per chunk, which may leave inconsistent number of individuals avail able per chunk.
plink --bfile "${chrom}.step2" --hwe 1e-6 --geno 0.02 --make-bed --out "${chrom}.step3"
plink --bfile "${chrom}.step3" --freq --make-bed --out "${chrom}.step4"

# Step 4-10. Run Oxford perl script and apply first 5 lines of its Run-plink.sh command, followed by '--update-name' command (update marker name s)
HRC-1000G-check-bim.pl -b "${chrom}.step4.bim" -f "${chrom}.step4.frq" -r "${hrc_sites_file}" -h
plink --bfile "${chrom}.step4" --exclude "Exclude-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step5"
plink --bfile "${chrom}.step5" --update-chr "Chromosome-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step6"
plink --bfile "${chrom}.step6" --update-map "Position-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step7"
plink --bfile "${chrom}.step7" --flip "Strand-Flip-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step8"
plink --bfile "${chrom}.step8" --a2-allele "Force-Allele1-${chrom}.step4-HRC.txt" --make-bed --out "${chrom}.step9"
plink --bfile "${chrom}.step9" --update-name "ID-${chrom}.step4-HRC.txt" --real-ref-alleles --make-bed --out "${chrom}.step10"

Make PRE_PHASING_QC.sh executable (chmod +x PRE_PHASING_QC.sh), and replace your job script with the following. This file won't be needed by your Nextflow process, but shows you how the above shell script can be invoked:
#!/bin/bash
#SBATCH --job-name=prephase
#SBATCH --account=p697_norment
#SBATCH --time=1:00:00
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=7600M
#SBATCH --cpus-per-task=7
#SBATCH --array=1-22

set -e

. /cluster/bin/jobsetup

export PATH="/cluster/projects/p697/github/norment/moba_qc_imputation/software:${PATH}"

module load plink

cd /path/to/hce-ec-eur-fin-pass-qc

/path/to/PRE_PHASING_QC.sh \
    /cluster/projects/p697/projects/moba_qc_imputation/resources/HRC/HRC.r1-1.GRCh37.wgs.mac5.sites.tab.gz \
    hce-ec-eur-fin-pass-qc \
    "${SLURM_ARRAY_TASK_ID}"

Copy PRE_PHASING_QC.sh and the Oxford Perl script, HRC-1000G-check-bim.pl, into your project's bin/ directory. You may also need to make the Perl script executable. Then your Nextflow process becomes:
hrc_sites_file = file(params.hrc_sites_file)

Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size:3 )
    .set { plink_data }

Channel
    .of(1..22)
    .set { chroms }

process prephasing {

    input:
    tuple val(prefix), path(bfiles) from plink_data
    each chrom from chroms
    path hrc_sites_file

    """
    PRE_PHASING_QC.sh "${hrc_sites_file}" "${prefix}" "${chrom}"
    """
}

Once the above runs successfully, you can declare your process outputs. Please note that none of the above was tested; hopefully not too many typos, but YMMV.

Thanks @Steve . I think you understand my explanation and your answer makes sense. I have updated the question with new code, error and expected output. Hope it helps to debug

@zillurrahman: You get the error Missing value declared as output parameter: chr because you've defined 'chr' in your output block, but this value wasn't defined in your input block (or globally in your script block). You didn't show what your input block looks like. I suspect it should look something like: tuple val(chr), val(num), val(steps), path(bfiles). I don't know what you job script does, but if your job script references absolute paths to files, you'll need to refactor it so that they can be passed in as variables on the command line.

My job script takes bfiles as input (e.g plink --bfiles input_bfiles --other argument --output). I just tried to modify this as : plink --bfiles $1 and tried to use input channel from nextflow. """  PRE_PHASING_QC.job "${bfiles}" """" In this case nothing happens, no error or success.

From our previous discussions, I thought plink took a prefix for it's '--bfile' option? Note this isn't the same as '--bfiles'. If I'm correct, I think your first task will be to work out how to supply a prefix on the command line so that it gets passed down to your plink command. In your NF process, I think you'll want: PRE_PHASING_QC.job "${prefix}". And in your script I think you'll want: plink --bfile "$1"

Thanks @Steve . I tried nextflow run prephasing1.nf --hrc_sites_file "path/to/hrc_sites.tab.gz" "path/to/hce-ec-eur-fin-pass-qc" "1" It give me Missing 'bind' declaration in input parameter . Should we use tuple(val),val(input_prefix),path(bfiles) from some_channel ? How could I declare input channel in this particular case? It seems I could run chromosome by chromosome? That if I put "1" in nf run parameter, it will only run for chromosome 1? Should we declare it in in_prefix=file(params.input_prefix) and params.chr=1 then we could change/call the arguments in command line?

Thanks @Steve So far I declared input_prefix and chrom as hrc_sites_file was declared in your script and specify them in command line. I got output for chromosome 1, chr1.step1a, chr1.step1b but in the next step it shows that no such file as chr1.step1b but it exists in my work/ directory. Any thought? When I use chrom 2, I got chr2.step1a, chr2.step1b and same error in the next step.

No worries @zillurrahman. I just updated the NF code to process all chromosomes like what we have done previously. I think you'd rather avoid running each chromosome one at a time. In your first comment above, you get the Missing 'bind' declaration in input parameter error since you're missing parameters for the "path/to/hce-ec-eur-fin-pass-qc" and "1" options.

Ohh! There was some typo in the PRE_PHASING_QC.sh script. I ran without the slurm, identfied those typos, edited them and now everything is Ok. Now, How could I declare the output files so that I get all the outputs in a single directory? So far, I tried tuple val("${chrom}"), val(steps), path("${chrom}.${steps}.{bed,bim,fam}") into qc1_bfile and tuple val(prefix), val("${chrom}"), path("${prefix}-${chrom}.step4-HRC.txt") into qc1_hrc But got error:  Missing value decalred as output parameter: steps How could I decalre steps so that all outputs from each step gathered in a folder?

Got it! I did as following: tuple val("${chrom}"), path("chr${chrom}*.{bed,bim,fam}") into qc1_chr and it works. Only thing is, file transfer taking too long. After finishing the process Waiting files transfer to complete (2 files) it is taking too much time. If I check output folder, I can see all my desired outputs but the waiting files transfer to complete message is still on the screen and I can not do anything in this tab. Any suggestion? However, Thanks for helping all the way!

Well done getting it working @zillurrahman. Is the issue that the file copying is taking too long, or is the issue just that while you're waiting you cannot use the current terminal? If it's the latter, it's not really a problem since the job can be run non-interactively (e.g. submitted as a job to your Slurm cluster), in the background (e.g. using screen or nohup) or can be backgrounded after i'ts started (e.g. using , then typing 'bg') to background the process.

Thanks @Steve. I have managed to run it using sbatch. It is ok now.

Hi, This process creates multiple output files (chr1-22.step1-10.bed|bim|fam) which are saved into qc1_chr  output channel. Now, say I need only those files which has string step10. How could I get that subset into a new channel? I need it into a new channel because I need to merge them with elements of other channels. Did my explanation make sense?

@zillurrahman you could use separate output channels for each step, but you might not need the other step outputs. Another way would be to select just the files you need from the qc1_chr output channel using the filter operator. There's lots of ways to do this - for example: qc1_chr.filter { bfile -> ["bed", "bim", "fam"].any { bfile.name.endsWith("step10.${it}") } }.view(). HTH.

Thanks @Steve. It gives me No such variable: Exception evaluating property 'name' for java.util.ArrayList, Reason: groovy.lang.MissingPropertyException . If I remove name from the above snippet, I got Unknown method invocation 'endsWith' on ArrayList type. Any idea?

Sorry @zillurrahman I didn't see your qc1_chr declaration above. I think you want: qc1_chr.filter { name, bfiles -> bfiles.collect { it.getBaseName() }.every { it.endsWith("step10") } } HTH.

Thanks @Steve. It seems ok now, getting no error. Another issue. I merged other channels with it (each channel for 22 files, 1 for each chromosome) shape_para=shape_in.merge(shape_map).merge(shape_ref).merge(shape_legend).merge(shape_sample) process phasing { input: each file(penta) from shape_para output: path("*step10*") into phased script: def (in, map, ref, legend, sample) = penta """ shapeit2 -T 32 --input-bed "${in}" --input-map "${map}" --input-ref "${ref}" "${legend}" "${sample}" --duohmm -O "${in}.phased" --seed 54321 """ } I am getting no output or error!

I'm not sure exactly @zillurrahman, sorry. Your input declaration doesn't look right, but I'd still expect some output or log on stderr if your process is executed. If you can ask a new question, I'll see if I can add an answer. Please also show the contents of shape_in, shape_ref, shape_legend and shape_sample. Cheers.

Thanks @Steve I posted a new question with code for each channel here; https://bioinformatics.stackexchange.com/questions/15680/merge-multiple-channels-each-emits-multiple-files-in-nextflow


text:
Combine two alleles to genotype locus?
I have a following fasta (produced by stacks populations):
>CLocus_10056_Sample_934_Locus_12529_Allele_0 [sample_934; groupI, 49712]
TGCAGGCCCCAGGCCACGCCGTCTGCGGCAGCGCTGGAAGGAGGCGGTGGAGGAGGCGGCCAACGGCTCCCTGCCCCAGAAGGCCGAGTTCACCG
>CLocus_10056_Sample_934_Locus_12529_Allele_1 [sample_934; groupI, 49712]
TGCAGGCCCCAGGCCACGCCGTCTGCGGCAGCGTTGGAAGGAGGCGGTGGAGGAGGCGGCCAACGGCTCCCTGCCCCAGAAGGCCGAGTTCACCG

as both alleles are represented separetely, I want to combine both genes using IUPAC nucleotide code.
Is there a tool for it or is it more useful to do it via a custom script?


text:
Do proteins move inside the cell?
Do proteins move inside the cell and how fast do they move or move? Please support your answer with articles or books.

I’m voting to close this question because while it may be vaguely related to biology, it isn't related to bioinformatics.

Proteins move in cells or get transported in the cell. You can go to NCBI (https://www.ncbi.nlm.nih.gov/), and start searching for protein dynamics within cells. Also there is a technique called FRAP (https://en.wikipedia.org/wiki/Fluorescence_recovery_after_photobleaching) which is used to measure kinetics and dynamics.
Also, this differs probably a lot between the cellular system you are asking for. You can imagine that some cells have a more active transport system because stretched out (neurons) vs others with maybe less active transport systems, because of a more static nature.


text:
Toptable error, wont recognize condition
I am getting a rather strange error from topTable.    When I run my code I get Error in fit$coefficients[, coef] : subscript out of bounds from topTable as if it is not recognizing the condition I am asking it to return the top ten DEGs of. I have provided a bit more context and some descriptions of my variables below.  Let me know if you need any more information or if anything is unclear.
Variables:
dge is a DGEList object of RNA seq data.
dge$counts is the raw counts.
dge$samples$dge_meta is the metadata and log2 transformed counts.
> colnames(dge$samples$dge_meta)
 [1] "Lab.sample.ID"
 [2] "Accession.Number"
 [3] "Alias.used.in.some.figures"
 [4] "Developmental.stage"
 [5] "Infected"
 [6] "Batch"
 [7] "Trimmed"
 [8] "Number.of.reads.that.pass.Illumina.filter"
 [9] "Total.number.of.reads.mapped"
[10] "%.of.total.reads.mapped"
[11] "Reads.mapped.to.T..cruzi.Esmeraldo.haplotype"
[12] "Reads.mapped.to.hg19"
[13] "%.of.mapped.reads.belonging.to.T..cruzi.Esmeraldo.haplotype"
[14] "%.of.mapped.reads.belonging.to.hg19"
[15] "sample"
[16] "hpi"

Lastly keep_it is a binary vector used to filter our zero count genes.
Code:
mm <- model.matrix(~~hpi*Infected, data = dge$samples$dge_meta) #making our design matrix

more_than_1cpm_model <- voom(dge$counts[keep_it,], mm, plot = TRUE) #this one is filtered by the 25% greater than 1 cpm

cutoff <-  1e-05

empBayesfit <- eBayes(fit_gt1cpm)

topTable(empBayesfit,
         coef = c("Infected"),
         p.value = cutoff, number = 10)

Error message:
Error in fit$coefficients[, coef] : subscript out of bounds

To understand the context better I ran rlang::last_error and rlang::last_trace as well.
> rlang::last_error()
<error/dplyr_error>
arrange() failed at implicit mutate() step.
* Problem with `mutate()` input `..1`.
x Input `..1` must be a vector, not a `standardGeneric` object.
i Input `..1` is `samples`.
Backtrace:
Run `rlang::last_trace()` to see the full context.
> rlang::last_trace()
<error/dplyr_error>
arrange() failed at implicit mutate() step.
* Problem with `mutate()` input `..1`.
x Input `..1` must be a vector, not a `standardGeneric` object.
i Input `..1` is `samples`.
Backtrace:
     x
  1. +-dplyr::arrange(x, samples)
  2. +-dplyr:::arrange.data.frame(x, samples)
  3. | \-dplyr:::arrange_rows(.data, dots)
  4. |   +-base::withCallingHandlers(...)
  5. |   +-dplyr::transmute(new_data_frame(.data), !!!quosures)
  6. |   \-dplyr:::transmute.data.frame(new_data_frame(.data), !!!quosures)
  7. |     +-dplyr::mutate(.data, ..., .keep = "none")
  8. |     \-dplyr:::mutate.data.frame(.data, ..., .keep = "none")
  9. |       \-dplyr:::mutate_cols(.data, ...)
 10. |         +-base::withCallingHandlers(...)
 11. |         \-mask$eval_all_mutate(quo)
 12. +-dplyr:::abort_glue(...)
 13. | +-rlang::exec(abort, class = class, !!!data)
 14. | \-(function (message = NULL, class = NULL, ..., trace = NULL, parent = NULL, ...
 15. |   \-rlang:::signal_abort(cnd)
 16. |     \-base::signalCondition(cnd)
 17. +-(function (e) ...
 18. | \-rlang::abort(...)
 19. |   \-rlang:::signal_abort(cnd)
 20. |     \-base::signalCondition(cnd)
 21. \-(function (cnd) ...

The error is telling you that there isn't any Infected column on the model.matrix.
Check the column names of mm. Also not sure that having two ~ will work well.

I did notice the ~~ mistake thanks for pointing it out.


text:
Why do molecular generation models maximize “penalized logP” as a measure of drug-likeliness?
I found that Lipinski's rule of five states that Log P (octanol-water partition coefficient, lipophilicity measure) usually should not exceed 5.
Many papers about drug discovery machine learning models tell about maximization of "penalized logP", but the paper they refer in the end does not contain any information on it. Let me show that.
Let's check 1, here we read about penalized log P:

The second task is to produce novel molecules with desired properties. Following (Kusner et al., 2017), our target chemical
property y(·)is octanol-water partition coefficients (logP) penalized
by the synthetic accessibility (SA)score and number of long cycles
y(m) = logP(m)−SA(m)−cycle(m)

where cycle(m)counts the number of rings that have more than six atoms.

Then we go into Kusner et al., 2017:

For the second optimization problem, we follow (Gómez-Bombarelli et al., 2016b) and optimize the drug properties of molecules. Our goal
is to maximize the water-octanol partition coefficient (logP), an
important metric in drug design that characterizes the drug-likeness
of a molecule.As in Gómez-Bombarelli et al. (2016b) we consider a
penalized logP score that takes into account other molecular
properties such as ring size and synthetic accessibility (Ertl &
Schuffenhauer, 2009).

Here is Gómez-Bombarelli, 2016 and it seems to contain nothing on the topic (I can be mistaken, please point me then), except

The objective we chose to optimize was 5×QED−SAS, where QED is the Quantitative Estimation of Drug-likeness (QED),37 and SAS is the
Synthetic Accessibility score.36Thisobjective represents a rough
estimate of finding the most drug-like molecule that is also easy to
synthesize.

Other examples of the papers also leading to the same papers are 2, 3, etc.
What does penalized logP mean for drug-likeliness and why it should be maximized?

A positive logP is good (hydrophobic), but high results in ADMET problems (cytochromes break overly hydrophobic compounds) and aggregation (gunk, tar, crud are words happily thrown around). However, in most cases the logP is not greater than +3 because good binding is result of hydrogen bonds etc —some bad algorithms will suggest a polyaromatic hydrocarbon that makes lots of pi-pi interactions but these are few.
Confusingly, a high synthetic accessibility score is bad (hard). 3 is the average and 10 is basically something impossible. This score is how hard it is to make (i.e. pricier), which is one of the major drivers.
The number of cycles is number of rings. Even if you get pi-pi, cation-pi and sulfur-pi bonds you also get aggregation due to stacking (a crud). So is somewhat bad.
So basically that score is
high=good minus high=bad minus high=bad

If difference in relative Gibbs energy of binding is thrown it, negative ∆∆G is good, so that would be subtracted. As its atom count dependent ligand efficiency is better, but that ranges from 0 to -1.2 kcal/mol. So would need scaling.

You were trapped by arXiv versioning convention that shows you the latest version of the preprint as default. It seems that this experiment was removed from the final publication.
If you have a look at the v1 (https://arxiv.org/abs/1610.02415v1) instead of the latest version (v3), you find the explanation:

As a simple example, we first attempt to maximize the water-octanol
partition coefficient $logP$, as estimated by RDkit. [43] $logP$ is an
important element in characterizing the druglikeness of a molecule,
and is of interest in drug design.  To ensure that the resulting
molecules to be easy to synthesize in practice, we also incorporate
the synthetic accessibility [44] (SA) score into our objective. Our
initial experiments, optimizing only the $logP$ and SA scores, produced
novel molecules, but ones having unrealistically large rings of carbon
atoms. To avoid this problem, we added a penalty for having carbon
rings of size larger than 6 to our objective. Thus our preferred
objective is, for a given molecule m, given by:
$J(m) = logP(m) − SA(m) − ring-penalty(m)$,          (1)
where the scores $logP(m)$, $SA(m)$, and $ring-penalty(m)$ are normalized to
have zero mean.


text:
How to share a collection of sbol designs in SynBioHub
We currently work on a project that will generate a collection of recombinant plasmids targeting over 3500 genes of the cyanobacterium Synechocystis sp. PCC 6803. We would like to share these designs in SynBioHub but we are not sure how best to go about it. The plasmids have been designed and now they are starting to be produced in batches, checked and finally transformed into Synechocystis to induce a single gene deletion by homologous recombination.
During the project, we will have:

Stage 1 – the initial plasmids designs, waiting for production. They
are worth sharing as the community could influence the production
queue.
Stage 2 – the assembled plasmids. In some cases they will have
different structure than initial designs from stage 1. If nothing
else their description should reflect the status change to
'produced'.
Stage 3 – plasmids are verified by sequencing (also to
check the barcodes are correct) and sequences attached.
Stage 4 – plasmids will be transformed into Synechocystis and checked by PCR to
see if i) the target gene has been knocked out and ii) if the strain
is segregated (i.e. all copies of the genome in the strain contain
the gene deletion).  This information should also be recorded with
the description.

Our initial plan was:

Create a collection with initial designs (component definitions) Stage 1
As stage 2 progresses add new versions of udpated designs
Stage 3, 4 update plasmids descriptions (synbio desc or notes) and attach additional files (e.g. from sequencing).

However, our current understanding is the public collections are frozen so stage 2,3,4 cannot be accommodated.
We are seeking advice how to plan for sharing this resource and some clarifications.

Is there a way of making collection readable by others without freezing its content?
What operations are possible on public component definitions, can the free text descriptions/notes be updated or file attached?
Can new component definitions or their versions be added to a published collection
For 2) 3) who can perform these operations, the collection’s creator or SynBio curators?

Looking forward to input from SynBioHub community

Welcome. Typically we ask that questions only contain a single question. Can you please consider condensing your post down to a single question and/or asking the remaining question in a separate post. Thanks.

Knowing something about SynBioHub, I think this is well-answerable as one question, because it's all focused on its publication workflow. I'm not sure of the answers myself, but I'm going to try to get one of the folks who does know to answer.

There is a balancing act between allowing all edits, which defeats the purpose of having URIs point to static objects, and allowing no edits, which defeats the purpose of providing some updated information about an object.  SynBioHub resolves this issue by allowing a limited amount of edits after publication.  To answer your questions specifically:

Yes, you can share objects by generating a share link OR adding another owner to an object.  Sharing can be done by the original owner of the object.
There are some mutable fields that any owner of the object can change, these include the description, notes, source, and references in the Details section.  You can also add attachments.
Yes, when you publish (Make Public) a new object or collection, you have the option of making it public into an existing collection OR a new collection.  A person with curator privileges on the SynBioHub instance you are working on can perform this operation.
Answered above.


text:
Removing duplicate FASTA sequences based on headers with Bash
I used the following command to remove duplicate FASTA sequences based on the header sequence:
paste -d $'\t' - - <fastaFileWithNoLinebreaksInSeq | sort -t $'\t' -uk1,1 | awk 'BEGIN{FS="\t";OFS="\n"}{print $1,$2}'

After this command, the output result is shown below:
>_10|_10
AAAGCAAAACACAGAATTGACGGAAAGACTTACGTTATTAAACGTGTTAAATATAATAACGAGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTC
>_10|_11
ATTATGATCCTGAGACCAGTGATGATTCTCTTGAGAGCAGTGATTATGATCCTGAGAACAGCAAAAATAGTTCAAGGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACACA
>_10|_12
TTTGGAACTCTTTGAACAAATAACAAAAGGGGTGGATTATATACATTCAAAAAAATTAATTCATAGAGATCTTAAGGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACACA
>_10|_13
ACTTGTAACATCTCTGAAAAATGATGGAAAGCGAACAAGGAGTAAGGGAACTTTGCGATACATGAGCCCAGAACAGGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACACA
>_10|_14
AACGAATTTCTTCGCAAGACTATGGAAAGGAAGTGGACCTCTACGCTTTGGGGCTAATTCTTGCTGAACTTCTTCATGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACAC

This is my actual FASTA file: I don't want to lose header information from my FASTA sequence after removing duplicate sequences based on header information.
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGCATGAAC
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGAAGTTG
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGGTGGGG
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGTAATTA
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGTCATGG
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGTTCCTTG
>ENST00000642616_8|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTATTTGGTGCCTATCCTAAAACTTCCAGTGGGTTTTCTTGCCTTCCAGTGGCCATATAAAGGGCGTCCTGAAAAC
>ENST00000642616_7|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTCATGAACAAGGCTCAGTTATGATGTTTGAGCACCTTATTAATTCTAATGGAATTAAGCCTGTCATGGAACAATA
>ENST00000642616_1|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGCGCCTGCGCGCGGGTCCGGCGCCGAGGTTCTTGACTGCTGTGCCGGACGCCAGGTGTAGCCATGCAGCGAGCC
>ENST00000642616_9|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGGACTGCCATCATCTTGGAATCCAAAATAATTTTGATTACAAGCGCTTTATTAAGTTTGCCCGTGTCTGTGAAG
>ENST00000642616_5|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGGTGGGGTATCTAGCAGGATGTCTAGTTCACGCACTGGGTGAAAAACAACCAGAGCTGCAGATAAGTGAACGAG
>ENST00000642616_4|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTAATTAATGATCCTATCCATGGCCACATTGAGCTCCACCCTCTCCTCGTCCGAATCATTGATACACCTCAATT
>ENST00000642616_6|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTCATGGGCCATTTTCTCACATGTTTGATGGACGATTTATTCCACTTGCTCGCCCGGAGGTGAAATGGACGATT
>ENST00000642616_6|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTCATGGGCCATTTTCTCACATGTTTGATGGACGATTTATTCCACTTGCTCGCCCGGAGGTGAAATGGACGCAT
>ENST00000642616_6|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTCATGGGCCATTTTCTCACATGTTTGATGGACGATTTATTCCACTTGCTCGCCCGGAGGTGAAATGGACGGAA
>ENST00000642616_3|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTTTCCTTGGGGGAGAGGAAGAAGCTGCTTAGTTATATCCAGCGATTGGTTCAAATCCACGTTGATACAATGAAGA
>ENST00000642616_3|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTTTCCTTGGGGGAGAGGAAGAAGCTGCTTAGTTATATCCAGCGATTGGTTCAAATCCACGTTGATACAATGAAGC
>ENST00000642616_3|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTTTCCTTGGGGGAGAGGAAGAAGCTGCTTAGTTATATCCAGCGATTGGTTCAAATCCACGTTGATACAATGAAGG
>ENST00000462861_2|ENST00000462861_3
AAAAAACTCTTCTACAGAAATTACTCTCAAAGAAACCTGAGGATCGACCTAATTTTTCACAGACCTACGGGATGGCATCATCTCAGATATATTTGATAAAAAAG

Can you please clarify in what form/language you need the solution? I see that the solutions provided so far are not in the same form/language as your example.

I don't understand what you actually have and what you need. What output are you expecting from this example? Did the command you used work or not? I tried running your command and it seems to work fine, I don't get the output you show in the first example. What operating system are you using?

If you want to avoid using extra libraries for any reason, you can just use a simple Python script (version 3.6 and above) to do this:
fr = open("dup_test.fasta", "r")
fw = open("dup_edited.fasta", "w")
seq_dict = {}
curr_header = ''

for line in fr:
   line = line.strip()
   if line[0] == '>':
        if line not in seq_dict:
          seq_dict[line] = ''
          curr_header = line
   else:
        if seq_dict[curr_header] == '':
          seq_dict[curr_header] = line
          fw.write(f"{curr_header}\n{line}\n")

fr.close()
fw.close()

You can use seqkit for this purpose.
seqkit rmdup -n seqs.fa -o seqs_without_duplicate.fa

Very cool kit, thanks for sharing!

@Bioinfomrat you're welcome :) Your Script is also a quite nice solution.

One way using awk:
awk '/^>/ { f = !a[$0]++ } f' seqs.fa

Explanation:
The above is basically the same as:
awk '/^>/ { f = !($0 in a); a[$0]++ } f' seqs.fa

That is, for each header line (lines that start with a '>' character), set a flag, 'f' to true if the current line is not in an array, 'a'. Then add the (header) line to the array. Note that in this form, it's not necessary to increment the value for each element - the following would be sufficient:
awk '/^>/ { f = !($0 in a); a[$0] } f' seqs.fa

The default action is to print, and we print lines only when 'f' is true. This is the same as:
awk '/^>/ { f = !($0 in a); a[$0] } f==1 { print $0 }' seqs.fa

Thanks for the explanation! This is a nifty solution. Any idea on its limitations? I can see that it will work as long as the array can be managed in memory, which is quite a small limitation.

No worries @RamRS. Nothing comes to mind. It should work on most platforms. All bets are off if the input is badly formatted like in the example. It's worth mentioning also that AWK arrays are associative, so only unique keys will be stored in memory. This might be an issue if the input contains a large number of sequences and the number of duplicates is low.

Also, note that only the header lines are kept in memory, not the sequence, meaning that it is vanishingly unlikely that you will ever come across a file with enough sequences not to fit in RAM on any modern machine.


text:
Understanding ViennaRNA RNAdistance scoring table
I'm trying to compare the output of 2 different algorithms of RNA structure prediction
(my implementation of Nussinov vs RNA-mfold algorithm) using the RNAdistance algorithm that is part of ViennaRNA package.
I'm getting a distance score of 38 between Nussinov prediction and the actual NMR structure, and a score of 14 between the m-fold algorithm and the NMR structure and I'm trying to understand whether it is a big difference.
In order to do that I have to look at the scoring table that is used by RNAdistance and found the following score table in the source code:
RIVATE CostMatrix  UsualCost =
{

/*    Null,       U,        P,        H,        B,        I,        M,        S,        E,        R     */

   {        0,        1,        2,        2,        2,        2,        2,        1,        1, DIST_INF},   /* Null replaced */
   {        1,        0,        1, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF},   /* U    replaced */
   {        2,        1,        0, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF},   /* P    replaced */
   {        2, DIST_INF, DIST_INF,        0,        2,        2,        2, DIST_INF, DIST_INF, DIST_INF},   /* H    replaced */
   {        2, DIST_INF, DIST_INF,        2,        0,        1,        2, DIST_INF, DIST_INF, DIST_INF},   /* B    replaced */
   {        2, DIST_INF, DIST_INF,        2,        1,        0,        2, DIST_INF, DIST_INF, DIST_INF},   /* I    replaced */
   {        2, DIST_INF, DIST_INF,        2,        2,        2,        0, DIST_INF, DIST_INF, DIST_INF},   /* M    replaced */
   {        1, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF,        0, DIST_INF, DIST_INF},   /* S    replaced */
   {        1, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF,        0, DIST_INF},   /* E    replaced */
   { DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF, DIST_INF,        0},   /* R    replaced */

};

I also found the enter link description herefollowing explanation, but I do not understand what is the meaning of the table letters: U,P,H,B,I,MS,E,R ?
What is the meaning of the score of [H,B] = 2 for example?
If anyone can explain the table it will be a great help in my project, thank you!

Ok, the answer was right here and I missed it, posting it in case anyone else will do the same mistake as me in the future.
U = unpaired
P = paired
H = hairpin
B = bulge
I = internal loop
M = multiloop
S = stem (or stack)
E = external loop


text:
File format of substitution matrix in clustalw
I need to set the substitution matrix used by command line CLUSTALW when comparing DNA sequences to:
 0 -1 -1 -1
-1  0 -1 -1
-1 -1  0 -1
-1 -1 -1  0

from my understanding I need to pass the location of a file with the matrix to the -dnamatrix parameter, however, I cannot find out what format the file is supposed to be in. Any idea?

This example seems to have information on this, looks like it should have this format:
   A   G   C   T   U   *
A  1   -1  -1  -1  -1  -1
G  -1  1   -1  -1  -1  -1
C  -1  -1  1   -1  -1  -1
T  -1  -1  -1  1   -1  -1
U  -1  -1  -1  -1  1   -1
*  -1  -1  -1  -1  -1  -1

However, the page I linked is specifically about a bug regarding ClustalW ignoring custom made matrices. It's a pretty old post so it might have been fixed since. However, if you encounter that bug, that same page thankfully does contain a workaround!


text:
Get Salmon mapping/alignment summary
With HISAT2, after the alignment of fastq files you get an alignment summary like this:
20000 reads; of these:
  20000 (100.00%) were unpaired; of these:
    1247 (6.24%) aligned 0 times
    18739 (93.69%) aligned exactly 1 time
    14 (0.07%) aligned >1 times
93.77% overall alignment rate

Is it possible to get this information with Salmon?

The salmon output log contains a line with [info] Mapping rate which allows to calculate how many reads went mapped or unmapped. Is that sufficient? It should be in the logs folder of the directory salmon creates per sample.


text:
fastq file format unknown
I have paired-end fastq files some of which seem to be in a weird format (from a collaborator, not a public database). When I cat the file I get what seem to me to be binary syntax but not sure?

However when I cat into the other pair-end the output seems fine

Any help/suggestion/solution would be greatly appreciated.
It does seem to be compressed with a missing suffix:
file pair1.fastq
   pair1.fastq: data.

When I run file pair1.fastq on the problem fastq file I get the following output: pair1.fastq: data. When I run it on the other pair-end file which doesn't have any issues I get pair2.fastq: gzip compressed data, extra field

gzip compressed data, extra field is probably produced using bgzip. If it were simple gzip, the extra field part wouldn't be there.

This indeed looks like some sort of binary file. I have seen bunch of those when people renamed .fastq.gz to .fastq without actually uncompressing them. However, as you say, this is not your case, as file pair1.fastq would identify the file as gzip compressed data.
Looking at possible return values of file comments here, data basically mean the command could not figure out anything about it beside being some sort of binary data file, but I suspect it probably somehow broken file. Honestly, at this point I think you should check the source of the data. If it's from a public database, re-download the data. If it's your own sequencing, contact the sequencing provider and check the data with them. I don't see any way to rescue this file (which does not necerasrily mean there is nothing you can do).

thanks again for the quick response. I was afraid that that might be the case. unfortunately the file came from a collaborator not a public database. I will have to bring that to their attention.


text:
How can I use statistics to compare microbial phenotypes?
Note: this question has also been asked on Biostars
I am currently trying to create a theoretical argument that a microbe's phenotype can affect gene expression in their host. I have 5 species of microbes, each with a different COG (Cluster of Orthologous Genes) profile.  I've been using python and a large database (DualSeqDB) of significant changes in gene expression to compare gene expression in a host and the COG profiles of my 5 species.
This is an example of the dataframe I'm working off of:

columns are counts of COGs present in a microbial species and rows are species
I was wondering if there is any way to compare the species and their COG profiles. My end goal is to draw parallels between how prevalent a certain phenotype is and how much that species can affect a host's gene expression.
If this seems infeasible, would a better approach be to BLAST the species and compare their genome sequences phylogenetically? (this would change the nature of my argument, so I'm hesitant)
I'm not the best with statistics so I'm sorry if there's a really obvious answer I'm not thinking of. If anyone knows of any resources I could use to read up on this stuff more on my own, I would really appreciate that as well.


text:
Calculating most abundant transcript from RNA-Seq data
vcf2maf uses VEP to annotate variants, and I believe selects the default Ensembl transcript to use for annotation. Sometimes the transcript that VEP selects is not the transcript I'm interested in, usually because the selected transcript is not the most highly expressed transcript in my tissue of interest (skin). vcf2maf allows you to provide a transcript override list so that VEP annotates the variant using the specified transcripts instead.
I have several skin samples sequenced with bulk RNA-Seq. I want to estimate the average abundance for each transcript across all samples and then use these abundances to rank transcripts from most to least abundant. Then I will use the most abundant transcript as the default VEP transcript. I plan to use salmon or kallisto to quantify transcript abundance. Should I use TPM or normalized counts to calculate average expression?
My initial thought is to use normalized counts (generated by DESeq2 from raw counts). Are there any problems with this approach? GTEx displays transcript abundance with average TPM, but I thought TPM was inappropriate to use across samples because it doesn't account for between sample differences.
Update: I forgot to mention I also tried using TPM ranks like @ATpoint describes. I haven't fully compared how this compares to transcripts identified by normalized counts, but the initial genes I checked showed good concordance between methods

Transcript abundance quantification is a tricky topic since a read often could belong to several transcripts, so any "count" is a best guess as to which transcript it actually originates from. That being said, there are tools that can help you here:

salmon (as you mentioned) to quanitfy. Run it with --numGibbsSamples 50 (or higher if your computer has the hardware. Also run with -d for the next tool
terminus is a tool that also originates from the COMBINE-lab, It's relatively experimental, but it applies well to your use case. It will take the salmon output and assign counts to transcript groups. If a set of transcripts can be clearly quantified above a certain confidence, that transcript will stay unique. If it is confounded with other transcripts, those will be placed in a transcript group and their counts will be summarised. More info here

You are right in that TPM, FPKM etc. are not appropriate measures to compare across samples. TMM, or DESeq2's median ratio method (what you get with counts(dds, normalized=TRUE) are a step in the right direction, but are still inappropriate across samples. I'd suggest this often posted blog post in the topic.
There are no perfect solutions here since RNA-Seq is inherently a relative quanitifcation method so this will always be an issue. IMO, use DESeq2::varianceStabilizingTransformation(), DESeq2::rlog() or limma::voom() after terminus and rank your transcripts that way.
Just keep the inherent transcript incertainties in mind.

I did not know about terminus, it looks really interesting. Could you explain why TMM/median ratio method are still inappropriate across samples? I thought they were designed for between sample normalization.

Don't get me wrong, I think both TMM and median-ratio are good enough for your approach, I just don't like to throw around the term normalization lightly. IMO normalization implies aligning dsitributions to allow for statistical testing, which those two don't fully do. If there is a huge disparity in the number and magnitude of expressed genes between conditions, or lots of zero inflation, those methods do not produce a clean normalization either.

If I read you correctly, you want to find the gene with the highest expression?
With RNA-seq data from multiple samples, I'd find this by ranking each gene by TPM/FPKM and taking the mean rank across samples. This method will tend to mitigate any large differences in a gene's expression between samples.
As you are comparing (RNA-seq) expression across genes you must take into account transcript length.

OP seeks to find the most highly-express transcript, not gene. The idea with the ranks though should be considered, and one could rank by the TPMs that salmon returns for each transcript. That avoids the use of between-sample normalization. So one would quantify with salmon, take the TPMs, rank transcripts per gene, and then use the transcript with the minimal aggregated ranks per between all samples.


text:
Nucleotide Frequency - HOMER?
I'm trying to find the frequency of UA dinucleotide in a bed file. I want to convert the bed file to FASTA, count all the UA dinucleotides up, then divide by the total length of the FASTA file.
Also, I probably need to compare the value I get above to the background, but I'm not sure how to go about this?
I understand homer should do the counting for me and also create a background distribution for me. It looks like I need to build a basic motif file and give homer a bedfile of locations and motif file to find. However, I'm unsure of how to set this up.
Other suggestions are welcome.

Is it UA per BED interval (so per row), or total entries, or per chromosome/contig? Please provide reproducible in- and output for this question to make sense. I fear otherwise it will quickly get closed for lack of details and clarity.


text:
Solutions to reference bias issue
I've recently learned about the reference bias issue - inability to properly map NGS reads in certain genomic regions caused by the fact that our standard genomic references are linear and do not fully represent the genomic variation present in the population. According to what I read this has multiple consequences: missing structural variants or variants in highly variable regions like HLA for humans.
Now I wonder what is the best way to address this. Of course, "best" is context-dependent but I don't have any specific context so general pros and cons will be highly appreciated. I know about two approaches but probably there are more. The first one is an alt-aware mapping that's used in DRAGEN by Illumina. As far as I understand, they just keep multiple linear alternatives for some regions and try doing a "smarter" mapping for those regions. The second approach that I've found is using graph references where the reference is not linear anymore but rather a graph and mapping is taking into account this graph structure.
Thanks in advance.

Hi @igor, is there any chance you could add some more personal context to this? Bioinformatics frequently offers multiple solutions to the same problem, so answers here are much higher quality if they can be targeted to your specific use case.

Thanks for your comment @gringer.
I don't really have context beyond the described one. Just curious about the approaches that address this issue.

Answer from @throckmorton, converted from comment:
Here is a recent approach to this problem, using multiple linear references:
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02229-3

Slightly older than @gringer 's reference, but still from Richard Durbin himself
"Variation graph toolkit improves read mapping by representing genetic variation in the reference"
https://pubmed.ncbi.nlm.nih.gov/30125266/


text:
How to estimate the phenotypic variation explained by top SNPs from a GWAS study?
I have conducted a large-scale GWAS study and got a few significantly associated SNPs. I used GEMMA with -lmm 1 options to run the GWAS and obtain the beta and standard-error estimates. I want to estimate the percent phenotypic variation explained by each of the significant SNPs. I used the following procedure for estimating the variance explained in R:
fit <- lm (Phenotypic_value ~ SNP_data, data = a)
summary(fit)$adj.r.squared

Here, the datafile a contains three columns namely, sample_ID, Phenotypic_value for each sample, and the biallelic SNP_data. I got a value which is 0.43 meaning 43% phenotypic variation explained by the SNP.
Again, I used another formula which is: 2*f*(1-f)*b.alt^2. Here, f is the minor allele frequency and b.alt is the effect size i.e. beta estimate obtained from GEMMA. This gives me a value of 0.03 meaning 3% variation explained which seems reasonable to me.
My question is that which of the above method is correct? Why do I obtain such a huge difference between the two approaches? Is there any other way to estimate the percent variation explained?
Alternatively, from the GEMMA google group, I have got this formula pve <- var(x) * (beta^2 + se^2)/var(y). But I do not understand how can I obtain the value of var(x) and var(y).
I borrowed the first and the last formula from GEMMA google group discussion. The first one is obviously a normal linear model equation in statistics. The second formula is from link. Also, the data structure is three columns, 1 is with individual, 2 is the genetic information (individual carrying A or T) and 3 is the phenotype.
In the last formula, var(x) is the variance of the genotype vector as they said in GEMMA group and the beta is the effect size estimate and se is the standard error from GEMMA output.
It will be great to receive some feedback on this. Thank you.

Answer from @maximilian-press, converted from comment:
At a guess, the GEMMA linear mixed model is controlling for population structure and other SNPs, which your R example is not (it only fits 1 SNP?), but it's hard to be sure with this information. It would be helpful to see data samples. There is probably more context and usage information available in those resources beyond just the equations. It is especially difficult to evaluate with no terms defined as in your last example (defining the terms in the earlier ones is a good start).


text:
How can I extract sequences from FASTA file using OTUs ID?
I have one FASTA file with sequences and headers, and one OTUs table.
What I want to do is replace the current headers in the FASTA file with just the OTU IDs from the .txt file. The files look like something like this:
FASTA file:
>9d4544fa-322a-4fac-bcf5-f8ca410c7e8a runid=b48dcb2029ef589b8cdb55fe529d93b3baadbc7c sampleid=28012021COVID3 read=57150 ch=310 start_time=2021-01-28T18:23:35Z
TCAGTTACGTATTGCTGGTGCTGAAGAAAGTTGTCGGTGTCTTTGTGTTAACCTTTCTGTTGGTGCTGATATTGCAGTCGCCCCCACCTCCTTGGATCCGCCAGGTTAAACACCCAAGCAGACGCCGAAGATAGAGCGACAGGCAAGTAGGTTAACACAAAGACACCGACAACTTTTCTTCCAGCACCAGCAGCCCGTAACT
>a8d321b9-97f1-42aa-a17a-fa1c4609e60e runid=b48dcb2029ef589b8cdb55fe529d93b3baadbc7c sampleid=28012021COVID3 read=46034 ch=83 start_time=2021-01-28T18:24:16Z
CGTTCGGTGCGTATTGCTGGTGCTGAAGAAAGTTGTCCGGTGTCTTTGTGTTAACCTTTCTGTTGGTGCTGATATTGCGGCGTCTGGCAGGTGTTAACCTGGCCTCGAGAGAGTTTGATCCTGGCTGGGATGAACGCTGACGTGCCTAATACATGCAAGTAG
>6d97cca8-6795-4c59-b4d1-05880735b957 runid=b48dcb2029ef589b8cdb55fe529d93b3baadbc7c sampleid=28012021COVID3 read=46363 ch=227 start_time=2021-01-28T18:24:08Z
GTTCAGTGCATATTGCTGTGGTGCTGCGAAAGTTTGTGGTGTCTTTGTGTTAACGCTTTCTGTTGGTGCTGATATTGCGGCGGCGTCTGCGCGGTGTTTAACCTGGCGGATCAAGGAGGTGTTCAGCCGCACCTTCCGTTTTGACTACGCCATTACGACTGCCCAATCTGTTTTACCCTAGGCCCGATCCTTGCGGTCACGGACTTCAGGCACCCCCGGCTTTCATGGCTTGACAGGCGGTGTGTACAAGGCCCGGGAACATTCACCCGCGCCATGGCTGATGCGCGATTACTAGCGAATCCCAGCTTCACGGGTCGGGTTGCA
>772e0bef-fad4-4767-b0e6-d88e6f16b2ed runid=b48dcb2029ef589b8cdb55fe529d93b3baadbc7c sampleid=28012021COVID3 read=49763 ch=349 start_time=2021-01-28T18:24:10Z
CGTTCAGTTACGTATTGCTGGTGCTGAAGAAAGTTGTCGGTGTCTTTAGGTTAACCTTTCTGTTGGTGCTGATATTGCGGCGTCTGCTTGGGTGTTTAACCTGGCCTCGAGAAGTTTGATCCTGGCTCAGGATGAACGCTAGCTTCTGAGCTTAACATGCAGAGTCCAGGGCAGCATGGAAGAAACTGCTTCTTCTGATGGCGACAACGCACGGGTGCATGCGCGTATCAAACCTGCCTCATA
>2b673c48-69c9-48a4-995b-4ad630e23502 runid=b48dcb2029ef589b8cdb55fe529d93b3baadbc7c sampleid=28012021COVID3 read=71966 ch=124 start_time=2021-01-28T18:24:11Z
TTCGTTCAGTTACGTATTGCTGGTGCTGAAAGTTGTCGGTGTCTTTGTGTTAACCTTTCTGTTGGTGCTGATATTGCGGCTCTGCTTGGGTGTTTAACCTGGCGGATCAAGGAGTGTTAGCCGCACCTTCCGGTACCAACTACCTTGTTACGACTTAGCCCC

OTUs table:
OTU ID  IA1 IA2 IA3 IA4 IA5 IS6
otu500  7   10  3   0   5   0
otu502  7   8   4   0   11  1
otu504  4   5   6   2   4   2
otu505  4   4   6   2   8   3

and my desired output is that the .fasta file would look like this:
>otu500
tggggaatcttagacaatgggcgcaagcctgatctagccatgccg
>otu502
tgaggaatattggtcaatggaggcaactctgaaccagccatgccttggtcaatggaggc
>otu504
gcgaacaggattagataccctggtagtccacgccgtaaacttggtcaatggaggc
>otu505
ttggtcaatggaggcttggtcaatggaggctaccctggttaccctggt
>otu506

There is not enough information here to label the fasta sequences with their OTU assignments. What program did you use to produce the OTU table? The program should have a way to list each sequence's OTU assignment. For example, in the mothur program, the cluster command outputs a list file.

There seems to be some missing info - I can't see any logical mapping between reads and OTU table - do you have an OTU sequence file?


text:
Get frequency of index sequences of fastq file
I recently did some sequencing on a MiniSeq and unexpectedly the Undetermined_S0_L001_R1_001.fastq.gz file contained quite a lot of data (about 40% of total). I would like to know which indices are present in this file and ideally get a frequency of all index sequences in this file. How can I do this?
This is the output of head:
@MN00148:128:000H2YCW2:1:11102:12116:1065 1:N:0:0
GNAGAGGATTCAGGGAGCCAGCGCGCATATATCAG
+
F#FFFFFFFFFFFFFFFFFFA=/A///////////
@MN00148:128:000H2YCW2:1:11102:11354:1065 1:N:0:0
GNACCCGTGCGGCTGGAGAGACCACTAGAAAGAAT
+
A#AFFFFFFFFFFFFFFFFF///////FF=////6
@MN00148:128:000H2YCW2:1:11102:4103:1066 1:N:0:0
ANGTTGGTGGTACATCAGGGACGCGGAGATCTCAG

I do not recall how the formatting of the header exactly is, but the index read sequences should be in there and it is rather trivial to parse and then uniq -c them.

How did you demultiplex? With bcl2fastq? You can generate the index reads like I1

I actually went back to the output folder on the MiniSeq and found this file: DemultiplexSummaryF1L1, which contains demultiplexing information of the whole run. Half-way through this file I found what I was looking for:
### Most Popular Index Sequences

### Columns: Sequence ReverseComplement HitCount

Index
TGACCA  TGGTCA  2455964
CGATGT  ACATCG  2320579
GGCCCC  GGGGCC  2023652
CGCGGG  CCCGCG  1983605
CCGTCC  GGACGG  1641274
CAGATC  GATCTG  1514252
CCGGCC  GGCCGG  1319383
CCGCGC  GCGCGG  1109940
ACAGTG  CACTGT  891344
CCCGGG  CCCGGG  706083
TGCCCA  TGGGCA  168294
CAGCTC  GAGCTG  164091
.... [data continues]

I don't know that you can get them from that file.  I've noticed that when the Miseq is left to run its version of fastq generation software, there are no indices in the undetermined file.  I usually run bcl2fastq myself if I want that information.
They recently came out with a new version of software for the MiSeq, so things might be different now.


text:
Bruker MALDI-TOF bacteria species identification scoring algorithm
Just wondering whether anyone can point me to a research paper which describes how the scoring values are generated by the Bruker software for species identification of bacteria in MALDI-TOF MS. For example, there are countless papers describing the Bruker output score as green (species-level confidence), yellow (genus-level), or red (neither species- nor genus-level confidence). However, none of the papers that I see are discussing how these numbers are actually calculated within the Bruker system. Even the Bruker brochure glosses over the point: "After the acquisition of the spectral data has been completed, a report is generated. The result for each sample is clearly listed under ‘Organism (best match)’ accompanied by the resulting score and an appropriate ‘traffic light’ color scheme." Can anyone point me in the direction of a paper or at least some online information such as wikipedia describing how Bruker came up with these numbers? It's not just an arbitrary number, right?
Another generic description of the scores is in this presentation: "Higher the log (score), higher the similarity between mass spectrum of isolate & the database entry in the reference library."

From the 5th page of the brochure: "Unknowns are then compared to the MSP library using a superior pattern-matching approach. This includes peak positions and intensities, ensuring the highest possible levels of accuracy and reproducibility across the complete range of microorganisms". It seems like it's their secret ingredient and they're not really keen in sharing it with the world :/

This paper as well as this other paper both cite this other paper in their methods section when they talk about the parameter settings for the MALDI Biotyper, and the paper itself doesn't give much detail. That's all I've found so far, and I'm probably showing you something you've seen already

"For identification scoring, the software correlates signal intensities of matched signals of mass spectra. The three scores obtained from such a procedure are multiplied and normalized to a value of 1,000 and its common logarithm. Log scores over 2 are considered ...  whereas log scores over 1.7 generally ... Log scores of 3 are obtained when spectra are matched with themselves."
Sauer, S., Freiwald, A., Maier, T., Kube, M., Reinhardt, R., .., & Geider, K. (2008). Classification and identification of bacteria by mass spectrometry and
computational analysis. PloS one, 3(7), e2843.

For identification scoring, the software correlates signal
intensities of matched signals of mass spectra. The three scores
obtained from such a procedure are multiplied and normalized to a
value of 1,000 and its common logarithm. Log scores over 2 are
considered as reliable identification of a bacterial species (Table
4), whereas log scores over 1.7 generally indicate reliable
identification of bacterial genera. Log scores of 3 are obtained when
spectra are matched with themselves.

Sauer, S., Freiwald, A., Maier, T., Kube, M., Reinhardt, R., Kostrzewa, M., & Geider, K. (2008). Classification and identification of bacteria by mass spectrometry and
computational analysis. PloS one, 3(7), e2843.
Going further, from a random thread on ResearchGate:
https://www.researchgate.net/post/Does_anyone_know_how_is_the_p-value_calculated_by_Bruker_Biotyper_during_MALDI-MS_fingerprinting_of_bacteria

Number of peaks matched in your unknown sample divided by the number
of peaks in the unknown sample, this is given a score (out of 10)
Number of peaks matched in your unknown sample divided by the number
of peaks in the reference sample, this is given a score (out of 10)
Score of relative intensities of matching peaks (out of 10)
Total score out of 1000
Log10 score out of 3
There is no relationship of the score with a p-value nor you can
calculate the statistics because the bruker database is not open and
you can not know the values for what the score was calculated to do
the statistical analysis


text:
Combining dataTracks in Gviz by averaging or subtraction
I would like to average (and maybe also additional operations such as subtraction) 2 or more DataTracks and plot the result using the Gviz (https://bioconductor.org/packages/release/bioc/html/Gviz.html) package.
Say I have:
chr="chr2"
gen="mm10"
from=1000
to=5000

dTrep1 <- DataTrack(range = path/to/rep1.bw, genome = gen,
             name = "bigWig",chromosome = chr)
dTrep2 <- DataTrack(range = path/to/rep2.bw, genome = gen,
             name = "bigWig",chromosome = chr)

Instead of:
plotTracks(list(dTrep1,dTrep2),type="polygon",from=from,to=to)

I would like to:
plotTracks(mean(dTrep1,dTrep2),type="polygon",from=from,to=to)

Or:
plotTracks(dTrep1-dTrep2,type="polygon",from=from,to=to)

See posted also at: https://support.bioconductor.org/p/9137211/
Thanks.

It might help to post this on the BioConductor support forum and add a link to that post here, either as a comment or by editing your post and adding it in the body.