text:
How do I call variants using BCFTools without piping mpileup file?
I need to call variants on a large number of reads using both VarScan2 and BCFTools. Both of these variant calling tools use mpileup file as input to call the variants. VarScan2 command is straightforward, as we can specify the mpileup file as a command parameter. The problem is, BCFTools expects mpileup to be piped from another command and I cannot find an input  parameter to specify a mpileup file that already exists (the one used with VarScan2). It seems horribly inefficient to run the same command twice instead of reusing the file obtained from previous step (especially because generating mpileup takes a substantial amount of time on my computer).
Is there a way to reuse an existing mpileup file with BCFTools?
I am using Nextflow to run my pre-processing pipeline, so could it be possible that I can pipe an input channel with mpileup into the BCFTools's call command?
I also tried doing
cat file.mpileup | bcftools call ... but this does not seem to work.
Edit: cat file.mpileup | bcftools call ... does work, which would be a solution to this question, however, answers offer other, even more elegant, solutions.

I am voting to close this since the OP clarified that cat file.mpileup | bcftools does work, and there was a typo in the first attempt.

No close: There are misconceptions and nuances addressed in the answers that are of educational interest to other readers.

Why not use tee? Please see this post.
You can write mpileup to file while also piping it to another process. For example, samtools mpileup ... | tee file.mpileup | bcftools call ...
the bcftools will run while generating the file.mpileup. After it's done, you can use file.mpileup for your VarScan2.
You can test it with the following:
touch helloworld.txt
ls | tee file.txt | grep helloworld    # output helloworld.txt
cat file.txt | grep helloworld         # output helloworld.txt

Thanks @Law for the help. It turns out cat file.mpileup | bcftools call ... worked after all. It was a silly typo mistake from my side.
But still, I cannot understand how to use tee here. In your example, you have samtools mpileup ... | tee ... | bcftools call .... In my case the first part of that pipe needs to be a file instead of command.

In your original post, you said The problem is, BCFTools expects mpileup to be piped from another command and I cannot find an input parameter to specify a mpileup file that already exists (the one used with VarScan2). From what I understood is that bcftools expect only a piped input. You also said you didn't want to run the command twice of generating the mpileup, hence my answer. You can use tee to pipe mpileup to bcftools while also writing it to a file (aka to stdout and to file). Glad that it was only a typo then.

Bcftools and samtools both follow the very common Unix filter pattern: they read from all input files specified (as non-option plain arguments), or if there are none they read from standard input.
bcftools call is no exception:
$ bcftools call

About:   SNP/indel variant calling from VCF/BCF.
         To be used in conjunction with bcftools mpileup.
         [etc]
Usage:   bcftools call [options] <in.vcf.gz>

So you can write either of the following equivalently:
cat file.mpileup | bcftools call [options]
bcftools call [options] file.mpileup

This assumes that the file is VCF/BCF as produced by bcftools mpileup, and not actually textual “mpileup output” as suggested by the filename. If it is the latter that VarScan2 is expecting, you won't be able to share the input files between the two after all.
See this workflow example, which demonstrates bcftools call both with a pipe and a previously made input file.


text:
mpiblast alternatives?
I am looking for a sequence alignment software capable of pairwise alignment of a large number of protein sequences (10^7).
I tried to go for mpiblast: https://wiki.canterbury.ac.nz/display/UCHPC/Installing+mpiBLAST
but seems their website was down.
Do you know any alternatives to mpiblast which is also highly efficient in protein sequence alignment?
Thanks!

The Diamond protein aligner is a very popular tool for pairwise protein alignment which is much faster than blastp while having comparable sensitivity. Diamond also supports the use of using multiple compute nodes (similar to mpiBLAST). You can view the details of the new Diamond 2.0 release in this recent Nature Methods paper (source of the image as well).

...we aligned all 281 million protein sequences from the NCBI nr database against the UniRef50 database, which consists of 39 million sequences, using DIAMOND (v2.0.7) in ultra-sensitive mode...

Throckmorton thanks a lot for the suggestion! I'll have a look at this aligner now.

source: Stack Exchange

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

sort: votes


text:
Site Mutations in a 3D protein structure to locate the core of mutations
I have a pdb file for a protein (3D structure predicted by a tool) and around hundreds of known point mutations in the protein. I want to visualize the protein with all the mutations to see the hub of or core where the mutations are mostly. Is there are tool that does this? I would give the structure in a pdb file and a list of point mutations myself.

Mutations can have different effects and severity. "Visualising" them is very powerful for discussion purposes, but quantifying them may be better, see this past Q. You can open your protein in PyMOL and select the residues in the command line part (select interested, resi 1+2+3+5) and then show them (show sticks, interested)... but the problem then the format you have the mutations in influences how you make the concatenated indices...

Also this would be showing the wild type residue, but showing all mutations in a single structure makes no sense in most cases.

@MatteoFerla Clustering mutations could reveal lots about the protein providing one has the functional impact of each mutation (this cluster are benign, this cluster causes misfolding etc). Indeed a region of no mutations could indicate constraint.

Clustering mutations could reveal lots about the protein providing one has the functional impact of each mutation (this cluster are benign, this cluster causes misfolding etc). Many tools are emerging to explore this and there are plenty of established tools that can handle this.
VarMap is a tool that is able to map genomic variant coordinates to the best available structure.
PyMol can open PDB files and you can programmatically or manually mutate each position using the mutagenesis tool.

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


text:
Assembling all transcripts for an individual gene? (using single sequence to seed the assembly)
Let's say I have a candidate gene and I believe that in an individual sample, the genome sequence differs from the reference which then interferes with alignment.
Is there a way for me to do a "transcriptome assembly" with a specific read or sequence as a starting point? And only do the k-mer assembly from that single sequence?
I could probably code something myself but I feel like a tool to do this must already exist?
The goal here is to avoid the significant overhead of assembling a whole transcriptome when I'm only interested in one gene.
NOTE: In this case the data is RNA-seq

What did you mean by genome sequence differs? Does it differ because of SNV or indel, or because of different splicing?

Just assume it could be any one or combination of those features.

In essence what you're looking for is targeted assembly. There aren't a lot of tools out there for that that I'm aware of, but on of them is Mapsembler. The general idea is to input a seed region and then let the program perform iterative extension of it. That sort of strategy should largely match your needs.

Yeah this looks like exactly what I would be looking for! Will check it out and get back to you.


text:
creating a dictonary in Python
From the FASTA header I get the following information:
gi|1543850138|ref|NC_040192.1| Desmodium mottle virus isolate UG5 segment DNA-B, complete sequence
['gi', '1543850138', 'ref', 'NC_040192.1', ' Desmodium mottle virus isolate UG5 segment DNA-B, complete sequence']

Is there a better way to create a dictonary as I did with the below script?
from Bio import SeqIO

organismInfo = {}
for seq_record in SeqIO.parse("data/landplantsviruses.fasta","fasta"):
    print(seq_record.description)
    temp = seq_record.description.split('|')
    print(temp)
    organismID = temp[3]
    organismName = temp[4].lstrip()
    print(organismID, organismName)
    organismInfo[organismID] = organismName
print(organismInfo)
```

You may notice that Uniprot and NCBI have different standard: so your solution is not universal, which is why BioPython does not split the description for you —it could know different 'dialects' of fasta files but that is a tedious minefield.
About coding, one write something in a code golf format, but this is problematic as it is hard to interpret and adapt, i.e. it's a bad practice. Having said that if made a list of keys you could use zip with your list of elements as the list of keys may be more intelligible.
header = 'gi|1543850138|ref|NC_040192.1| Desmodium mottle virus isolate UG5 segment DNA-B, complete sequence'
keys = 'db|id|type|id2|wordy.'.split('|')
organism_info = dict( zip(keys, map(str.strip, header.split('|')) ))

Also, for variables in Python, camel case is for Classes, not instances. Even if you have to send data to JS, using a constistent style, i.e. the PIP guideline one, avoids errors.


text:
Edit end column in GFF3 file
I am trying to edit the "end" of my sequences in a GFF3 that contains 10000 sequences (so I do not want to do it manually). I want to add 30 residues to every single fragment. How can I do this?
Currently, I have a file that looks like this:
##gff-version 3
sp_P45452_MMP13_HUMAN   SignalP-5.0 signal_peptide  1   19  0.994300    .   .   .
sp_Q8R2H5_PHLD_RAT      SignalP-5.0 signal_peptide  1   23  0.986089    .   .   .

I would like to edit it so that it looks like this:
##gff-version 3
sp_P45452_MMP13_HUMAN   SignalP-5.0 signal_peptide  1   49  0.994300    .   .   .
sp_Q8R2H5_PHLD_RAT      SignalP-5.0 signal_peptide  1   53  0.986089    .   .   .

Thank you!

That's what awk excels at! Simply set the input (-F '\t') and output (-v OFS='\t') record field separators to tab, and add 30 to the value of the fifth field on lines that don't start with a #:
$ awk -F'\t' -v OFS='\t' '/^[^#]/{$5+=30}1' file.gtf
##gff-version 3
sp_P45452_MMP13_HUMAN   SignalP-5.0 signal_peptide  1   19  30.9943 .   .   .
sp_Q8R2H5_PHLD_RAT  SignalP-5.0 signal_peptide  1   23  30.9861 .   .   .

yep this is the good stuff

Thank you so much! I am new to awk and this was super helpful!


text:
snakemake multiple parameters for multiple input and single output in snakemake. ConbineGVCFs gatk problem
I have written a rule for CombineGVCFs in gatk4. The rule is as follow
all_gvcf = get_all_gvcf_list()

rule cohort:
  input:
    all_gvcf_list = all_gvcf,
    ref="/data/refgenome/hg38.fa",
    interval_list = prefix+"/bedfiles/hg38.interval_list",
  params:
    extra = "--variant",
  output:
    prefix+"/vcf/cohort.g.vcf",
  shell:
    "gatk CombineGVCFs -R {input.ref} {params.extra} {input.all_gvcf_list} -O {output} --tmp-dir=/data/tmp -L {input.interval_list}"

all_gvcf is the dataset for all gvcf files which will be combined. But the problem is I need to add --variant parameter before every input. Command I'm getting right now is as follows
gatk CombineGVCFs -R /data/refgenome/hg38.fa --variant /data/prjna644607/vcf/SRR12165216_HC.g.vcf /data/prjna644607/vcf/SRR12165217_HC.g.vcf /data/prjna644607/vcf/SRR12165218_HC.g.vcf /data/prjna644607/vcf/SRR12165219_HC.g.vcf -O /data/prjna644607/vcf/cohort.g.vcf --tmp-dir=/data/tmp -L /data/prjna644607/bedfiles/hg38.interval_list

The command I want to achieve is as follows
gatk CombineGVCFs -R /data/refgenome/hg38.fa --variant /data/prjna644607/vcf/SRR12165216_HC.g.vcf --variant /data/prjna644607/vcf/SRR12165217_HC.g.vcf --variant /data/prjna644607/vcf/SRR12165218_HC.g.vcf --variant /data/prjna644607/vcf/SRR12165219_HC.g.vcf -O /data/prjna644607/vcf/cohort.g.vcf --tmp-dir=/data/tmp -L /data/prjna644607/bedfiles/hg38.interval_list

How can I add this extra "--variant" tag before every input? I have added it in the get_all_gvcf_list() function. But then snakmake gives me inputfiles not found problem.

Found out the problem. Turns out I can write a lambda function as follows
params:
        extra=lambda wildcards, input: ' -V '.join(input.all_gvcf_list)

and add '-V' before {params.extra} in the "shell" sction. That solves the problem

i am facing a similar issue can you please paste the whole code for your solution?

That is the whole rule/code actually.  get_all_gvcf_list() will give the absolute path of the files in a dataset. You just need to see if the problem really matches with you. If u still need any more help then maybe you can write the entire problem and tag the question with me.

well I came up with the solution. I think it is a little unconventional.
    rule gatk_CombineGVCFs:

        input:
            vcf_dummy = expand("gvcf/{sample}.g.vcf.gz", sample = SAMPLES), # a dummy vcf to connect this rule to gatk_HaplotypeCaller
            fa = "data/genome.fa",

        output:
            gvcf = "gvcf_combined/all.g.vcf.gz",
            #index = directory("gvcf_index")                                 #index generated without seperate output
        params:
            list = "gvcfs.list"
            log: "logs/CombineGVCFs.log"
        shell:
            "ls gvcf/*.gz > gvcfs.list && gatk --java-options -Xmx7g CombineGVCFs -R {input.fa} --variant {params.list} -O {output.gvcf} "
             #Two bash comands running sequentially to make a gvcf.list as during my testing the combineGVCFs tool takes lists as input.

You came up with a solution, not the solution (by your own admission).

The solution to your problem should be described elaborately. Either describe your answer elaborately else add your question separately and then answer

The comments in the code describe the solution pretty elaborately. The problem was in executing the the combineGVCF to which my solution was to add  --variant in the shell and first make a list of vcfs with ls and pass it as a list to combineGVCF.


text:
Variant vs Allele vs SNP
Coming from a CS background. Reading through the wikipedia page, these all sounds like the same thing: Variant, Allele, and SNP.

Variant/Allele/SNP: Some gene locus that differs from the idea human.

For example, if 99% of humans have a T at some locus, and you have a G, C, or A, then you have a Variant/Allele/SNP at that locus.
Is there intuitive way to explain the difference?

Variant = a locus in the genome where there are differences between individuals
Allele = one of the possible bases/sequences that can occur at the variant
SNP = single nucleotide polymorphism, a type of variant. Here the alleles are just single bases, but other variants can be insertions/deletions of one or more bases.

First some basic biology, taking human as an example. Humans have 22 chromosomes (image source):

As you can see (if you squint) in the image, we have two copies of each of the 22. With the exception of the sex chromosomes, X and Y (seen in the bottom right of the image above), you have the same complement of genes on each of the copies of the chromosome. Traditionally, the term "allele" refers to the sequence of one of the two copies.
We inherit one copy of each chromosome from each of our parents. So, if you have inherited a mutation at position $Ψ$ of chromosome 21 from one of your parents, you will be "heterozygous" for that mutation. This means that of the two copies of chromosome 21 in your cells, one will have the mutant allele at position $Ψ$, and the other the wild-type (i.e. "normal", although that's possibly a misnomer) allele.
SNP stands for "single nucleotide polymorphism". This means that it describes a difference of a single nucleotide between the genome you are looking at and the reference genome. For example, if the reference genome has a G at position Ψ of chromosome 21 and you have a T, then this is a SNP.
The term SNP is also commonly used to refer to polymorphic positions. This means specific, single base positions of the genome where variation is common in the human population. Having differences from the reference genome in such polymorphic sites is usually not problematic. The specific use of the term varies a little between subfields. If you're working on population genomics, you are likely to use SNP only for such known polymorphic sites, while others may use SNP to refer to any single nucleotide variant.
Which brings us to variant. This is a blanket term encompassing any and all forms of variation in a genome. SNPs are variants, but so are larger insertions or deletions (indels), copy number variants, large transposed sequences, the list goes on. You can basically think of variant as the catch all term used when you don't want or need to use anything more specific.
To sum up, we have:

Allele: one of the two copies of a genomic sequence present in (diploid) cells.
Variant: a site where the genome you are looking at differs from a reference genome of that species.
SNP: a single nucleotide variant often, but not always, in a position of the genome where a significant proportion of the species' individuals show variation.

Just to add on to previous answers, it is important that allele is one of the possible versions or sequences that can happen at a given position.
Take SNP. Most SNPs have two alleles, meaning that in the population we only see two possible nucleotides at this position. For instance, everybody either has an A or a G at this position. This type of SNPs are called biallelic.
In rare cases, a SNP can have three different alleles, like A/G/C. These are called triallelic. So already you can see that allele is different from SNP. They are actually completely different concepts.
Also, the word allele has been used for a long time. Scientists first used this to refer to a certain version of genetic variation at one place in a genome. This happened even before we know what is DNA. On the contrary, SNP was only discovered after we can sequence human DNA, which is basically in the past 20-30 years. Therefore sometimes it is hard to distinguish the subtle differences between these terms, because they can carry both original meaning based on primitive understand of heredity as well as newer meaning in the era molecular genetics. It might be helpful to think about what these term are referring to in different time and contexts.


text:
Is it possible to do DEG analysis without replicates?
I have three bulk RNA-seq samples (test1, test2, test3) without replicates.
And I noticed that DEG analysis tools such as DESeq2/edgeR/etc cannot be applied for data with no replicates.
So, I just drew a heatmap with expression value (TPM) of genes of interest, using NMF (aheatmap).
But now I want to perform DEG analysis across the three samples. Is there anything I can try for the analysis?
My TPM data (.csv) looks like below.

Thank you so much.

Basically, you can’t do DEG analysis with no replicates in a way that will yield useful results.

According to DESeq2, you cannot analyze your data without having replicates.
In short, you have an experiment where you compare A to B because, in your theory, A is more expressed than B. However, you do the experiment only once. Without repeating the same experiment again under the same circumstances, there's no way to confirm you will have the same results and that your theory is correct.
You can still compare the ratio of A/B but that will not be meaningful in any way since there's not enough statistical power behind it. In your case, doing the ratio of test1/test2, test2/test3 and test1/test3.

Thank you for the answer. Yes, right. So can I do the log2FC of "TPM value" of test/test2, test2/test/3 and test1/test3?

You can't really apply a Log2FC here since these values are not fold-changes. If you mean calculating the Log2 on the ratio of these values, you can do it, but it's only another way of representing your data. You have to be careful about the 0 since you cannot do Log2(0) -> undefined. Obviously, you also can't divide by 0 when performing the ratio. So you have to decide what to do with rows where there are 0. If this answer satisfies your question, don't forget to accept the answer please.

(just quick idea...) How about adding 1 to all the values to calculate Log2FC with TPM values?

In this situation, you can apply the strategy in section 2.12 in the edgeR manual which basically "makes up" a dispersion estimate and then run the normal DE strategy. This is obviously neither reliable, nor publishable but at least it gives you a list of genes you can use for validation. Treat results with care, focus on genes with decent counts (not low on both conditions of the contrast) and select those with large effect sizes (fod changes). Alternatively, you can run vst from DESeq2 and then rank them by fold change, again focusing on genes with decent (non-low) counts. Both is not really reliable, so be careful with the results. As you have no information on the variability of the groups the fold changes are of limited value, but if you are forced to go with these data it is better than nothing.


text:
Defining a bedfile for chromosome X without pseudo-autosomal regions
I would like to define a bedfile for chromosome X, to target only basepairs outside of pseudo-autosomal regions (and going forward other regions, too).
First, I wasn't able to find an already prepared bedfile in any of the online repositories (but let me know if it is available somewhere).
Second, I compiled the bedfile for all basepairs in chromosome X as a starting point (by looking at the latest assembly at genome browser):
chrX    0       156040895       # as defined at https://genome.ucsc.edu/

The genome browser is 1-based and the last base is inclusive, while the bedfile is 0-based and last base is not inclusive, therefore I had to subtract 1 from the bed start.
Afterwards, I wanted to exclude the pseudo-autosomal regions as defined at this page. However, I am not sure if basepair "Start" and "Stop" are 0- or 1-based and inclusive/not inclusive. I might use the same convention used on the genome browser, but I am not sure if that's correct in this case. I wasn't able to find this information on their website or online.
Assuming that the NCBI "Start" and "Stop" field follow the UCSC parameters (1-based and inclusive), the resulting bed (PAR regions excluded) might be the following:
chrX    0       10000
chrX    2781479 155701382

Edit: the section on the centromeres has the "Length" parameter, from which I could infer that "Start" and "Stop" base-pairs are both inclusive. At the same time, the "Stop" base-pair of PAR#2 in chromosome X is base-pair 156,030,895, which is also the last base-pair number if we were considering the UCSC genome browser notation (1-based and inclusive). Therefore, it seems to me clear that the same notation is used, and the proposed bed-file should be correct.

The UCSC Genome Browser renders data visually with a one-based, closed index. The underlying data are often stored and presented as files worked with on the command-line, using zero-based, half-open coordinates.
This difference in indexing can be important when using the website to view data or when using the download tools to grab a dataset for offline calculations, such as files in BED format.
There is a good overview of indexing on Biostars.
To solve this in a general way:
The UCSC Kent Utilities fetchChromSizes tool can be used to retrieve the sizes of chromosomes for your genome assembly of interest, e.g. hg38.
Using this assembly, to build a sorted BED file for set operations:
$ fetchChromSizes hg38 | grep -v '_*_' | awk -v FS="\t" -v OFS="\t" '{ print $1, "0", $2 }' | sort-bed - > hg38.bed

Replace hg38 with your assembly of interest.
To then filter this for chrX:
$ awk '($1=="chrX")' hg38.bed > hg38.chrX.bed

The intervals in hg38.bed and hg38.chrX.bed use a zero-based, half-open coordinate system.
To then get regions of chrX without autosomal regions, you could define a sorted BED file containing those regions, which also use zero-based, half-open coordinates:
$ echo -e 'chrX\t10000\t2781479\nchrX\t155701382\t156030895' > pseudoautosomal.chrX.bed

I am assuming that the coordinates on Ensembl and NCBI use one-based indexing. Converting to a zero-based index requires subtracting one base from the start positions of these intervals.
It would be ideal if these groups were to clarify the indexing they use. I am making this assumption based on how Ensembl defines regions on chrY, which start at position 1 and not 0, as one would expect to see with a one-based coordinate system.
The numbers reported on Ensembl and NCBI are otherwise identical, so I also assume that NCBI's numbers use the same indexing and require the same correction to the start base.
In any case, to calculate the difference between these regions and chrX:
$ bedops --difference hg38.chrX.bed pseudoautosomal.chrX.bed > answer.bed

This is what my output looks like:
$ more answer.bed
chrX    0   10000
chrX    2781479 155701382
chrX    156030895   156040895

These coordinates are in BED format and so use a zero-based, half-open coordinate system, by convention.
To turn these into one-based, closed coordinates — to query signals over exact regions in the UCSC Genome Browser, for instance — you would add one base to the start positions of each genomic interval.
(Edit: Fixed PAR#2 start coordinate per comment.)

Thank you, great answer! I will look more into that, to understand why our results differ. I could spot one mistake on mine, I missed the 10k bp at the end.

Another caveat if you work with a mix of data from UCSC and Ensembl is that they pick different naming schemes for chromosomes. For instance, chrX is X in Ensembl world. If you use set operation tools or awk, Python scripts, etc. to process BED files, the name scheme you use should be consistent over files, either one or the other, in order to get correct results. This is easy to fix, but can also be easy to overlook.

Is there a problem with your command: echo -e 'chrX\t10000\t2781479\nchrX\t155701383\t156030895' > pseudoautosomal.chrX.bed? Why in the first line the start is 10000 (Start in PAR#1, - 1) and the second line is 155701383 (Start in PAR#2, without -1)? To summarize, why you subtracted 1 from the first line "Start" and not from the second? Thanks! (I am assuming you got the "Start" numbers from https://www.ncbi.nlm.nih.gov/grc/human)


text:
Unable to execute DSSP algorithm from BioPython
I have the following settings in PyCharm:

BioPython is properly installed:

dssp_test.py
from Bio.PDB import PDBParser
from Bio.PDB.DSSP import DSSP
p = PDBParser()
structure = p.get_structure("16PK", "16pk.pdb")
model = structure[0]
dssp = DSSP(model, "16pk.pdb")
a_key = list(dssp.keys())[2]
print(dssp[a_key])

Output
C:\ProgramData\Miniconda3\python.exe C:/Users/pc/source/repos/DSSP_algoritm_test/dssp_test.py
Traceback (most recent call last):
  File "C:\Users\pc\source\repos\DSSP_algoritm_test\dssp_test.py", line 6, in <module>
    dssp = DSSP(model, "16pk.pdb")
  File "C:\ProgramData\Miniconda3\lib\site-packages\Bio\PDB\DSSP.py", line 437, in __init__
    dssp_dict, dssp_keys = dssp_dict_from_pdb_file(in_file, dssp)
  File "C:\ProgramData\Miniconda3\lib\site-packages\Bio\PDB\DSSP.py", line 231, in dssp_dict_from_pdb_file
    p = subprocess.Popen(
  File "C:\ProgramData\Miniconda3\lib\subprocess.py", line 947, in __init__
    self._execute_child(args, executable, preexec_fn, close_fds,
  File "C:\ProgramData\Miniconda3\lib\subprocess.py", line 1416, in _execute_child
    hp, ht, pid, tid = _winapi.CreateProcess(executable, args,
FileNotFoundError: [WinError 2] The system cannot find the file specified

Process finished with exit code 1

Why am I unable to execute DSSP algorithm from BioPython?

Just in the offchange, you install the DSSP binary (https://swift.cmbi.umcn.nl/gv/dssp/) as described in the BioPython documentation, right?

@MatteoFerla, Yes. Check the newly added screenshot.

Installing biopython does not magically install DSSP.

This looks like a python problem and not a “biopython”problem.
Check your file structures. Try giving a full path for 16pk.pdb
These might help. https://stackoverflow.com/questions/17658856/why-am-i-getting-a-filenotfounderror
https://stackoverflow.com/questions/18757127/windowserror-error-2-the-system-cannot-find-the-file-specified-is-not-resol

Up-voted to offset a highly incorrect down-vote.

Ok so this might not provide a solution, but here's a couple things someone might find useful;

First, Biopython does not automatically contain DSSP, which you have to download as mentioned in the comments below (from here swift.cmbi.umcn.nl/gv/dssp)
Secondly, DSSP within Biopython, for whatever reason, accepts the path to the file, not the file itself. I see you have that correctly (dssp = DSSP(model, "16pk.pdb")), but just pointing it ou
Third, DSSP originally had dssp executable, but the newer version has mkdssp executable, so you might need to do dssp = DSSP(model, "16pk.pdb", dssp="mkdssp"), so it knows what to associate it with. This was the problem I had, because without specifying the dssp, it is not able to launch DSSP.


text:
Generating sequences from HMM
I encountered this question in a textbook (source: https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.721.2540&rep=rep1&type=pdf) and was trying to workout the solution for it. Was wondering for the last question, wouldn't the probability be 0?

My working:
START -> M1 -> M2 -> M3 -> M4 -> M5 -> M6 -> M7 -> END
Probability Vector [A, T, C, G]
M1 : [0.6, 0.4, 0, 0]
M2 : [0.6, 0.4, 0, 0]
M3 : [0.4, 0.4, 0.2, 0]
M4 : [0.6, 0.4, 0, 0]
M5 : [0.4, 0.6, 0, 0]
M6 : [1, 0, 0, 0]
M7 : [0, 0, 0, 1]
To generate CTAAAAG = 0 * 0.4 * 0.4 * 0.6 * 0.4 * 1 * 1 = 0

[Comments from other post migrated as answer]
The parenthetical statement in the first bullet point says without corrections. In which case, if an unseen nucleotide has a probability of zero, the probability of the sequence is zero. If the unseen nucleotides had the probability of NaN, then the sequence has a probability of NaN, which is a more thorough solution.
The correction would be adding a small iota value to the probabilities to be non zero —fudging against the zero outcome you otherwise find, which may be the result of the training being too small. You have 5 training sequences, so the iota would be <1/5. Let's hypothetically assume that if you were to have picked a 6th MSA sequence it would have had the missing value, so you could set the iota at 1/6. The assumption that you may find that element were you to sample twice as many is more statistically sensible, so iota to 1/10.
I used the letter iota, which is common parlance for something small, but iota is a really bad in maths as it looks like a 1 with English handwriting or in print: you can call the fudge value anything.


text:
How to perform a meta-analysis using data consisting of paired-end and single-end reads generated from Illumina and Ion Torrent?
So basically I have RNA-seq reads that were generated from Illumina and Ion Torrent platforms for yeast species. I have seen an article where they compared liver cells of a rat that were sequenced with these two platforms, but they made use of PORT bioinformatics tool for the normalization of the data. This requires input data to be in the form of the Ensembl database. This database only has Saccharomyces cerevisiae annotation files available on it, whereas I have other species that I need annotation files for as well, which are available in RefSeq format from NCBI.
So, I am still brand new to doing RNA-expression analyses and have never done anything remotely at this scale. Please could anyone give me advice on the steps involved in analyzing the paired-end and single-end reads together as comparable data sets and then carrying out a complete differential expression analysis on this data. Any advice or guidance would really be appreciated, I am honestly not sure how to approach this at all.

link to the various studies? Is single-end IonTorrent only or is there SE data for Illumina too? Why do you have both read types? Are you trying to compare the two technologies, or are you trying to compare different samples with the only available data which happens to be IonTorrent in some samples and Illumina in others? I worry about comparing results confounded by technology, as this could lead to misleading conclusions,.

Hi there @MaximilianPress, thank you for your answer. Only 1 of the datasets is publicly available, it is one of the paired-end libraries generated by Illumina. The idea is to compare the samples from the different studies. There is only single-end available from IonTorrent sequenced samples and only paired-end for the Illumina sequenced samples. I spoke with my supervisors though, they said to just compare the counts data for each independent study to each other. Starting with a PCA plot of the data from the different studies.

I worry about comparing across studies with different technologies, there are likely to be significantly different sample prep, instrument, etc. issues. Sometimes it makes sense to compare across studies in this fashion but usually you want to explicitly control in some way. PCA is likely to just show that the different studies/technology explain most variation, rather than biologically meaningful variables. Again, it is hard to say much without more information.

@MaximilianPress so basically the data is assessing the effects of physical contact between cells of different species versus when these species are grown as mono-cultures. So we are interested in the genes that are differentially expressed as a result of contact with, or the presence of, the other species. I am not sure if this makes any more sense, but my supervisors have also said there obviously will be differences in how these datasets were prepared and so forth. But, they maintain that you should still be able to detect biologically significant information from these comparisons.

Fair enough, I can't really provide guidance here as it sounds like it will necessarily be a little complex and I'm not sure that I understand how the different datasets will be used. Possibly update your question with some of this information?


text:
Why does the SARS-Cov2 genome has letter t
ATTAAAGGTT TATACCTTCC CAGGTAACAA ACCAACCAAC TTTCGAT... is part of the 5'UTR of genome of an RNA virus SARS-Cov-2.
RNA contains letters C+G+A+U.
DNA contains letters C+G+A+T.
Why are there letters T, in the SARS-Cov-2 genome, instead of U?

This was also answered here: https://bioinformatics.stackexchange.com/questions/11353/why-does-the-fasta-sequence-for-coronavirus-look-like-dna-not-rna/11358#11358

We sequence and therefore typically report assemblies as DNA sequences, even if they're actually RNA.


text:
Mapping ICD9 codes in UKBB data to PheCodes using PheWAS package
I am working on the UKBB data and trying to specify cases and controls by converting ICD9 codes to PheCodes. I have input data (icdt) where the first column (f.eid) is patient IDs and the remaining columns are ICD9 codes. I am using CreatePhewasTable function in R for the conversion. The code seems working well, But I am not getting enough 'TRUE' values for the PheCode 250.1 and 250.2, which represent PheCode for T1D and T2D. I am not sure what I am overlooking because I am getting most of the 'FALSE' values for T1D and T2D columns (250.1 and 250.2). I am sharing my R script
library(data.table)
library(dplyr)
library(PheWAS)
icdt <- fread("icd_codes.tab",header=T,sep="\t")
icdt <- as.data.frame(icdt)

for (i in 215:261) {
  icd9 <- data.frame(matrix(nrow = dim(icdt)[1], ncol = 2))
  icd9[,1] <- icdt[,1]
  colnames(icd9)[1] <- 'f.eid'
  colnames(icd9)[2] <- 'DiagnosisCodes'
  icd9[,2] <- icdt[,i]
  icd9 <- icd9[!rowSums((is.na(icd9))),]
  icd.codes <- icd9 %>% mutate(DiagnosisCodes = (paste0(substr(DiagnosisCodes,1,3),".",substr(DiagnosisCodes,4,nchar(DiagnosisCodes)))))
  icd.count <-  as.data.frame(count(icd.codes, f.eid, DiagnosisCodes))
  phenotypes <- createPhewasTable(icd.count, min.code.count = 1, add.exclusions = T, translate=T)
  phenotypes <- as.data.frame(phenotypes)
  assign(paste('Phecodes',i,sep=''),phenotypes)
  myfile <- paste0( "ALL_Phecodes", i, ".txt")
  write.table(phenotypes, file = myfile, sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE, append = FALSE)
}

Thanks


text:
AttributeError: 'FeatureDB' object has no attribute 'strand'
It appers when a gene has a negative strand then I need to reverse the sequence. Unfortunately, I got the following error:
print("1" + '"' + str(fastaIdx) + '"' + '"' + str(mrna_id) + "'" + db.strand())
AttributeError: 'FeatureDB' object has no attribute 'strand'

with the following code:
import gffutils
import pyfaidx
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC

def getPeptideSeq(fastaIDX, db, mrna_id):
    print("Petide Sequence")
    return Seq(getCdsSeq(fastaIDX, db, db[mrna_id]), IUPAC.unambiguous_dna).translate()

if __name__ == '__main__':
    db = gffutils.FeatureDB('annotation.db', keep_order=True)
    fastaIdx = pyfaidx.Fasta('assembly.fasta')
    print("1" + '"' + str(fastaIdx) + '"' + '"' + str(mrna_id) + "'" + db.strand())

        seq = getPeptideSeq(fastaIdx, db, mrna_id)

What did I miss?
Thank you in advance,
Best wishes,

Presumably the field isn't called strand. Can you just print all the available fields in your db object? Also, what input are you giving? Can you provide a minimal, full example that we can try?

gffutils databases don't typically have a single strand as they contain many features with different strands. It is typically the Feature objects within the database that have a strand attribute: https://pythonhosted.org/gffutils/autodocs/gffutils.Feature.html. Was it supposed to be db[mrna_id].strand or something like that.

Thank you, it worked.

If something worked, can you please write it up as an accepted answer? This keeps the StackExchange answer-fishing bots happy.


text:
Help with Limma-model
I'm starting to use Limma to find differential expressed sites within my data and would like to ask if my approach was correct.
The data has a control and a treatment set, each containing data of 48 patients, which were measured before treatment (t0) and after treatment (t1). So I have Control-t0 and Control-t1 as well as Treatment-t0 and Treatment-t1, whereas Control-t0 and Treatment-t0 should not show any effect.
Now I would like to know if there is a significant change in Treatment compared to Control. Therefore, I think, I have to adjust for the baseline (t0), right?
The code below is probably the most verbose option Limma provides, but I wanted to use the long method first before I use the shortcuts.
So would the Interaction-coefficient I defined within my cont-matrix be the difference between Control and Treatment adjusted for the baseline?
Also, I had to turn of adjusting for multiple testing here as my example data is random

#Data.frame with columns representing samples and rows different sites/genes
data_df <- data.frame(replicate(8, sample(100,1000, rep=TRUE)))
colnames(data_df) <- c("Patient_01_t0", "Patient_01_t1",
                       "Patient_02_t0", "Patient_02_t1",
                       "Patient_03_t0", "Patient_03_t1",
                       "Patient_04_t0", "Patient_04_t1")

#Metadata for data_df
meta <- data.frame(Sample= colnames(data_df),
                   ID= c(rep("Patient_01",2),
                         rep("Patient_02",2),
                         rep("Patient_03",2),
                         rep("Patient_04",2)),
                   Time= c(rep(c("t0","t1"),4)),
                   Group=c(rep("Control",4),
                           rep("Treatment",4))
)

#Create levels for model-matrix
time <- factor(meta$Time, labels=c('t0','t1') )
group<- factor(meta$Group, labels=c('Control','Treatment') )
lev=paste(group,".",time,sep="")
lev

design <- model.matrix(~0+lev)
design

#Create contrast matrix with all factors of interest
cont_matrix <- makeContrasts( t1_VS_t0_Ctrl = levControl.t1 - levControl.t0,
                              t1_VS_t0_Treat = levTreatment.t1 - levTreatment.t0,
                              Ctrl_vs_Treat_at_t0= levControl.t0 - levTreatment.t0,
                              Ctrl_vs_Treat_at_t1= levControl.t1 - levTreatment.t1,
                              Interaction = (levControl.t1-levControl.t0)-(levTreatment.t1-levTreatment.t0),
                              levels=design)

#Calculate different expressed sites
fit <- lmFit(data_df, design)
fit <- contrasts.fit(fit, cont_matrix)
fit.eBayes <- eBayes(fit, robust=TRUE)

summary(decideTests(fit.eBayes, adjust.method="none"))

It looks good to me. One think I could suggest is to see the distribution of your data/model (with limma/voom).
After you get the list of your genes with:
top.table <- topTable(fit.eBayes, sort.by = "P", n = Inf)

Maybe you can check your multidimensional scaling (MDS):
plotMDS(your_data, col = as.numeric(group))

Then you can also check with voom
y <- voom(data, model.matrix, plot = T)

I hope this helps

Thank you for your answer. Checking the distro might be very useful, indeed!


text:
Demovir produces an empty output
I am new to bioinformatics and am working on viral metagenomic (virome) datasets and so would like to use Demovir* for doing taxonomic annotations on my viral contigs.
First, I installed all prerequisites in my macOS but I found that Usearch 32-bit version could not run on macOS Big Sur, so I installed all and re-ran demovir.sh my.viral.contig.fasta threads in ubuntu VM. I could therefore successfully obtained trembl_ublast.viral.u.contigID.txt as an input for R, but I would always get an empty output after performing the demovir.R command of demovir.sh in ubuntu VM.
Thus, I took these back to my macOS, used conda to create a virtual environment for R-3.5.1 and installed R using conda install -c r r=3.5.1, then ran Rscript demovir.R trembl_ublast.viral.u.contigID.txt. Again, I got an empty output with only showing the headers of five columns.
I also tried to modify the demovir.sh script by removing all previous parts I've completed and retaining the last command Rscript $DIR/demovir.Rthen ran modified.demvir.sh my.viral.contigs.fasta threadsagain, but still got the same result.
Hopefully anyone can do me a favour to check and resolve this issue.
Many thanks.
Ernie

Quoting the GitHub repository:

"Democratic taxonomic classification of viral contigs to Order and
Family level
When performing metagenomic sequencing of Viral-Like Particle (VLP),
the majority of returned sequences often bare little to no homology to
reference sequences - Viral Dark Matter. Frequently it may be useful
to know which viral taxonomic group these novel viruses are likely to
belong to as this will give information about nucleic acid type, size
and behaviour."

maybe you can open the demovir.R file with Rstudio and execute line by line to check where the failure is. - @zorbax

Thanks for your suggestion. I eventually got a correct final output already.
The key issue was Prodigal would add extra strings after my original contig IDs, so it would cause R not recongnising my viral contigs. Therefore, I added a flag -trunclabels in Usearch (after e-value flag) to truncate my contig names. Finally, I re-ran the demovir.sh in ubuntu VM and it has worked well.
Ernie


text:
How to optimize the number of amplicons ordered for a PCR wet experiment with several genomic ROIs?
Imagine we have an experiment for which we would like to minimize the number of purchased primers. Let's look at the layout of a single primer below:

Assume we have a precomputed set of all necessary amplicons and we are trying to reduce the costs and find the minimum set of amplicons that tile all regions of interest (ROIs) with some flanking.
The algorithm seems to work, but I am not sure it is ideal, so I would like to know how to assert it finds the minimal covering set of amplicons and how to improve the runtime complexity.
First, we define the underline representation of a primer/amplicon (I  use the terms interchangeably), and an auxiliary function for maintaining the same colors among different plotting sessions.
from dataclasses import dataclass
import numpy as np

def rand_color():
   return "#" + "%06x" % np.random.randint(0, 0xFFFFFF)

@dataclass(eq=True,frozen=True)
class Primer:
    left_end:  int
    right_end: int
    left_roi: int
    right_roi: int
    color: str
    level: int

We then create a random set of primers:
N = 25; max_coor = 2_000
min_primer_len, max_primer_len = 200, 500
min_required_flank_len, max_required_flank_len = 10, 50
roi_len = 20
primers = []
for i in range(N):
  left_end = np.random.randint(0,max_coor-max_primer_len)
  primer_length = np.random.randint(min_primer_len,max_primer_len)
  required_roi_flank = np.random.randint(min_required_flank_len,max_required_flank_len)
  left_roi = np.random.randint(left_end+required_roi_flank,left_end+primer_length-required_roi_flank-roi_len-1)
  right_roi = left_roi + roi_len
  primers.append(Primer(left_end, left_end+primer_length, left_roi, right_roi, color=rand_color(), level=i))

The algorithm for finding the minimal primer covering set is as follows:
primers.sort(key=lambda p:p.left_end)

def does_a_cover_b(a,b, flank=0):
  return a.left_end< b.left_roi-flank and a.right_end> b.right_roi+flank

p_tmp = primers.copy()

cover = set()
for p in p_tmp:
  add_to_cover_set = True
  for c in cover:
    if does_a_cover_b(c,p):
      add_to_cover_set = False

  if add_to_cover_set:
    cover.add(p)

Now cover is supposed to be a minimal subset of primers that covers all ROIs. To visualize, we use plot_primer_set twice:
import matplotlib.pyplot as plt
from shapely.geometry import Polygon

def plot_primer_set(primer_set):
  plt.figure(figsize=(4,4), dpi=120)
  for c in primer_set:
    p_all = Polygon([(c.left_end,c.level),(c.left_roi,c.level),(c.right_roi,c.level),(c.right_end,c.level)])
    x,y = p_all.exterior.xy
    plt.plot(x,y, markersize=7, marker='o', linewidth=5, color=c.color)
    plt.yticks(np.arange(0, N, 1))
    axes = plt.axes()
    axes.set_xlim([-100, max_coor+100])

plot_primer_set(primers)
plot_primer_set(cover)

Is this the optimal and correct minimal set?
My complexity is $O(N^2)$, but I expected it to be $O(N \log{N})$, what am I missing?

I ended up implementing the following greedy algorithm that should work fairly good for any binary relation even if it is not transitive nor symmetric.
Assume you have the following binary relation (you may plug in any function you want):
def does_a_cover_b(a,b, min_dist=30, max_dist=130):
  return min_dist < b.left_roi - a.left_end < max_dist or min_dist < a.right_end - b.right_roi < max_dist

Using the auxiliary functions to define the binary relation, draw the graph, and build the graph:
import networkx as nx

def does_a_cover_b(a,b, min_dist=30, max_dist=130):
  return min_dist < b.left_roi - a.left_end < max_dist or min_dist < a.right_end - b.right_roi < max_dist

def draw(G):
  pos = nx.spring_layout(G,scale=2)
  nx.draw(G, pos, with_labels=True, font_weight='bold', node_size=750, width=2)

def build_a_graph(primer_list):
  G = nx.DiGraph()

  for i,p_i in enumerate(primer_list):
    G.add_node(i, data=p_i)

  for i,p_i in enumerate(primer_list):
    for j,p_j in enumerate(primer_list):
      if i == j:
        continue
      if does_a_cover_b(p_i, p_j):
        G.add_edge(i,j)
  return G

we may implement the following algorithm:
p_set = set(primers)
cover_set = set()
covered_primers = set()
while len(covered_primers)<N and len(p_set)>0:
  G = build_a_graph(p_set)
  nodes =sorted(G.out_degree, key=lambda item: item[1], reverse=True)
  node = nodes[0][0]
  absorber = G.nodes[node]['data']
  p_set.remove(absorber)
  cover_set.add(absorber)
  covered_primers.add(absorber)

  for edge in G.out_edges(node):
    primer = G.nodes[edge[1]]['data']
    p_set.remove(primer)
    covered_primers.add(primer)

print(f"Used a cover set of length {len(cover_set)} primers to cover the total primer set with {len(covered_primers)} primers")
plot_primer_set(primers)
plot_primer_set(cover_set)


text:
Why do SAM and BAM use different coordinate systems?
BAM files are, at least as far as I know, simply binary compressed versions of SAM files. They have the exact same information and are used in the same way. Why then does the SAM format use a 1-based coordinate system while BAM uses a 0-based one?
This is the relevant section of the SAM format specification (emphasis mine):

1-based coordinate system A coordinate system where the first base of
a sequence is one. In this coordinate system, a region is specified by
a closed interval. For example, the region between the 3rd and the 7th
bases inclusive is [3, 7]. The SAM, VCF, GFF and Wiggle formats are
using the 1-based coordinate system.

0-based coordinate system A
coordinate system where the first base of a sequence is zero. In this
coordinate system, a region is specified by a half-closed-half-open
interval. For example, the region between the 3rd and the 7th bases
inclusive is [2, 7). The BAM, BCFv2, BED, and PSL formats are using
the 0-based coordinate system.

The two coordinate systems are a very common source of error in bioinformatics, so there must be some reason why these two largely equivalent formats use different systems. For that matter, why do the VCF and BCF formats, whose relationship is equivalent to that between SAM and BAM, also do this? I am guessing there is a good reason and it is equally applicable to both pairs of sister formats.
Note that the samtools manual states that "all position coordinates are 1-based" when giving regions, so it must handle the conversion internally. Indeed, I can confirm that this command returns identical md5sums for both inputs:
md5sum <(samtools view -b sample.bam chr1:36931633-36931871 | samtools view) \
       <(samtools view sample.bam chr1:36931633-36931871)

Arithmetic on a zero-based coordinate system is less complicated than that on a one-based system, so it appears zero-based is often (not exclusively) used for binary data formats, like BAM or bigBed, or text formats like BED, where computers are used to more efficiently calculate lengths of or set operations on intervals.
A more complete answer on the benefits of zero-based accounting is available here, which references the original note (EWD831) on this subject from Dijkstra.
When reading other textual data formats like SAM or VCF, a one-based index is perhaps more "natural" for humans used to counting digits on their hands. If you visually inspect an annotation on the closed interval [3,7] on a line that denotes some abstract chromosome and starts at position 1, it is clear to you where those data points are positioned. Colloquially, you understand that the interval starts on the third position and ends on the seventh. Or, perhaps, people would typically label the initial item in this set the first item, and not the zeroth item.
A one-based index therefore seems to get used often with formats used for visual inspection or presentation of genomic intervals or sequences, such as in genome browsers — or via CLI tools like samtools, when exporting reads to human-readable SAM.
The language used for working with and organizing data can determine choice of indexing, as well. For programs written in C and offshoot languages, including an API like htslib, binary data is typically worked with as chunks of memory accessed via pointers or arrays, and it is convention to iterate over that data in increments that use a zero-based index. Contrariwise, R's roots in FORTRAN and focus on linear algebraic operations (and matrix indexing) may perhaps explain why that language uses one-based indices for its lists, matrices, and vectors.
The designer of the BAM and SAM formats, Heng Li, comments below(*):

I think 1-based is easier for biologists and
mathematicians while 0-based is easier for C/C++/Java programmers.
SAM/VCF is for biologists to read. BAM/BCF is for a small number of
programmers. I wrote the initial specs of sam, bam and bcf2. PS: when
sam and vcf are stored in memory by htslib, everything is 0-based in
accordance with C conventions. bam and bcf are largely direct dumps of
their in-memory representations.

Note: Original comments were unfortunately deleted from answer by someone, but the quotes were taken verbatim at the time they were written.

That's a fair point about human readable formats, but calling SAM human readable is a bit of a stretch, really. Yes, it's text, but it is designed to be consumed by machines, not humans.


text:
ViolinPlot x-axis
I'm trying to set the x-axis in my Violin Plot to go by each patient rather than it being scattered.
VlnPlot(obj, features = c("classical_basal"), pt.size = 0) + NoLegend() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

The above code gives me the x-axis but there is no pattern or order.
I would like the x-axis to be grouped with the HTB, for example ALL of the HTB2867 first, followed by HTB2876, etc...

Seurat does violin plots using ggplot, which orders based on the factor levels. You need to re-factor your grouping category so that the levels matches your preferred order:
# default level sort order is alphabetical
obj[["grouping"]] <- factor(unlist(obj[["grouping"]])) # refactor
VlnPlot(obj, group.by = "grouping", features = c("classical_basal"),
        pt.size = 0) +
  NoLegend() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

How do I re-factor the grouping category?

It is difficult to say exactly what you would like to achieve without an example.
But if you would like an individual axis based on a factor like patient maybe you want to use facet_wrap for the Patient variable like in the example below?
library(ggplot2)

data <- iris

ggplot(data, aes(x=Sepal.Width, y=Sepal.Length, fill=Species)) +
  geom_violin() +
  facet_wrap(~Species, scales="free_x")
```

I have attached an image as well as a better description.


text:
Is there any database about the abundance of proteins in urine of healthy patients?
I have a dataset about adundance of urine proteins from cancer patients. I would like to compare these abundance levels to healthy patients' urine proteins. However I don't have any data with healthy patients. Is there any database about it?

If you have access to the UK Biobank dataset, then they have urine data for healthy individuals.
https://www.nature.com/articles/s41588-020-00757-z


text:
Using SBOL for lineages of recombinant DNA molecules and strains
I wonder about using The Synthetic Biology Open Language (SBOL) for a project involving the development of a web application to track the lineage of recombinant DNA molecules in a laboratory collection. I am wondering whether SBOL can support this use case, and I could really use some help from community members.
Use case
The ultimate goal is to develop a web application to keep an audit trail of recombinant DNA in a laboratory collection, where sequences are traceable through cloning steps up to their entry point into the collection.
For instance, we can imagine we want to document the generation of a new plasmid by combining fragments of two existing plasmids by restriction and ligation. The web application would provide an interface where you can select the input sequences, select the restriction sites and the assembly, and would generate the child sequence. This information would then be stored in a relational database, and would also be exportable to a text file. Essentially, something like the history of a sequence in SnapGene, but using a relational database instead of the file system, and having an open standard to export the information.
Implementation (roughly)
This requires the following steps:

Generating a schema to document genetic engineering steps: The most abstract version of a step would just have inputs and output, but there would have to be a specific implementation for each operation. For the restriction and ligation described above, there would be three steps

Digestion plasmid A: which enzyme was used, and where it cut.
Digestion of plasmid B: which enzyme was used, and where it cut.
Ligation of pieces: Order of the pieces of the assembly

Plasmid A -[restriction]-> Restriction Fragment A
                                                |-[Ligation]->Final product
Plasmid B -[restriction]-> Restriction Fragment B

Building an API that reads the inputs and description of a step and returns the output: For this, I had thought of using the library pydna, since it seems to cover the initial use cases I could think of.

Building a visual interface to generate this documentation: A mockup with minimal functionality can be found here

Developing a web application managing a database based on the schema mentioned above

Will SBOL ever cover this use case?
As far as I understand it, SBOL currently provides the most abstract layer of the schema described in the previous section, supported by the <prov:> tags, so each entity could be related to its ancestors through an Activity that represents the ligation, the restriction. etc. However, from looking at the documentation there seems to be no schema for any concrete implementation of Activity. For example, to describe a restriction operation. Restriction is relatively simple, but for PCR there would be the need to also define a schema for the alignment of primers, with potential missmatches etc. Provenances should also support importing entries from public resources (such as genome databases). These concrete Activity schemas should be readable by an API that would take the inputs and return the output. For cloning steps, this would be performed with pydna or something similar. For sequences imported from public resources, the provenance would have the url of a web request, and the API would request the sequence, and parse it.
From the documentation of prov:Activity: to capture how an entity was derived, it is expected that any additional information needed will be attached as annotations. From this, I wanted to ask whether specific schemas representing the cloning steps will ever be part of the SBOL, or should be provided as external files.
How to approach implementation?
My original plan (before realising about <prov> tags in SBOL) was to do the schema in json, since it's the easiest way to go for a web application. I suppose that the json schema would be easily convertible into xml to eventually be used in SBOL.
I still think it would be best to start with json, since I don't think for the nature of this task I can leverage any of the features of SBOL beyond the <prov> tag. The sequences themselves will be stored in the original import format, and the only new aspect that I will introduce would be the schemas for cloning steps. Once I have a minimal prototype of the schemas and I have tested the application with users, I could come back to these forums and discuss the schema and potential changes that can be made to it, if it is to be included in SBOL.
If you have some thoughts in any of these points, I would really appreciate any feedback.

SBOL 3 is well suited for an application like this, though you may wish to consider some of the specifics (e.g., whether to use a conventional database or an RDF triple store), and has been used for similar (though not exactly the same) workflows in other projects.
To some of your specific questions:

In addition to tracking the provenance relationships, SubComponent relations can be used to track the inclusion of parts within your final build. The sourceLocation field can be used to express the difference between the original part and a post-assembly trimmed version when there are changes in the reaction.
The assembly plans and reactions can be represented not just with activities, but using SBOL Interaction objects.
SBOL Implementation objects can be used for tracking specific aliquots as well as the designs that are intended to be included in them.
If you still need to enhance Activity records, you can either use the base extension capabilities in pySBOL3 or (better) use SBOL Factory to create extension classes via an ontology.
With regards to serialization format, since SBOL 3 is pure RDF, you can use any RDF serialization you want, including JSON-LD.

There's a lot to unpack here, so I would suggest that if you want to dig deeper into any of these specific responses that you ask new questions focused on individual aspects.
Here is an example in prototype work on protocols:
https://github.com/SD2E/paml/blob/main/paml/paml.ttl


text:
NMF clustering not giving results
I'm trying to run the Non-Negative Matrix Factorization(NMF) for my Gene Expression dataset which was originally in matrix form. But it throws errors as given below:
expr <- read_excel("./expr.xlsx")
expr<- as.matrix(expr)

#removing the first column containing gene symbols before clustering
expr_rm <- expr[,-1]
class(expr_rm)
[1] "matrix" "array"

So converted these expr_rm to numeric:
expr_rm <- as.numeric(expr_rm)

#removing negative values for running NMF
eps <- .Machine$double.eps
expr_rm[expr_rm<=0] = eps

estimate <- NMF::nmf(expr_rm, 3, nrun = 10, method = "lee", seed = 123)

Error: NMF::nmf - 10/10 fit(s) threw an error.
Error(s) thrown:
run #1: 'x' must be numeric
Timing stopped at: 1.95 0.14 4.4*

cluster <- cbind(names(NMF::predict(estimate)), NMF::predict(estimate))

Error in h(simpleError(msg, call)) :
error in evaluating the argument 'object' in selecting a method for function 'predict': object 'estimate' not found

result <- list(estimate, cluster)

Error: object 'estimate' not found

For clarity, my expr_rm has row indices as 1, 2, 3, 4,.. which indicates different genes and column headings as samples.
Please suggest a solution for this soon!
Additional information (requested from comments):
>summary(expr_rm) - before conversion
 Sample_C44_0001.genes.results Sample_C44_0002.genes.results Sample_C44_0003.genes.results Sample_C44_0005.genes.results
 Length:145   Length:145 Length:145  Length:145  Class :character

 > summary(expr_rm) - after numeric type conversion
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  0.000   1.808   4.264   4.220   6.333  16.246

If in expr_rm is data.matrix() also, it gives the error like:
> median_centered_data<- expr_rm - med
Error in expr_rm - med : non-numeric argument to binary operator

Answer from @maximilian-press, converted from comment:
Sometimes casting to as.numeric() will convert from two-dimensional matrix to one-dimensional vector. Based on the output you put into comments (which would be helpful in the question), I wonder if this has occurred and caused an issue. For data.matrix, I suggest reading docs for that function- it should have returned only numeric, such that your med variable is likely not numeric, or there is a problem with your data. Possibly the tidyverse tibble created by read_excel might be the issue. I don't know enough tidyverse to say what exactly is going on.


text:
How to copy only certain counts for genes in tsv file to new file in linux
Hi there I have generated a counts table of samples I need to compare by differential expression analysis.
The layout of the counts table is as follows:

Gene_id
Sample_A_r1
Sample_A_r2
Sample_B_r1
Sample_B_r2
Sample_AB_r1
Sample_AB_r2

KLTH001
0
0
2
5
8
3

KLTH002
0
0
5
3
9
7

YL99812
6
4
0
0
12
14

YL99813
9
7
0
0
15
19

So here, for the differential expression analysis, I only want to compare for example the genes starting with "YL-". So in a new file I want the rows including the first row which is "Gene_id" and the sample names e.g. "Sample_A_r1" and then only the rows with genes starting with "YL-" in their names. This will allow me to only compare the changes between Sample_A alone versus when it was combined with Sample_B, or in the table Sample_AB.
If anyone can please give me an indication of how to do this in linux I would really appreciate it. I have seen other similar posts, involving grep and awk I am just not sure of exactly how to do this. I am still quite new to linux. Any help will be appreciated.

Would something like this work?
awk ' $1 YL* ' counts_table.tsv > edited_counts_table.tsv

Would this preserve the header/first line?

Something like this should do job:
# creating a new file with the first line of the original
head -1 foo.txt > bar.txt

# ">>" appends to the file on the right hand side
# "^" dictates that the search term should be at the beginning of a line
grep "^YL" foo.txt >> bar.txt

Thanks very much, will give it a try now

One way using awk:
awk 'NR==1 || /^YL/ { print $0 }' file.tsv

The default action is to print the line/record, so all you need is the condition:
awk 'NR==1 || /^YL/' file.tsv

Thanks I will try this out, it seems the previous command didn't copy all of the occurrences of "YL-" across. Could you maybe suggest how I could just remove any gene names starting with "KLTH-" from my original counts file as well? I would just like to be able to check that I got all records from the original file in the separated files I create.

@Justin1609 You can just negate the pattern match expression: awk 'NR==1 || !/^KLTH/' file. This will just print the first line, or lines that don't start with KLTH.

okay awesome thanks @Steve I appreciate it, will give it a try

Alternatively, you can use grep like so to remove all lines starting with KLTH:
grep -v '^KLTH' in_file > out_file

For example:
echo "Gene_id\nKLTH001\nYL99812" | grep -v '^KLTH'

Prints:
Gene_id
YL99812

Here, GNU grep uses option:
-v : Print lines that do not match.
In the regular expression (regex):
^ : Match the beginning of the line.


text:
Gwas analysis studies using the result of differential analysis
I want to do a GWAS analysis using the differentially expressed genes that I got from a differential analysis. The purpose is the identification of biomarkers.
It's my first time using GWAS, can you please suggest to me a pipeline to use?

If you are looking for biomarkers, why not look at all genes? There will be a lot more information there.

ok thank u,  so we have to use all genes as phenotypes.
it's my first time using GWAS, can you suggest me a pipeline to use.

@MaximilianPress

I am not an expert on recent packages but I would suggest using Matrix eQTL for purposes of the statistical analysis.
As mentioned in the comments I think that it is worth testing all transcripts for eQTL associations. This will additionally give some idea of the false discovery rate. It is possible that testing "candidate gene" differentially expressed genes will lead to a lot of false positives relative to the background level of associations.


text:
Why do we count reverse compliments when counting kmers in a single DNA strand?
When looking for repeat patterns in a DNA sequence, we can look for the pattern and its reverse compliment with up to d mismatches. However, why do we look for reverse compliments if we're analyzing a single DNA strand?
I thought about it and my guess is that this lets us analyze both strands of DNA. For example, if we find the reverse compliment at position x, that means that the other strand of DNA has the pattern at position x - pattern.length.
Is that correct?

When we sequence genomes, the individual reads come from both strands, and that also applies for assembled sequences* (contigs / scaffolds). Therefore for most of the genome assemblies out there, strandeness of individual scaffolds is actually an arbitrary property, and reverse complementing any sequence will generate an equally valid representation (with reverted coordinates).
The coordinates of the reverse complementary k-mer matches depend on the notation. Usually, the notations are still in the coordinates of the reference, for example, blast gives for reverse complimentary strand matches the from subject coordinate greater than the to coordinate. Some others are denoting "the leftmost" position and strand encoded in a different column (e.g. sam).
*perhaps with the exception of chromosomal level assemblies

Do you have any suggested reading for insight on composing a single sequence from both strands? I feel like that's a big mental roadblock for me.

To be completely honest, I don't! I would be looking out for assembly literature (especially of short reads). I will also try to keep my eyes open and let you know if I run into something.


text:
Given two ADN sequence, can there be multiple optimal alignement? And how to optimally find them?
Given two sequences, I know we can determine optimal alignment of sequence with the Needleman–Wunsch algorithm. But is it possible that two alignments have multiple optimal alignment (I guess it would mean that multiple alignement have the same score)?
And if possible, how would you find them via an algorithm? Just looking for a general idea on how to go about it, I can't see how you could do it beside by checking for all possible combinations.

As Needleman-Wunsch is a Dynamic Programming algorithm you do need to complete the entire matrix to find the optimal solution and it is possible to have multiple best alignments.

The algorithm assigns a score to every possible alignment, and the purpose of the algorithm is to find all possible alignments having the highest score.

https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm


text:
How can I create a no-frills view config file for HiGlass?
Is there a basic example anywhere that explains how to create a no-frills view
config to view the contact map for a non-human genome assembly using the
server-based viewer at http://higlass.io?
I believe I have a valid .mcool file for the contact map. But looking at the
description of view config files it seems clear those are not created by hand.
I haven't been able to find any place that describes how to create one.
Any help would be greatly appreciated.

I agree that it would seem nice to have one of these in an easily findable place. There is a minimal JSON in a repo that might be related but it's not at all clear that it's the same thing. But have you looked at the tutorial which seems to have information about getting a quick start on some data? Presumably you have also looked at the viewconf docs.


text:
Statistical or performance test for genomic intervals between 2 samples
Is there an appropriate statistical or performance test to perform when I want to compare genomic intervals between 2 samples. Sample A is the the benchmark of comparing the intervals and sample B is the sample in question. What I am really looking for is not only how many intervals from sample B did indeed overlap with sample A but additionally what intervals from sample B has intervals that don't overlap with sample A and also those intervals in sample A that don't overlap with B and being able to measure the performance of how accurate sample B intervals overlap with sample A. This in some sense, can be translated into the following:

True positive: intervals from sample B that overlap with sample A.
False positive: intervals from sample B that don't overlap with sample A.
False negative: intervals from sample A that don't overlap with sample B.

If I consider the above schema, then I will not be able to do any of the performance and statistical test (ROC curve or accuracy) because it is not possible to know what the true negative and the true negative are needed parameters. Then, what would should I be using to be able to incorporate the above schema and still apply performance and statistical testing.
I highly appreciate any help or pointers in advance.

Could you please give more context around your question? What samples? What gene? Can you provide any examples of what you're interested in?
Abstract questions about hypothetical situations are not easy to answer.

Are you looking for some coding solution to get which A interval fall into B interval ?

Please consider @MaximillianPress's for upvoting/accepting. Accepting an answer does help with site stats

I would recommend looking into BEDtools, assuming that you are talking about genomic intervals, which I think you are based on the question. Possibly a diagram would help.
BEDtools has a suite of statistical tools for evaluating overlaps, including Fisher Exact Tests and also approaches based on the relative distance of 2 sets of intervals. So as long as your genomic intervals are defined by BED files and a genome file, this should be quite easy, e.g.
bedtools fisher -a A.bed -b B.bed -g genome_file.genome
# will print p-values for left, right, and two-tailed tests


text:
If I have data from a full genome sequence, how can I search for a known pathogenic variant of a known gene?
I have a CS background, but am a bioinformatics neophyte
I did a full genome sequence which provided me with ~100GB of files (SNP VCF, Indel VCF, BAM, Indel TBI, SNP TBI, BAM BAI, CNV VCF, CNV TBI, CV VCF, SV TBI, FASTQ R1, and FASTQ R2).
I'd like to search the genome for a pathogenic variant of a known gene. What tools do I need to do this, or where do I start?

Hi and welcome @ensnare
Is only the gene known?  or is the specific variant known as well?  What format is the data you are searching for?

Well, does the vcf have gene annotation in it?  or rsIDs?

I'm going to try and keep this answer from going overboard while still helping you find your feet. The good news is that it sounds like you have started off with some of this work already done for you, and based upon the files you have above, I will try to indicate what you have already and also where you might be able to go with this. The term you should be entering into Google as the generally accepted term for this pipeline is variant calling or possibly variant curation for the downstream steps.

This process starts off with one or more files containing raw observed sequences from a sequencing device (often, but not always, Illumina)  This files will usually be in the FASTQ format, but historically I have seen another format called QSEQ. You will likely never see this format and I feel old just bringing it up, but if you ever encounter it, just convert your files to FASTQ immediately.  (Files entering here: FASTQ, possibly QSEQ if it was generated well before Carly Rae Jepson suggested that we call her, maybe)
Although not a hard requirement, I recommend running something like fastqc or multiqc or some other form of quality control on your reads and making sure that the reads are of satisfactory quality with no anomalous metrics on them. You may or may not decide that you need to do some clean-up on the reads if you determine that there might be adapter contamination or that some reads are partially/entirely garbage quality. (Files entering here: FASTQ, reference genome in FASTA. Files generated: a quality metrics report and possibly a cleaned up FASTQ)
Once you have your raw reads and consider them to be of satisfactory quality, they will need to be aligned to a reference genome. Read alignment is a critical, computationally intensive process that is required for many bioinformatics pipelines, meaning that there are many aligners out there that have been optimized for many different situations. If you are running a variant calling pipeline on Illumina short reads from genomic DNA, a current version of BWA should be a solid choice for you here.  (Files entering here: FASTQ.  Files generated: SAM and/or BAM files with aligned reads)
Now that alignment is complete, we can start getting reads for the real variant calling process by doing any preprocessing needed. The most popular tool(s) for this historically were Picard Tools and GATK, but those two tool sets were used together for so long that they were just begging to be merged. The recent versions of GATK have made it so and Picard tools are now rolled into and called directly from GATK. If you are using an older GATK version, at a minimum you will want to recalibrate the base quality scores to try correcting biases that your sequencer had when assigning confidence scores to individual base calls (they all have this, trust me) and realigning reads around detected insertion/deletion (indel) sites. Newer versions of GATK don't need the indel realignment steps due to the haplotype calling algorithm going even further and reassembling reads around these complex variants. Those are just the minimal preprocessing steps, you may also need to assign read groups if they aren't already in the BAM files and possibly some other stuff (most likely GATK will let you know what's missing as you progress).  (Files in: BAM, reference genome and "gold standard" variants in VCF format. Files out: preprocessed BAM)
Actual variant calling can happen now, and this is likely going to happen using HaplotypeCaller (mutect2, if you're dealing with cancer/somatic changes). If this is your intention, you can look up the latest and greatest practices at https://gatk.broadinstitute.org/hc/en-us/articles/360035535932-Germline-short-variant-discovery-SNPs-Indels- for the most current guidance. (Files in: preprocessed BAM, reference genome.  Files out: variant calls in VCF, possibly an active sites BAM file to visualize reassemblies around active indels)
Filtering of variants should now happen. At this stage, you may be very application dependent, but essentially this process should assign some variants to being high confidence and others to being potentially technical artifacts. I have seen people skip this step either due to lack of interest or concern that it will filter out a real variant of interest, but I don't recommend that approach and have seen more people missing real variants due to being inundated with low confidence variants that should have been otherwise filtered out.
Variant curation can now happen where one takes their high confidence variants and figures out (or usually has a computer program figure this out) which genes (if any) are being affected by the variant, what the chance will be to the protein coding as well as if the variant might affect splicing or some aspect of gene regulation or structure (such as an out of frame indel). It is often also useful here to compare to some database of mutation frequencies to filter out common SNP variants as potential cause for rare conditions.
Identification of candidate mutations based upon the curation results. Identifying the mutations that are in genes that could plausibly control the observed phenotype, are uncommon, and are potentially causing pathogenic changes to the gene or its product (protein).

It looks like you already have VCF files, so at least the initial phase of variant calling has been finished for you already.  I would recommend checking out something like Annovar or Varient Effect Predictor (VEP) that can run over your (hopefully filtered) VCF and assign a predicted gene and effect to a variant. Check both the SNP and Indel VCF files, as indels, while more rare, can have profound effects on a gene.
From your description, it sounds like you wanted to check a few genes in particular, and after the variant effect preduction, you should have the needed information to do this.


text:
Why is bcl2fastq2 taking so long to calculate stats?
Our lab has been using bcl2fastq v2.20.0.422 to demultiplex RNA-seq data sequenced on an Illumina Novaseq machine on a beefy EC2 instance and we've run into the strange problem: namely that while demultiplexing is very fast, generating the stats files are unbearably slow.
In a recent example, we demultiplexed one lane of a Novaseq run (total Fastq.gz files were 464 GB) on a m5.24xlarge instance (with 96 vCPUs) with the data being read and written to the same EFS mount. The demultiplexing took only ~70 minutes, while generating the Stats files took over 4 hours (14,557 seconds to be exact)!
Running with log level INFO, the demultiplexing progressed quickly as expected, all Fastq files were generated, and the bcl2fastq output the message:
INFO: Created InterOp file '"./InterOp/IndexMetricsOut.bin"

and sat for hours using <1% CPU writing to IndexMetricsOut.bin very slowly. Tracking the size of the IndexMetricsOut.bin file showed that it was being continuously written to at the rate of about 4kb per second. It did eventually finish (final size 49M) and bcl2fastq reported that it completed successfully, but it seems crazy to me that this stage could take so much longer than actually demultiplexing almost half a TB of data.
Does anyone have any ideas why the stats files are taking so long to generate, and what I could do about it? Additionally, since the fastq files seem to be completed, is it safe to just move on to alignment without waiting for bcl2fastq to finish?
NOTE: this question was crossposted on Biostars: biostars.org/p/9491612

cross-posted:https://www.biostars.org/p/9491612/

@swbarnes2 That's not against the rules, is it? Since this site is still pretty small, I figured crossposting would increase the odds that someone would be able to answer it

@divibisan it's not against the rules, but it's considered polite to at least make it clear your crossposted it in both posts, otherwise it runs the risk of people wasting their time answering a question which has already been answered elsewhere.


text:
How to find differential expressed genes within pseudotime trajectory with Seurat cluster?

Did anyone know how to find the differential expressed genes within pseudotime trajectory with Seurat cluster? In monocle tutorial, we can use BEAM to find the differential expressed genes between the two branches even with the states difference. But how can I find the differential expressed genes with seurat culster enven with the genotype diffence? Thank you so much!

You probably want approaches such as tradeSeq.
For Seurat to tradeSeq see also this issue.

Not sure if it helps, but there's a tutorial for how to switch between Seurat and monocle here. The basic approach is to convert it to a cell_data_set object, then process that object with monocle:
erythroid.cds <- as.cell_data_set(erythroid)
erythroid.cds <- cluster_cells(cds = erythroid.cds, reduction_method = "UMAP")
erythroid.cds <- learn_graph(erythroid.cds, use_partition = TRUE)

Pseudotime values can be returned back to the seurat object after processing with monocle via AddMetaData:
erythroid <- AddMetaData(
  object = erythroid,
  metadata = erythroid.cds@principal_graph_aux@listData$UMAP$pseudotime,
  col.name = "monocle_pseudotime"
)

Thank you for your reply. So I need to do the code you sent from erythroid.cds <- as.cell_data_set(erythroid) ...to ...col.name = "monocle_pseudotime"), then use the findmarker function of Seurat to find the differential expressed genes? It's that correct?

It's unclear what you're asking, but what you describe here will produce a list of differentially expressed genes based on the cluster definitions of monocle. As far as I'm aware, these clusters don't depend on the pseudotime values.


text:
Dictionary and index of vcf for base recalibration step
I need suggestions in creating index and dictionary with vcf files. For the base recalibration step, I downloaded Homo_sapiens_assembly38.known_indels.vcf.gz from the given link: https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0;tab=objects?prefix=&forceOnObjectsSortingFiltering=false
I am using GATK version 4.2.2.0
When I try to create index and dictoionary of the vcf file, I get the following error: no suitable codecs found. I also tried to check the vcf file using ValidateVariants and I get the same error: no suitable codecs found.
Can anyone suggest which indels vcf file is used for masking sites in the base recalibration step. Has anyone come across this above error and how did they resolve the issue.
Thanks

Can you show which commands are you trying on the file?

I have also encountered this error before in cases where there was something wrong with the VCF file I was supplying.

Double check your command to ensure that you are feeding it the file you believed you were and verify that file is a vcf (gzipped or otherwise)
If you know how, please compare MD5 or SHA checksums of your copy vs. the original to ensure you have an uncorrupted copy.
If this is something you haven't dealt with before, a less effective version would be to check and make sure that the file size on the server is the same as the file size on your machine.
Additionally, you may want to just use the less command to look inside the file and see if anything appears to be corrupted (you may need to decompress it for that purpose).
It may also be worth trying to gunzip or gzip -d commands to decompress the file (assuming you have the drive space).  This will both ensure that decompressing the file on the fly is not causing any grief to GATK (which is unlikely) and that the file decompresses without error (which would also be an indication of corruption).
Lastly, you could try downloading a fresh version.
Please let me know if this helps and you discover the source of the error.


text:
Concatenate multiple fasta files into one file for MLST
I have seven fasta files, one per gene, with more than 400 fasta entries per file, like this:
Input:
Gene1.fasta
>1721_1 gene name
ATG
>10_14 gene name
GTT

Gene2.fasta
>1721_1 gene name
TGA
>10_14 gene name
GAT

Output:
>1721_1
ATGTGA
>10_14
GTTGAT

What I want to do is create one fasta (continuously) file containing seven genes per sample. As you can see, the sample name is not the same before the _ in the fasta headers.
Is there any way to do this with awk or other command-line tools? Thank you very much.

cat Gene*.fasta > all_genes.fasta ?

@TimurShtatland Thank you. But it should contain all genes sequences of each sample. Each sample should be separate files.

I would suggest using samtools faidx to pull the sample's gene out of each file, and then removing the header.  It will ignore anything after whitespace so gene name after space should not be problem. You could iterate across sample names in a bash loop. Of course, it requires that every sample have every gene.

okay. I will try

The concat command from seqkit could concatenate sequences with same ID from multiple files.
seqkit concat [file1] [file2] ...

I think it can be done using python dictionary. You can split the file using the header and then make the fasta id as your key and the sequence as your value. Then using simple loop and conditions you can concatenate them.

Okay. I will try

Hi and welcome to the site! Please try and post complete answers. This isn't really answering, it is only giving a general direction in which an answer could be found.

Assuming your files are small enough to allow you to fit one of them into memory, you can do this:
$ awk '{
        if(/^>/){ name=$0 }
        else{
            if(NR==FNR){seq[name]=$0 }
            else{
                print name;
                printf "%s%s\n", seq[name],$0
            }
        }
      }' file1.fa file2.fa
>1721_1 gene name
TGA
>10_14 gene name
GTTGAT

If you cannot fit the files in memory, you can try this:
$ join -o1.1,1.2,2.2 -t $'\t' <(FastaToTbl file1.fa ) <(FastaToTbl file2.fa) | sed 's/\t//2' | TblToFasta
>1721_1 gene name
ATGTGA
>10_14 gene name
GTTGAT

Note that this approach requires the FastaToTbl and TblToFasta scripts I have posted here: https://bioinformatics.stackexchange.com/a/2651/298.
Or, if brevity is more your thing, the same approach but more compact:


text:
Alternatives to CAP-miRSeq for microRNA processing
I am involved in a Drosophila transcriptomics project in which we have RNAseq data for small RNAs (for microRNAs) in addition to total RNA (for mRNAs).

The total RNA reads have been analysed satisfactorily using the Tuxedo pipeline using TopHat, Cufflinks, Cuffquant and Cuffdiff. (I run this locally on a fast desktop Mac.)
Initially we applied the Tuxedo pipeline to the short RNA reads, but abandoned this because of an incidence (1 in 10) of anomalously low or high replicates that was inconsistent with the bigwig plots from the bam stage, displayed in the UCSC browser.
We there switched to analysing the short RNA reads using the CAP-miRSeq pipeline on Unix, which initially — though very slow — gave consistent results. However recent attempts to process later samples by this method failed as the pipeline can no longer find some dependencies it requires.

I see no point in spending further time on either of these methods.
Question
Has anyone experience of a more modern analysis method for microRNA RNAseq data, preferably one that would run on MacOS?
I have no trouble in running the Tuxedo pipeline and get good results with total RNA. The problem there is not my computer. Unless someone who has analysed microRNA using this package understands the problem I have encountered, there is no rationale basis for trying this again. As regards Cap-miRSeq, it does not run on the Mac, only on Unix. A colleague set it up on a virtual box. We are our support. Is there anyone working on RNAseq of microRNA?

Okay, thanks why not then simply use Bioconductor microRNA, which is under active/recent development?

@M__  Thanks. Didn’t know about that. I’ll have a look at it.

@M__ I had a look at it. It seems to be a tool for identifying microRNAs. However that is not the problem. The reference database for the Drosophila genome includes the microRNAs.

Okay thanks, @llrs will know how to do this, but really I think you are into R-based coding. Between Bioconductor and Seurat there must be a replica of your original program

Answer from @supertech, converted from comment:
This can be frustrating I understand. In my experience it's easier to read the paper and install the modern versions of the tools compatible with your computer (paper was published in 2014). Even more convenient you could open a free account at https://galaxyproject.org/support/account-quotas/ (gives 250GB free space), and do your analysis there. Tophat aligner can be very slow. STAR is better from mapping accuracy to speed; I would recommend that you use that if possible. If you are doing  analysis with large genomes your laptop computer is unlikely to have necessary RAM to run the aligners.

Answer from @m, converted from comment:*

Why not then simply use Bioconductor microRNA, which is under active/recent development?
Looks like STAR has been used for microRNA mapping for ENCODE:
https://www.encodeproject.org/microrna/microrna-seq/
From the ENCODE pipeline developer (@user3274289):
The ENCODE microRNA pipeline uses STAR for mapping and quantification and should be runnable on mac (you might need one on the beefier side though to have enough memory but i have run some workflows on a new macbook pro). We do not have drosophila STAR index available so you need to generate that using the workflow that is provided. See instructions on github.


text:
How do you search for new enzymes that are more stable for handling, immobilization?
Noob here. I get that I should  deduce what characteristics the ideal new enzyme should meet, and then use tools such as PDB and blast to compare to the old enzymes, and use other tools such as pymol, but I have no idea on how to approach this. Where should I look for tutorials?

Please consider the answer and upvoting/accepting. It's good for site stats and its a good answer IMO.

From your title, you are after better immobilisation, i.e. your enzyme behaves poorly when bound.

Poor thermostability: Engineering enzymes for better thermostability, involves predicting mutations that improve its Gibbs free energy (∆G, negative potential = better). Finding mutations that improve a protein by going through the mutational landscape is called scanning.
Rosetta pmut_scan is a command line tool to do this, while Pross is a really nice webapp that also does it.
Poor catalysis: if your attachment is close to the active site, that is a design flaw and you need to attach on the other terminus or within a loop. Circular permutation may be an option.
Aggregation: your enzyme aggregates due to pI close to 7 (see protparam) or has a face that is hydrophobic (see the APBS plugin shipped within PyMOL). Rosetta supercharge is great at suggesting mutations to decrease the pI, but if doing it crudely by changing Gln -> Gln & Asn -> Asp, do check for surface exposure.


text:
Convert VCF file to mpileup.txt
I  am working on an iterative analysis that uses orthologous pipelines that require mpileup.txt files as input for a visualization step. This requires me to convert VCF files to mpileup.txt.
This experiment uses error-prone-pcr and follows a somatic variant discovery workflow that was designed by a 3rd party. I have written and implemented several orthologous pipelines with small differences from the original workflow in order to produce a truth set (variants identified consistently by all of the pipelines). The 3rd party used a workflow that calls Samtools to produce mpileup.txt files as the final output and necessitated the development of a custom visualization script. My orthologous pipelines use mutect2 to produce vcf files that need to be compared to the mpileup.txt files.
If I could convert my vcf files to mpileup.txt files without loss of information it would save me significant time cost associated with writing a similar visualization script for vcf files and prevent bias involved in converting the 3rd parties output to vcf format. There are plenty of examples for using samtools and bcftools to produce an intermediate mpileup file that is then converted to a finalized fitered vcf file. But I can't find much direction for converting a vcf file to an mpileup. Can anybody provide direction for this task?

Cross-posted here: https://www.biostars.org/p/9492155/

Going VCF to mpileup is not really something one does or can do. The mpileup should be generated from a BAM or SAM file or something else that has raw, unfiltered read alignments in it.  The VCF just has data on variant sites, and usually just variant sites that passed a certain likelihood threshold at that.  With the VCF, the best you can do is generate a very limited pileup based upon the VCF-reported counts of reads at each variant locus with no information in between reported variant loci. This would not be a true mpileup and would most likely cause an error with the downstream process or worse, not cause an error and give a result that is divorced from the reality of your sample.


text:
Meaning of order of alleles in VCF GT field for unphased genotypes
Order of alleles in genotype field (on VCFs) can be used to infer haplotypes
However, what about unphased genotypes? I have a VCF file where I only have unphased genotypes:
0/1
1/0
1/1
1/2
2/1

However, I have, for example, both 0/1 and 1/0. Does the order of alleles have some different meaning with unphased genotypes, or genotypes with same alleles (e.g., 0/1 and 1/0) can be considered synonymous?

The fact that they are unphased necessarily means there is no order to them. However, the convention is usually to place the smaller value first, like 0/1 or 1/2 etc. But in reality it doesn't matter.

In the case of unphased VCF 0/1 and 1/0 are exactly equivalent and the order has no meaning.
Throughout the VCF specification, it is always encouraged to use numerical precedence in order to order the genotypes, e.g. 0/0, 0/1, 1/1 would be the preferred representations for diploid bi-allelic genotypes.
However it is not strictly enforced in the VCF specification as for now and you could have 1/0, 2/0, etc, and bcftools will not complain parsing them. The presence of two representation in the same VCF might indicate that the data comes from multiple variant callers or has been through manual editing, since for example variant callers will always use one representation only (numerical precedence usually)
Also, the order does impact how you read some of the sample information fields, for example the AD (allelic depth), for which the values are read following the order present in the GT field.
I had a look at the VCF 4.2 and 4.3 spec, but indeed is not so precise there on these details.


text:
Genotype quality in VCF for haploid calls
I cannot find in the VCF specs what is the default value of Genotype Quality GQ for haploid calls. In such case, is the default the missing value .?
From the specs:

If any of the fields is missing, it is replaced with the MISSING value. For example if the FORMAT is GT:GQ:DP:HQ
then 0 | 0 : . : 23 : 23, 34 indicates that GQ is missing. If a field contains a list of missing values, it can be represented
either as a single MISSING value (‘.’) or as a list of missing values (e.g. ‘.,.,.’ if the field was Number=3). Trailing
fields can be dropped, with the exception of the GT field, which should always be present if specified in the FORMAT
field.

From this I assume that there is no default value, GQ might be dropped from the FORMAT field, correct?

Can you give an example? Why would there be a default quality?

Thanks @terdon, I probably misread the specs. I updated the post with a part of the specs that indicates the GQ field can be dropped from the FORMAT field, if I understood correctly.

Thanks, but I'm still a bit confused. You are already quoting the relevant section that explains that yes, . is what you put when one of the samples or alleles in the file is missing a GQ value.  Are you asking if GQ is a required field for a valid VCF file or if you can have a VCF without GQ?

Sure. My doubt derives from the fact that I haven't understood if, for haploid calls, the GQ equals the QUAL field of the variant, or if it's dropped because we only have one allele.