BLAST(n): No hits found I am currently exploring the BLAST program, just for testing purpose i generated two FASTA files, containing two genes A and B, such that B is just a motif of repeated 'G's that occurs in A.

file A.fna:

>A
ACCAACAGACGGGGGGGGGGGGGAGAAAACCGGCCCCCGCAGGGACGCCACCCAGAGGGAAGGG

file B.fna:

>B
GGGGGGGGGGGG

So by my understanding BLAST should find the local alignment of B in A easily? But the program returns 'No hits found'? Maybe I am doing something wrong here?

I executed the following:

makeblastdb -in A.fna -dbtype 'nucl'
makeblastdb -in B.fna -dbtype 'nucl'
blastn -db A.fna -query B.fna

output of blastn:

BLASTN 2.7.1+

Reference: Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb Miller (2000), "A greedy algorithm for aligning DNA sequences", J Comput Biol 2000; 7(1-2):203-14.

Database: A.fna
           1 sequences; 64 total letters

Query= B

Length=12

***** No hits found *****

Lambda      K        H
    1.33    0.621     1.12

Gapped Lambda      K        H
    1.28    0.460    0.850

Effective search space used: 480

  Database: A.fna
    Posted date:  May 3, 2018  5:17 PM   Number of letters in database: 64   Number of sequences in database:  1

Matrix: blastn matrix 1 -2 Gap Penalties: Existence: 0, Extension: 2.5

Thanks for any suggestions (:

There are three possible problems that come to mind.

Masking
Blast will mask low complexity regions by default. Since your sequence is nothing but Gs, it is a safe bet that it is being masked, so no hits will be found for it.

Score/e-value thresholds
Another source of complication is that even if a match is found, that match will have very bad scores. Both the actual score of the alignment and the e-value will be very bad. Since this is such a simple sequence, it will always score badly.

Word-size
The way blast works is (simplifying a bit) by finding a match for N residues (the word size) and then extending that match if extending increases the score. If your query sequence is shorter than the word size, no match will ever be found.

With all that said, the specific problem in your case seems to be a combination of the word-size and masking. If you run your query setting the word-size to 11 and disabling masking, as shown below, you should get results.
blastn -db A.fna -query B.fna -word_size 11 -dust no

That should return two hits (two because the word size is one less than the query size, so it can align perfectly at two locations):
> A
Length=64

 Score = 23.3 bits (12),  Expect = 5e-05
 Identities = 12/12 (100%), Gaps = 0/12 (0%)
 Strand=Plus/Plus

Query  1   GGGGGGGGGGGG  12
           ||||||||||||
Sbjct  11  GGGGGGGGGGGG  22

 Score = 23.3 bits (12),  Expect = 5e-05
 Identities = 12/12 (100%), Gaps = 0/12 (0%)
 Strand=Plus/Plus

Query  1   GGGGGGGGGGGG  12
           ||||||||||||
Sbjct  12  GGGGGGGGGGGG  23

Note that the default word size for BLASTn searches is usually 11. I don't know why the blast+ version (2.2.31+) I have on my Ubuntu machine didn't use that as a default. The manual says it should. But I got no hits when I didn't explicitly set the word size to 11.


5' and 3' bias in Rna-seq data I'm working with rna-seq samples. I see 5' bias and also 3' bias in the per-base sequence content plot. From this link I see that the bias at the start of the sequences appears to be the result of biased selection of fragments from the library. And this won't effect any downstream analysis. But In all of my samples I also see bias at the end of the sequences. What can be the problem? Samples are strand-specific sequenced with poly-A selection protocol.

In the Transcript coverage profile I see 5'-3' bias for all my samples is: 0.9

Showing the per-base sequence plot for one of the sample: For all the samples I see that 5' and 3' bias.

When I looked at Adapter content this is what I see:

No samples found with any adapter contamination > 0.1%

I expect there was a sequencing problem during the last base, where some of the reagents were running low on the sequencer. This won't pose any real problem, RNAseq aligners like STAR will just soft-clip the last base or two if they're mismatches.

It's common to see a bit of bias toward the 5' or 3' ends in RNAseq, mostly due to whether poly-A selection was done or there was a bit of degradation. The primary thing to ensure is that the bias is similar across samples/groups. If that's the case then you don't have to worry that your downstream analyses are going be affected by some sort of per-group bias.


Merging two files together Is it possible to merge two files together using awk? I have two files with matched  Chrom and pos  column

file1 has the columns

CHROM   POS ID  REF ALT QUAL    FILTER  INFO

file2 has the columns

CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2

I want the final merged file to have these columns

CHROM   POS ID  REF ALT QUAL    FILTER  INFO FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2

I'd recommend using R and either merge or one of the joins. Assuming your data is not too large, data.frame + read.delim should suffice to hold the two datasets. If not, try data.table + fread..


Why does this naive Python paired Fastq reader work? I wrote some throwaway code for ingesting paired-end reads from a Fastq file into memory.

def ingest(stream1, stream2):
    while True:
        record1, record2 = list(), list()
        for _ in range(4):
            record1.append(next(stream1))
            record2.append(next(stream2))
        yield record1, record2

I had initially put the for loop in a try/Except block, anticipating I'd need to explicitly break out of the endless while loop. However, the code seems to work correctly even after I removed it. I even "corrupted" the final read in one of the paired files, and the ingest function performs reasonably (the last read pair is never yielded).

My question is...why does this code work, since I never explicitly break out of the loop?

Here's an example of how the code is invoked.

>>> with open('seqR1.fastq', 'r') as fh1, open('seqR2.fastq', 'r') as fh2:
...     counter = 0
...     for r1, r2 in interleave(fh1, fh2):
...         counter += 1
>>> counter
5000
>>> r1
['@read5000/1\n', 'AGGTCTGGGGACAGGTGTTTTACAGGCTCCTTTAGGGGTGGCAAGTCTCCCGGAAACATGATTAGTGAAGACACGTTAACAGCGAAACTGACAGAGAACCTCAACATTCTCTGGTAGGAAAAGCCAGGATCCTGGCAGAGGAAGCAACGATTCGAGCAGGGAACTGTGCAGGACAGAGAACAATGCGGCCCCAATTTTGATCCAGGGGTAGTTAGCTTAGCTAAGGAGTTCCCCAGTCCCTAAGTGGAGAGGGGATCCAGAACCCTACAATTTATTATTATATCTTAGCGGATATATCTAT\n', '+\n', 'CCCCCGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG@GGGGGGGFAGGGCGGGGGFGGGGGGGGGCGGGGGDGGGGFGGGGGGGGGCGGGGGGGGGGGGGGGCGGGGGGGGGGGGGGGGG@G<GGGGGGGGGGGGFGGFGFF:DGGFFCGGGGCGGFEGGGGGGEGGGFB:GGGDG=FGGGG,GAGGGDGGFGGGCD3GFG+G?G*FGGG<GCGGGGC@CD<FFGF*:+G=FE<;GG,8C7GF=;9D=G<C5C66G6C>7)FD*6G7D)9;*C:)4D5+9>@G*=0):)*;1*504\n']
>>> r2
['@read5000/2\n', 'TGTGGGGTTCTGGATCCCCTCTCCACTTAGGGACTGGGGAACTCCTTAGCTAAGCTAATTACCCCTGGATCAAAATTGGGGCCGCATTGTTCTCTGTCCTGCACAGTTCCCTGCTCGAATCGTTGCTTCCTCTGCCAGGATCCTGGCTTTTCCTACCAGAGAATGTTGAGGTTCTCTGTCATTTTCGCTGCTAATGTGTCTTCACTAATCGTGTTTCCGGGAGACTTGCCACCCCTAAAGGAGCTTGTTGAACACTTGTCACCAGACCTAAGTCCTTAGGGATGAACCTGTGTCTCGGCAC\n', '+\n', 'CCCCBGC9FGCD-GGGFGGGEFEFF@CCGFGGGGFF+DG6,F69+FGG+FGF<FCC;GGEBFGGEDECFGGG6GDG7BFG9GCG9G9GDGGGG<GDFFFGEGC,<GC=GFG:,:A:EE+GD<:E,FGF=FCAGEFF++G,G,<EGGFGFFA>*FG9EGFCFCGC,B>43GFF@E@2>9ED@,568D=A;C+=5@E/F53D>,6F*9@+,E,C=*0+*@56<*9@:/F,1,:*)*14*0CC=*1)*0*,*)0;+)*))/**)**;0;//9)1+)10.()5/)070*)7*C**3/)(2*****\n']

Your code only runs prior to Python 3.7, now it raises RuntimeError: generator raised StopIteration. This is essentially a bug fix described in PEP 479 -- Change StopIteration handling inside generators. You can enable the new behaviour from Python 3.5 via from __future__ import generator_stop.

Note, it's not related to file-handles per se, you can reproduce it with a smaller example with a simple iterator:

it = iter([1,2])

def f(it):
    while True:
        yield next(it)

print(list(f(it)))

You should not write code relying on the bug


BWA: Detecting Variation between Reference Genome and Short-Read Sequences I need to identify all loci in the short-read sequence at which the number of microsatellite repeats (i.e. number of copies of "AA," "GTC," etc.) differ from the reference genome, as well as the locations at which the number of repeats is the same in both the reference genome and the short-read sequence.
I used the Burrows-Wheeler Aligner (BWA mem) to map a high-coverage short-read sequence (obtained from the NCBI's short-read sequence archive) to a reference genome. The output is in .sam format. I have also used a separate program to identify the loci in the reference genome at which microsatellites occur.
I would like to identify all the loci in the short-read sequence at which the microsatellite length and loci differ from the reference genome. Does anyone know any tools or packages I could use to read a .sam/.bam file of a short-read sequence mapped to a reference genome and identify specific loci at which the short-read sequence differs from the reference genome? I am using RStudio and have access to my university's supercomputer clusters.
For info on microsatellites, see here: https://en.wikipedia.org/wiki/Microsatellite#:~:text=A%20microsatellite%20is%20a%20tract,locations%20within%20an%20organism's%20genome.
My end goal is to determine the loci at which microsatellite alleles tend to differ within a species, as well as the loci at which microsatellite alleles tend to be identical between all members of the same species.

I haven't done any whole-genome STR analysis from NGS data myself, but are aware of others that have used lobSTR for this. There's also a recent paper [here] that compares a few different STR analysis packages (i.e. RepeatSeq, LobSTR, HipSTR, GangSTR). Here's the concluding paragraph:

In conclusion, all these tools are built to genotype STRs but
have different strengths and weaknesses. Based on our analysis there is
no clear overall winner. RepeatSeq and HipSTR are the best when
considering genotyping error rate even with low coverage. On the other
hand, GangSTR has an advantage because it is the only tool among them
that can call alleles longer than the read length but shows higher error
rate unless looking at only the enclosed class of reads, which in turn
would lose the GangSTR's advantage of picking up long genotypes. In
addition, GangSTR is the newest tool and so comes with reference files for
different reference builds that are periodically updated according to
the tool's webpage. The correct choice of a tool and the
subsequent filtering depends on the aim of the analysis, and might be
influenced by available hardware resources and time limit for running
tools.

In other words, without more information about the specific problem you're trying to solve, it's not clear what the best tool to use is (this is a common issue with Bioinformatics problems).


Nvidia Parabrick fq2bam pipeline error - No such file or directory I am trying to implement Nvidia-parabrick fq2bam pipeline to process my WES dataset. The reason for opting this pipeline is the huge number of samples (1004) which take many months to process for somatic mutation identification. THe GPU accelerated pipeline will be more helpful to get somatic mutations in a quick manner. After pulling parabrick docker image from here, the following command works fine and return the options :

docker run -it nvcr.io/hpc/parabricks:v2.5.0 -h

But whenever I give the following command to run fq2bam pipeline, I am getting an error of No such file or directory:

docker run -it nvcr.io/hpc/parabricks:v2.5.0 fq2bam --ref ucsc.hg19.fasta --in-fq sample_1.fastq.gz sample_2.fastq.gz --knownSites 1000G_phase1.indels.hg19.sites.vcf.gz --out-bam sample_dups_gpu.bam out-recal-file sample_recal_gpu.txt --tmp-dir /tmp

The error is as follows:
Traceback (most recent call last):
File "/parabricks/run_pipeline.py", line 7, in <module>
sys.exit(PB.pb_main())
File "PB.pyx", line 846, in PB.pb_main
File "/parabricks/pbargs.py", line 1287, in getArgs
return PBRun(sys.argv)
File "/parabricks/pbargs.py", line 392, in __init__
self.runArgs = getattr(self, args.command)(argList)
File "/parabricks/pbargs.py", line 447, in fq2bam
self.addToParser(bsmb_parser_sysgroup, sysOptionGenerator().allOptions)
File "/parabricks/pbargs.py", line 50, in __init__
PBOption(category="sysOption", name="--num-gpus", default=GetNumGPUs(), typeName=int, helpStr="Number of GPUs to use for a run")
File "/parabricks/pbutils.py", line 48, in GetNumGPUs
output = subprocess.check_output(["nvidia-smi", "-L"], universal_newlines=True)
File "/usr/lib/python2.7/subprocess.py", line 567, in check_output
process = Popen(stdout=PIPE, *popenargs, **kwargs)
File "/usr/lib/python2.7/subprocess.py", line 711, in __init__
errread, errwrite)
File "/usr/lib/python2.7/subprocess.py", line 1343, in _execute_child
raise child_exception
OSError: [Errno 2] No such file or directory

As I am using the docker image as it is,it is not clear to me which file is missing. Kindly suggest. Thanks.

My query was answered by NVIDIA Parabrick developer team. The installation instruction are available at https://forums.developer.nvidia.com/uploads/short-url/t5uXsZgdMXx0G2wDcYamFl9sO2x.pdf.


Interpreting quantitive outputs from maximum likelihood phylogenetic trees I ran a calculation in RAxML to determine the majority consensus phylogeny of a maximum likelihood bootstrap (How to show bootstrap values on a phylogenetic tree constructed with RAxML), and I got three output files:

RAxML_bipartitions.output_bootstrap.tre
RAxML_bipartitionsBranchLabels.output_bootstrap.tre
RAxML_info.output_bootstrap.tre

What is the difference between the

RAxML_bipartitions.output_bootstrap.tre file and
RAxML_bipartitionsBranchLabels.output_bootstrap.tre?

In summary,

RAxML_bipartitions.output_bootstrap.tre

Is the only file of interest. The reason this is true in this context is really complicated and you have to understand the statistics of likelihood and how they are interpreted within phylogeny to understand why. This file is simply the final output of a non-parametric bootstrap analysis performed by maximum likelihood.

What on earth is a non-parametric boostrap?

A non-parametric bootstrap is resampling each alignment position with replacement. Thus if we have alignment positions 1,2,3,4,5
A bootstrap resample for 2 replicates might be,

Replicate 1

1,1,3,5,2

Replicate 2

4,2,5,2,1

The ML algorithm will make trees of replicates 1 and 2 and find the consensus between them. If you think about it in any other context a bootstrap replicate is pretty meaningless because it no longer reflects the true biological sequence. Thus information on how the consensus was derived, are not really of interest to us providing we are confident this has been done correctly, viz. RAxML_bipartitionsBranchLabels.output_bootstrap.tre and RAxML_bipartitionsBranchLabels.output_bootstrap.tre

So why is this output of limited use?

There are situations to some investigators this information is useful, but assess the robustness of a tree topology its not needed. The only thing we want is a phylogram (bestTree) with the bootstrap values superimposed on them. We really don't need complicated stuff such as the tree to be represented for example as a polytomy (non-bifurcating tree) because we can just read the bootstraps to make that deduction (values >> 75%). In addition, there is not perfect consensus what boostrap value constitutes robustness, but generally most agree >80% is robust.

What output files have useful information in them?

The information that is important are the files associated with "bestTree", that was the single maximum likelihood tree performed on the intact native sequence. The "info" file for this contains 3 really important parameters:

-lnL ... very important!!
Gamma distribution parameter "alpha",
PINVAR, proportion of invariant sites,

-lnL is the highest log-likelihood (probability) of the phylogeny. It is usually a very small number for which where is an enormous amount of theory over it.

Alpha parameter of the gamma distribution this is the shape parameter of the mutation rate, if it is very low (<1) the distribution of mutations across the alignment is very tight clustered and approximates to a negative binomial distribution. This means some sites don't mutate at all and a small number of sites mutate alot. If it is very large >200 (which is never observed) it approximates to the Poisson distribution, meaning the mutation distribution is randomised across the alignment.

PINVAR this is a straight percentage/frequency and simply means the sites that don't mutate.

How are they calculated?

PINVAR and alpha are not emperically calculated, i.e. if you look at an alignment and say 'no mutations at that position', PINVAR would of course agree but may consider other invariant depending on the phylogeny. These parameters are calculated by maximum likelihood and you can begin to see why the calculation takes so long ... alpha and PINVAR affect the tree topology (which affects -lnL), but the topology affects alpha and PINVAR. Thus, is a multidimensional search of tree and parameter space.

So what stuff do I report in my Results?

Anyway reporting -lnL is good technique and shows the reader you've done maximum likelihood, citing PINVAR and alpha from gamma distribution helps ('Methods' parameters were calculated reiteratively under maximum likelihood). This is only useful for bestTree. The -lnL, PINVAR and gamma's alpha are also calculated for every single bootstrap replicate, but these values are of limited to use, because we have resampled the data, only the consensus tree counts... Obviously presenting the bootstrapped phylogram is extremely important.

Welcome to the technical world of phylogeny!

The amino acid matrix you used BTW .. LG is in vogue right now.

How do I do it?

When I do this stuff its via Biopython and ETE3, I capture the values within the pipeline and don't examine the output files of RAxML because I generate my own.


Download protein structure from gene name I have a list of gene name, I want to download and save each protein 3D structure corresponding to each gene name.

There may also be a misunderstanding, from your previous post. To reiterate, I don't use Ensembl, which is the source of your gene names.

On Genbank there is an entire "Structural" database for every structure that is on PDB. When Genbank hold the protein sequence that was used to generate the structure Genbank only use PDB codes. For example "5LBV_A" is both a PDB code and a Genbank code. So this means PDB entry 5LBV chain A. Thus its easy to download the crystal because you simply ask PDB for 5LBV. To prove it,

efetch -db protein -format gpc -id "5LBV_A" | xtract -pattern INSDSeq -group INSDSeq -element INSDSeq_sequence

OUT
    ircigvsnrdfvegmsggtwvdvvlehggcvtvmaqdkptvdielvtttvsnmaevrsycyeasisdmasdsrcptqgeayldkqsdtqyvckrtlvdrgwgngcglfgkgslvtcakfacskkmtgksiqpenleyrimlsvhgsqhsgmivndtghetdenrakveitpnspraeatlggfgslgldceprtgldfsdlyyltmnnkhwlvhkewfhdiplpwhagadtgtphwnnkealvefkdahakrqtvvvlgsqegavhtalagaleaemdgakgrlssghlkcrlkmdklrlkgvsyslctaaftftkipaetlhgtvtvevqyagtdgpckvpaqmavdmqtltpvgrlitanpvitestenskmmleldppfgdsyivigvgekkithhwhrsgstiggpfeddddkagwshpqfekgggsgggsgggswshpqfek

What I cannot easily do, is go from a generic NCBI Genbank protein code to PDB code.

However, PDB uses Uniprot codes for proteins, so you may be able to use a Uniprot code to get hold of the PDB code, or Ensembl codes may, like Genbank, use PDB accessions.

So in summary perhaps you need to give more detail?


What is the difference between a Bioinformatics pipeline and workflow? I want to understand the difference between pipeline systems and workflow engines.

After reading A Review of Scalable Bioinformatics Pipelines I had a good overview of current bioinformatics pipelines. After some further research I found that there is collection of highly capable workflow engines. My question is then based on what I saw for argo. I would say I can be used as a bioinformatics pipeline as well.

So how do bioinformatics pipelines differ from workflow engines?

Great question! Note that from a prescriptive standpoint, the terms pipeline and workflow don't have any strict or precise definitions. But it's still useful to take a descriptive standpoint and discuss how the terms are commonly used in the bioinformatics community.

But before talking about pipelines and workflows, it's helpful to talk about programs and scripts. A program or script typically implements a single data analysis task (or set of related tasks). Some examples include the following.

FastQC, a program that checks NGS reads for common quality issues
Trimmomatic, a program for cleaning NGS reads
salmon, a program for estimating transcript abundance from NGS reads
a custom R script that uses DESeq2 to perform differential expression analysis

A pipeline or a workflow refers to a particular kind of program or script that is intended primarily to combine other independent programs or scripts. For example, I might want to write an RNA-seq workflow that executes Trimmomatic, FastQC, salmon, and the R script using a single command. This is particularly useful if I have to run the same command many times, or if the commands take a long time to run. It's very inconvenient when you have to babysit your computer and wait for step 3 to finish so that you can launch step 4!

So when does a program become a pipeline? Honestly, there are no strict rules. In some cases it's clear: the 10-line Python script I wrote to split Fasta files is definitely NOT a pipeline, but the 200-line Python script I wrote that does nothing but invoke 6 other bioinformatics programs definitely IS a pipeline. There are a lot of tools that fall in the middle: they may require running multiple steps in a certain order, or implement their own processing but also delegate processing to other tools. Usually nobody worries too much about whether it's "correct" to call a particular tool a pipeline.

Finally, a workflow engine is the software used to actually execute your pipeline/workflow. As mentioned above, general-purpose scripting languages like Bash, Python, or Perl can be used to implement workflows. But there are other languages that are designed specifically for managing workflows. Perhaps the earliest and most popular of these is GNU Make, which was originally intended to help engineers coordinate software compilation but can be used for just about any workflow. More recently there has been a proliferation of tools intended to replace GNU Make for numerous languages in a variety of contexts. The most popular in bioinformatics seems to be Snakemake, which provides a nice balance of simplicity (through shell commands), flexibility (through configuration), and power-user support (through Python scripting). Build scripts written for these tools (i.e., a Makefile or Snakefile) are often called pipelines or workflows, and the workflow engine is the software that executes the workflow.

The workflow engines you listed above (such as argo) can certainly be used to coordinate bioinformatics workflows. Honestly though, these are aimed more at the broader tech industry: they involve not just workflow execution but also hardware and infrastructure coordination, and would require a level of engineering expertise/support not commonly available in a bioinformatics setting. This could change, however, as bioinformatics becomes more of a "big data" endeavor.

As a final note, I'll mention a few more relevant technologies that I wasn't able to fit above.

Docker: managing a consistent software environment across multiple (potentially dozens or hundreds) of computers; Singularity is Docker's less popular step-sister
Common Workflow Language (CWL): a generic language for declaring how each step of a workflow is executed, what inputs it needs, what outputs it creates, and approximately what resources (RAM, storage, CPU threads, etc.) are required to run it; designed to write workflows that can be run on a variety of workflow engines
Dockstore: a registry of bioinformatics workflows (heavy emphasis on genomics) that includes a Docker container and a CWL specification for each workflow
toil: a production-grade workflow engine used primarily for bioinformatics workflows


GSEA- parsing error for .cls file I am creating a .cls file for my GSEA analysis and see the below error. Can anyone suggest a solution?

I create the file using TeXworks. The file looks like this:

37 2 1
#CD69neg CD69pos
CD69neg CD69pos
CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos

Error Details

---- Full Error Message ----
There were errors: ERROR(S) #:1 Parsing trouble edu.mit.broad.genome.parsers.Par ...

---- Stack Trace ----

of exceptions: 1
------Number of items found in cls data 10 is not equal to the number specified on the header line 37------

edu.mit.broad.genome.parsers.ParserException:
Number of items found in cls data 10 is not equal to the number specified on the header line 37 at
edu.mit.broad.genome.parsers.ClsParser.doSanityChecks(ClsParser.java:388) at
edu.mit.broad.genome.parsers.ClsParser._parse_genecluster_style_categorical(ClsParser.java:252) at
edu.mit.broad.genome.parsers.ClsParser.parse(ClsParser.java:226) at
edu.mit.broad.genome.parsers.ParserFactory._readTemplates(ParserFactory.java:342) at
edu.mit.broad.genome.parsers.ParserFactory.readTemplate(ParserFactory.java:293) at
edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:753) at
edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:726) at
edu.mit.broad.genome.parsers.ParserWorker.doInBackground(ParserWorker.java:52) at
javax.swing.SwingWorker$1.call(SwingWorker.java:295) at
java.util.concurrent.FutureTask.run(FutureTask.java:266) at
javax.swing.SwingWorker.run(SwingWorker.java:334) at
java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1149) at
java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:624) at java.lang.Thread.run(Thread.java:748)

You seem to have forgotten the space after the #. It should be as in here (Just for future readers know where to find the specs). Also it should be a single line (Maybe this error was introduced while pasting here or it has been a duplicate because there are 39 samples in the 3 line and they should be 37).

Either adding 2 more samples to the initial count and removing the newline of the third line:

39 2 1
# CD69neg CD69pos
CD69neg CD69pos CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos

Or removing the third line:

37 2 1
# CD69neg CD69pos
CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos


How to create artificial or dummy phenotype data I have whole genome sequence (human) data (SNP data), but don't have phenotype data ready for any traits yet. Hopefully, I will get real phenotype data soon. To run GWAS (genome-wide-association study) or any other analysis, i want to create artificial/dummy data for any trait, let's say human height. Could anyone please recommend what is best option? Does plink generate any dummy phenotype data? or any other programme?

Using wakefield package it is pretty easy, here is an example from GitHub page. This will create 500 dummy variables for id, race, age, etc. :

library(wakefield)

myDummyPheno <- r_data_frame(
  n = 500,
  id,
  race,
  age,
  sex,
  hour,
  iq,
  height,
  died
)

# result
myDummyPheno

# # A tibble: 500 x 8
#    ID    Race        Age Sex    Hour           IQ Height Died
#   <chr>  <fct>     <int> <fct>  <S3: times> <dbl>  <dbl> <lgl>
#  1 001   Bi-Racial    27 Male   00:00:00     93.0   71.0 F
#  2 002   White        28 Female 00:00:00    104     64.0 T
#  3 003   Black        28 Male   00:00:00    111     67.0 T
#  4 004   Hispanic     35 Male   00:00:00     90.0   67.0 T
#  5 005   Hispanic     21 Female 00:00:00     94.0   70.0 T
#  6 006   Black        35 Male   00:00:00     74.0   73.0 F
#  7 007   Black        27 Male   00:00:00    102     68.0 T
#  8 008   White        28 Female 00:00:00     95.0   70.0 F
#  9 009   White        30 Female 00:00:00     80.0   71.0 F
# 10 010   White        30 Female 00:00:00    114     66.0 T


How do I rename output folder created using WDL's glob function, with cromwell execution engine? I made a wdl (genomics pipeline language by the Broad Institute) task for fastqc that outputs the zip files, of which there will be as many as reads that the sample I'm working on at that moment has. We create the output using glob:
  output {
    Array[File] fastqcArray = glob("*.zip")
  }

I would like to give the folder a name that is representative of what it contains, for example fastqc, rather than glob-xxxxxxx, so that it's clear what it is in my task's output folder. Is this possible, and if so, how can I achieve it?

Ok, thanks to Ruben Vorderman on the openwdl github page, I found a solution that does not involve renaming the glob folder (which I believe is not possible?). The only way is to write out the names of the files and pass them explicitly to the Array[File] output. This is how I did it (simplified), in case it helps anyone else in the future:
task fastqc_task {
  File forwardReads
  File reverseReads
  String fwdName = sub(basename(forwardReads),".f.*q.*$","")
  String revName = sub(basename(reverseReads),".f.*q.*$","")
  command {
    mkdir -p ./qc/fastqc_zip
    fastqc ${forwardReads} ${reverseReads} --outdir ./qc/fastqc_zip
  }
  output {
    Array[File] fastqcArray = ["./qc/fastqc_zip/${fwdName}_fastqc.zip", "./qc/fastqc_zip/${revName}_fastqc.zip"]
  }
}


Chromosome information for miRNA Apologies if this question is previously asked. I have a BAM file alignment to hg19 and miRNA gff3 file from mirbase for annotation of BAM file. I want to extract information in the following format (e.g.):
miRNA        chr    chr_start    chr_end    expression_count    sequence
hsa-mir-155  chr2       ...         ....         1020            ATGC....

I know that I can use htseq-count for BAM annotation and for expression counts but how can I get their chromosomal location. The reason for doing this specific way is that I want to show them in a CIRCOS plot.
Just counting is not the goal, their chromosomal location also needs to be extracted. So just subread/htseq-count is not sufficient for extracting genomic location information.
Previously, I had the same idea about processing the sam records. If I want to annotate chromosomal information with miRNA using .gff file, I have to check for each line in .sam records whether their chr_start and chr_end lie in the bandwidth of start and end of any miRNA in .gff file. This process can be very heavy and slow because .sam records are of more than 5 GB. That's why I want to know if there is any extremely fast and generalized way to do it.
I want to avoid processing sam records in python because my experience with extremely heavy data like fastq or sam in python is not good. I would rather prefer awk/bash or any low-level language.
Thanks.

you can do samtools view pathtobamfile, this will convert the file to .sam format and send it to stdout, the columns can then be parsed out which contain the information you are after.


Which of these enzymes is a better target for virtual screening? I am a MSc student in Applied Biochemistry. For my Masters Thesis, I will be doing some research about finding new potent inhibitors for a specific enzyme of my choice using molecular docking (virtual screening).

I have selected 5 enzymes that have all been approved by my advisor and now the choice is for me to make. I am very interested in neurological diseases as you can see:

Kynurenine 3-monooxygenase (KMO)

Glycogen synthase kinase-3 (GSK3)

Dopa decarboxylase (DDC)

HMG CoA reductase

Phospholipase A2 (PLA2)

Knowing that I plan to publish my work. I would be really really grateful, if you guys could help pick :)

The ultimate choice stems from your individual interests and not a scientific reason... but having said that some pointer of issues may be handy.

I am assuming that you want to make an inhibitor that messes with the activity and not a ligand that binds somewhere and you then plan to add a moiety to make it a protac to get it degraded or to make it fluorescent to mark it. If not, ignore any talk of activity.

Crystal

is a crystal structure available?
does the structure have the cofactor if any?
is a structure available (itself or homologue) with the native metabolite bound?
how is the active site? In the case of many cytochrome P450 enzyme the entranceway is really deep and is guarded by a flappy gate, both of which confer a great deal of specificity regardless of the residues in the active site, which are often just water coordinating. Your first enzyme, is that a NADPH/FAD or a cytochrome enzyme?
What is the mechanism and is there a classic type of inhibitor chemistry? I.e. catalytic triad with a cysteine nucleophile is a nice target.(1)

Methodology

Wet lab: A problem with computational chemistry/biochemistry is that someone needs to make the protein and physically test your compounds. Are these easy to express? This may sound like not your problem, but a lot of issues can come from this latter stage.
is there a close homologue you want to avoid targetting? If so, you'll have to do twice the docking (not a big deal and potentially more interesting). (2)
is the method you plan to use okay with covalently bound ligands? Dopa decarboxylase sounds like a PLP enzyme. If it is, making an inhibitor using classic functional groups that nucleophilically attack the pyridine ring is very straighforward, but not all docking methods are okay with covalent ligands.
Are any hits likely to follow Lipinski's rule of five? If you have a very charged active site, your hits will be metabolite-like (polar) and not drug-like (horrendously insoluble).

A minor thing is to make sure that nobody else is working on it or you may be left behind. Industry is all cloak and dagger, but the academic SGC has a list of [target enabling packages)[https://www.thesgc.org/tep] and their future ones (nominations).

Physiology

This was most likely in your lit search. But I am double checking.

You have chosen these enzymes because their unchecked activity is disease-causing, right?
And they do not have any other roles in other tissues?
Abolition of activity does what you want?

Other hits

Lastly, but not least. See which of these have many paper about barely active hits: this should mean it is hard. For example, bacterial alanine racemase has been a holy grail for a while, but progress is slow as it has been described as a "ship in a bottle" due to its long and unforgiving entranceway.

Addenda

A catalytic triad is acid-base-nucleophile: if you have an electrophile with a leaving group (including internal ones, such as epoxydes) you can virtually screen more complex small molecules that resemble your ligand, but with the difference that there is a leaving group, resulting in an irreversible inhibition. Whereas competitive inhibitors require nanomolar affinity to outcompete the natural ligand, an irreversible inhibition can even slightly worse and will work.
The problem of homologues.

Paralogues (homologous genes in the same organism): you want to target protein 1, but you do not want to target protein 2. This is a common problem with kinases. You therefore need to find a hit that works better in one than in the other, which with a big of rational inference is easy.
Homologues (homologous genes in different organism). You want to kill species A (say a fly), but not B (say a human). Same thing.

Twice the docking. It is more of a challenge, but makes for a better spin on the story. There is an explosion of people doing virtual screens that do not go anyway afterwards, so it is getting harder and harder to publish studies without extensive empirical backing. Hence the need for an interesting angle.


Visualization tools for 3C/Hi-C long-range interaction data? I have identified a genomic region (CpG island) that is enriched for a characteristic of interest in the vast majority of my samples.

After viewing it in the UCSC genome browser, I have confirmed that it does not overlap any annotated gene promoters, gene bodies, or retrotransposons, but am convinced that it is functionally interacting with something long range.

Is there a publicly available visualization tool where I can input my genomic region and get a list of topologically associated domains from 3C or other chromatin-chromatin interaction data?

I would also recommend two very recent Hi-C visualization frameworks (with some public data available in both): HiGlass and JuiceBox.


Deciding which samples go in which batch I have 370 samples to sequence, we probably will end up using only 96 samples per run (due to the barcode with the primers we'll use). This means running 4 batches. To minimize the batch effect I need to decide which samples are sequenced in which run.

I have several variables from these samples, some of them are categoric (sex, disease status, treatment stage, ...) and some of them are numeric (age, clinical score...).

To minimize the batch effect I wouldn't group all samples from the same category in one run, and I would distribute as uniformly as possible the samples according to the numeric variables. However, I could find many papers about adjusting, taking into account, but I couldn't find a paper about software to help designing the runs.

Is there any software to help designing each batch?

If there are several I would like to point out that I would define best as the one with easier user interface

I found the OSAT package which does almost what I want.

However it is specific for Illumina and for a limited number of plates: "plates that hold 2, 4, 8 IlluminaBeadChip chips and have 24, 48, 96 wells, respectively."
It doesn't take in consideration the numeric values:

library(OSAT)
inPath <- system.file('extdata', package='OSAT')
pheno <- read.table(file.path(inPath, 'samples.txt'),
                    header=TRUE, sep="\t", colClasses="factor")
pheno$AgeGrp <- runif(576, 0, 100)
gs <- setup.sample(pheno, optimal=c("SampleType", "Race", "AgeGrp"))
gc <- setup.container(IlluminaBeadChip96Plate, 6, batch='plates')
gSetup <- create.optimized.setup(sample=gs, container=gc, nSim=1000)

The minDiff package reported in Stats.SE.

It handles both numeric and categorical data, but it can only optimize for two nominal criteria.

I wrote experDesign which does everything I wanted

It can handle both numeric and categorical data for an arbitrary number of variables. It is not restricted to an even number of samples per batch.


Which gene I should select from this qqplot I have a qqplot of my whole genome sequencing data; A plot is for showing possibly significant driver genes. I tried to read about qqplot though but people only say about the skewedness while I want to know from these two genes which one are more likely to be driver

Given that TP53 is the most significant and is already known to have driver mutations in cancer it would seem to be the logical choice.


Normalizing RNAseq for PCA and CCA Usually the expression data is transformed to log space using either RPKM, FPKM or CPM, this is required when looking for differential expression because the data is tested against the normal distribution(limma) or the negative bionimal distribution (DESeq2).

In weighted gene co-expression analysis (WGCNA) there is also the recommendation to normalize the data (see FAQ 4) but says "RPKM, FPKM, or simply normalized counts doesn't make a whole lot of difference for WGCNA" because it is based in correlations.

For PCAs the data should be on the same scale to avoid missing the scale as a factor. See this for more references. But should it be normalized too to look like a normal distribution using logarithms?

How does this extend to other methods like CCA? Should the data look like a normal distribution for CCA? As it is based on correlation I don't expect much changes, so I am using raw counts (RSEM output), but the examples I've seen use normalized data.

This question comes after reading a comment

I've given an answer about using how I normalise for PCA here. PCA works best if the data are normal-ish, but as models go it is fairly robust to fairly large deviations from normal. Transforming is not essential, but if data can be easily transformed into something that looks a bit more normal, then conclusions derived from looking at the PCA should be a bit more robust.

I'm not familiar with canonical-correlation analysis (CCA), but because I see sums and square roots in the formulas, I would assume it has similar data cleaning requirements to PCA. Unfortunately the R documentation for cancor doesn't have any information about any assumptions made on the data.


Retrieve overlap widths with Bioconductor I want to retrieve all overlaps between two sets of genomic intervals (GRanges) using Bioconductor. I am interested in the widths of the resulting overlapping ranges.
In cases with discontinuous ranges, intersect() works, but I am interested in the case where both sets contain contiguous ranges (no gaps).
I am using the findOverlaps function as in the example below:
gr1 <- GRanges(seqnames="chr2", ranges=IRanges(c(1, 11, 21, 31), c(10, 20, 30, 40)))
gr2 <- GRanges(seqnames="chr2", ranges=IRanges(c(1, 5, 15), c(4, 14, 30)))
hits <- findOverlaps(gr1, gr2)

This code yields the following intersections in a "Hits" object:
Hits object with 5 hits and 0 metadata columns:
      queryHits subjectHits
      <integer>   <integer>
  [1]         1           1
  [2]         1           2
  [3]         2           2
  [4]         2           3
  [5]         3           3
  -------
  queryLength: 4 / subjectLength: 3

But I could not find a way to extract the intersections widths from the resulting object.
It this the wrong approach ?

Note: Using intersect(gr1, gr2) in this case yields a single contiguous overlap range spanning the whole chromosome, but I want each overlap to be returned independently, just like findOverlaps does.

Whilst intersect is sort of doing what you want, it merges adjacent overlaps together which is not what you want I suppose. There doesn't seem to be dedicated Grange function to do this out of the box, but you could go through all overlapping pairs resulting from findOverlap and calculate it's widths using pintersect:

punion, pintersect, psetdiff and pgap are generic functions that compute the element-wise (aka "parallel") union, intersection, (asymmetric!) difference and gap between each element in x and its corresponding element in y.

"pintersect"(x, y, resolve.empty=c("none", "max.start", "start.x"), ...) # x, y   IRanges objects.

hits <- findOverlaps(gr1, gr2)

width(pintersect(gr1[queryHits(hits)],
                 gr2[subjectHits(hits)]))
[1]  2  6  4  6 10


Obtaining HGDP project data in fasta format I need to obtain sample data from modern humans in fasta format. I just need some megabytes of data from every individual. I actually use a script that obtains the cram file from here (ftp.1000genomes.ebi.ac.uk) and then processes it to obtain the fasta file.
The problem is that cram files are large, slow to download and slow to process. It takes days to get the samples.
Is there a better way to get these samples in fasta format?

The script already makes use of samtools to retrieve only the part of the bam file it needs but doesn't help much. Cram files are still gigabytes large for only a few megabytes of data that I need.

I have the same problem with data from the 1000 genomes project.

You can download HGDP data in FASTQ format here: https://www.internationalgenome.org/data-portal/data-collection/hgdp


Reference genome for allele specific expression We are trying to sort out a pipeline for doing allele specific expression. Our plan is to call SNPs from RNA-seq data and combine with known SNP annotations. A well known problem in ASE is reference bias, where reads are more like to map if they have the reference sequence than the alt sequence.

We have come accross two possible ways to deal with this. In many papers (e.g. Zhuo et al, Sigurdsson et al) a new reference is generated, substitute N bases where there are SNPs using FastaAlternativeReferenceMaker from GATK and reads are remapped to that.
HISAT2 offers the possibility to provide SNP calls to the index builder, and build an index that accounts for these SNPs.

Does anyone have any advice on which of these might be better?

As you said, this is a well-known problem for ASE. In Degner et al. 2009, they state that, "Perhaps surprisingly, we found that masking known SNPs does little to eliminate inherent biases in read-mapping."

As such, substituting N bases is not likely to completely solve your issue. Depending on the aligner, these N bases may also negatively affect your ability to map reads.

Of your two suggestions, I believe HISAT2 will be the better approach. There are also more intensive methods you can use such as empirical Bayesian methods or more traditional variant phasing (ex. AlleleSeq).


Why are TPMs per 10k or 100k in many scRNA-seq studies? I noticed that many scRNA-seq papers normalize TPMs to 10k or 100k as opposed to 1M (as the abbreviation defines them). It doesn't really matter since you are just moving the decimal point, so why mess with an established convention?

Additionally, the values might be log-transformed. In that case, log(x+1) would affect the data less if the total is 1M.

Collecting all the question comments as a community answer (feel free to edit this to make it more readable):
There tend to be far fewer counts in single cell experiments. By reducing the scaling factor you are avoiding having massive expression values for your samples.
It isn't uncommon to scale single cell experiments by the median cell library size. That way small values typically won't be too small and large values won't be too large. People also choose 10,000 because that's roughly around the library size for a single cell.
I haven't seen an explanation yet. Here is an example paper from last week (reputable journal and lab): doi.org/10.1016/j.celrep.2018.11.056 . "... We also excluded outlier cells from the MPPs group which were likely caused by cell contamination. To correct for batch issues, the two experimental batches were centered to have the same mean log(TP100K+1) (hereafter referred to as log(TPM), where TPM stands for ‘transcripts per million’) per gene. PCA was done in R with scaled..."


The first column contains Ensembl gene identifiers, and the suffix is a version number that can be used to track changes to the gene annotations over time. From the Ensembl Stable IDs documentation:

  Ensembl annotation uses a system of stable IDs that have prefixes based on the species scientific name plus the feature type, followed by a series of digits and a version e.g. ENSG00000139618.1. The version may be omitted.

Following the first link you provided leads to a page with details for the file 2edcaaa7-63b4-40b4-abbe-5d7a84012e60.FPKM-UQ.txt.gz. The first thing that caught my eye about this filename was FPKM, or "fragments per kilobase of exon per million reads", which is a commonly used unit of RNA expression. Since these are not raw read counts, there's no expectation that these values should be integers.

The best explanation I've seen of FPKM comes from a blog post written by Harold Pimentel of kallisto and sleuth fame. From the blog post:

  The interpretation of FPKM is as follows: if you were to sequence this pool of RNA again, you expect to see FPKM_i fragments for each thousand bases in the feature for every N/10^6 fragments you’ve sequenced. It’s basically just the rate of fragments per base multiplied by a big number (proportional to the number of fragments you sequenced) to make it more convenient.

More generally though, even when FPKM is not the unit of expression abundance used, most quantification methods and associated units of expression will not produce integer estimates.


How do programs like STRUCTURE determine populations from alignments? If I were to run BAM/SAM files through a program like STRUCTURE, or ADMIXTURE described  here how would they determine which individuals belong to what groups?

What part of the files would these programs use to determine which individuals belonged to different populations? I believe they do this via determining allele frequencies for each individual and then putting them into a population based on allele frequency, is the only criteria that they use comparison of allele frequencies between individuals?

structure works with SNP/genotype information, so won't work directly off BAM/SAM files. You would need to convert to VCF (e.g. using mpileup/vcftools), then process using something like PLINK to convert it into a format that STRUCTURE recognises.

Here's the description I have about structure from my PhD thesis. It might be a bit dated now, but I expect most of it is still applicable:

A clustering program called structure was created by Falush et al. (2003) that is used to determine population structure using a set of unlinked genotypes. The expected use case for structure is in the analysis of population structure, but it can also be useful for determining the utility of a given set of markers for categorising (or quantifying) a given trait.

The main benefit of using structure over other methods is that it is able to derive estimated ancestry coefficients from the data, rather than a discrete yes/no for each individual. This allows for the use of cutoff values for the assignment of individuals (useful for tweaking false positive and false negative rates in diagnostic tests), and a quantifiable value where a classification is not useful or does not make sense (as in the case of a continuous trait with high heritability such as height).

The program begins by assigning each individual to a given population. This can be either given in the input data as a ‘popinfo’ flag, or determined randomly by structure. After this, the program constructs a probability distribution for the data, considering the population assignments given, and uses that to estimate population allele frequencies. The program then constructs another probability distribution for these data, and uses that to estimate the population that each individual is likely to have originated from. This process is repeated many times (such a process is known as a random walk), and the underlying theory suggests that the progression of population assignments through the random walk process can be used as an approximation of the actual populations of origin for each individual.

The documentation included with structure describes the general algorithm for each iteration (of maybe hundreds of thousands of iterations) as follows:

Sample (population allele frequencies for iteration m) from the probability of (allele frequency) given (genotypes of sampled individuals) and (population of origin of individuals for iteration m − 1).
Sample (population of origin of individuals for iteration m) from the
probability of (population of origin) given (genotypes of sampled individuals) and (allele frequencies for iteration m).

In more simple terms, the steps are similar to the following statements:

these populations probably have this allele frequency distribution
therefore...
these individuals probably came from this population

The structure program is able to accommodate admixture into its model
by considering the population of origin as a probability, rather than a single
discrete value.

The file that is output from structure is a text file with individual identifiers and predicted ancestry coefficients (Q values). While the program does provide means to display these data graphically, the visual representation of data (including sorting) is better controlled using an external program. The most commonly used program is Noah Rosenberg’s distruct (Rosenberg, 2004), a program that allows colour customisation and group re-ordering. For the purpose of this thesis, a custom R script, snpchip2structure.R, has been designed which also allows selection of different sorting methods, a scatter plot for two populations, error bars, and a few additional features (see figure 1.7). More details about this script can be found in the programs directory of the supplementary CD for this thesis.


How does one convert a log-likelihood substitution matrix into a probability matrix? I am trying to calculate the substitution probabilities for amino acids (i.e. the probability of one amino acid given another). In theory, you should be able to derive them from a log-likelihood matrix (also called log-odds), such as Blosum62.
Does anyone know how to calculate probabilities from log-likelihoods or have a good modern source for a substitution probability matrix?

UPDATE2
@chubs wants to compare different clades. This is possible within recent releases of ConSurf:

It is additionally quite possible to use ConSurf's default outputs to get per-position conservation scores that can be used for post-processing in any fashion, as well as prefab scripts that can be used with PyMol to make custom visualization:

For all cases, ConSurf creates the following outputs:

The sequence and MSA colored by ConSurf conservation scores.

A text file that summarizes for each position the normalized score calculated, the assigned grade (1-through-9), the reliability estimation (for the Bayesian method) and the amino acids/nucleotides observed in the respective MSA column.

The sequences selected for the MSA and the MSA constructed (unless those files were uploaded by the user).

A file with the frequency of each amino acid/nucleotide observed in each column of the MSA.

The evolutionary tree, which was calculated by the server or uploaded by the user, together with the MSA are shown using the WASABI platform (Veidenberg et al., 2015).  Moreover, using WASABI, a user can select a subtree containing a fraction of the homologous sequences and conduct a follow-up ConSurf analysis with these selected sequences. To refine a ConSurf analysis to a selected subtree, a user can choose any internal node on the phylogenetic tree and open a WASABI menu by using a right mouse click. Selecting the option 'run ConSurf on subtree' will issue a new window with the new ConSurf run for the selected subtree sequences (see example in Figure 3).

UPDATE
It seems that what @chubs really wants to do is map protein conservation onto a structure. There are various tools for doing this (here, here). ConSurf is the older and more established tool, using evolutionary conservation to rule out trivial apparent associations due to phylogenetic non-independence.
here is the  somewhat related 3DPatch workflow, which uses profile HMM information as a measure of conservation:

Here is the output that they get:

Likewise ConSurf gets similar looking outputs:

ORIGINAL ANSWER
I think that more information is needed here to answer this question properly (what is your data? what is your ultimate goal? what are the things you've already tried?), but we can say a few things right away.
I would strongly suggest looking into substitution models, which are literally substitution probability matrices. Note that these are probabilities per unit time, which is to say that they are calibrated to e.g. evolutionary rates along a specific phylogenetic tree. I would therefore strongly suggest incorporating a phylogenetic model of some kind.
If you have e.g. a multiple sequence alignment, you can estimate a phylogenetic tree by one of the many tree estimation methods, which will give you a substitution probability matrix "for free" as this estimation is taken care of in the process of tree estimation. These will again be substitution probabilities per unit time; you will honestly not be able to get rid of the time-dependence here, as that's how evolution works.
I am a little behind the times, but some tools you might look into for tree estimation are PhyML and FastTree. If you only have sequences and haven't aligned them yet, you can use a tool like MUSCLE, ClustalO, or MAFFT.
BLOSUM and related models are derived from specific slices of protein sequence identity (evolutionary history) calibrated to certain scopes, it's unlikely that they are appropriate to your case, unless your data fits inside those rather narrow slices.


Multiple sequence pairwise global alignment against reference sequence Background

I have a fasta file with protein sequences. The first sequence in the fasta file is the reference sequence and every subsequent sequence is something I want to align to the reference. The end goal is to identify every unique sequence and be able to compare the regions of variation.

What I have tried

I'm looking into different tools such as Clustal Omega and MAFFT, but there is no option to select the reference sequence against which all the other protein sequences will be aligned.

Questions

Q1: I want perform one pairwise global alignment of each sequence in the file against the first (reference) sequence. What software can I use to do this?

Q2: What specific commands can I use to achieve Q1? I'm looking for specific instructions on how to complete this, not just a suggestion of a tool, since I'm not able to find anything from a relatively comprehensive search of available tools.

Q1

From this phrase, "The end goal is to identify every unique sequence", I am not 100% sure of what you mean by "every unique sequence". It sounds like maybe you want to group alignments of homologs and visualize them together?

If the proteins are similar enough, you might also benefit from trying BLAST to look at similar regions.

Otherwise MAFFT, Clustal, or MUSCLE are all appropriate tools to generate the alignments that you need. However I also do not know of any software aside from Geneious that has the option to align many proteins to a reference. If you were working with DNA I'm sure the stock suggestion would be to use a read aligner like bwa mem or minimap2 and to use a bam file.

Since I don't know of any software to do exactly what you would like, I wrote a small python script to run MAFFT many times for each sequence in a fasta after the first sequence, the reference.

Q2 - MAFFT/Python

In any case, here is a bit of python code that performs a MAFFT alignment of all of the sequences in a fasta file against the first sequence in the fasta file, which you have designated as the reference.

Here is the fasta file from Wikipedia. I called it seqs.fasta:

>reference
MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVSLRDLFDRAVMVSHYIHDLSS
EMFNEFDKRYAQGKGFITMALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHL
VTEVRGMKGAPDAILSRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDED
ARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC*

>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
ADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDADGNGTID
FPEFLTMMARKMKDTDSEEEIREAFRVFDKDGNGYISAAELRHVMTNLGEKLTDEEVDEMIREA
DIDGDGQVNYEEFVQMMTAK*

>gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
EWIWGGFSVDKATLNRFFAFHFILPFTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
LLILILLLLLLALLSPDMLGDPDNHMPADPLNTPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
IENY

Here is a script, mafft_ref.py. Run it with python mafft_ref.py seqs.fasta:

from Bio import SeqIO
from Bio.Align.Applications import MafftCommandline
import tempfile

fasta = sys.argv[1]

reference = 0
count = 0
with tempfile.NamedTemporaryFile() as temp:
    for record in SeqIO.parse(fasta, "fasta"):
        if count == 0:
            reference = record
        else:
            SeqIO.write([reference, record], temp.name, "fasta")
            mafft_cline = MafftCommandline(input=temp.name)
            stdout,stderr=mafft_cline()
            print(stdout)
        count += 1

Here is the output of python mafft_ref.py seqs.fasta printed to the terminal. Use python mafft.py prots.pep > output.txt to save the output to output.txt.

>reference
------------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQ
VSLRDLFDRA--------------VMVSHYIHDLSS-----EMFNEFDKRYAQGKGFITM
ALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHLVTEVRGMKGAPDAIL
SRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDEDARYSAFYNLLH
CLRRDSSKIDTYLKLLNCRIIYNNNC
>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
ADQLTEEQIAEFKEAFSLFDKDGDGTITTKEL------------GTVMRSL---GQNPTE
AELQDMINEVDADGNGTIDFPEFLTMMARKMKDTDSEEEIREAFRVFDK---DGNGYISA
A---------------------------------------ELRHVMTN-----------L
GEKLTDEEVD-------EMIREADIDG-------------------DGQVNYEEFVQMM-
-----TAK------------------

>reference
------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVS----
--LRDLF--------------------DRAVM----VSHYIHDLSSEMFNEFDKRYAQGK
GFITMALNSCHTSSL-------PTPEDKEQAQQTHHE-------------VLMSLILGLL
RS--WNDPLYHLVTEVRGMKGAPDAILSRAIEIEEE--------------NKRLLEGMEM
IFGQVI-------------------P-------------------GAKETE-PYPVWSGL
PSLQTKDEDARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC
>gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGY----VLPWG----QMSFWGA
TVITNLFSAIPYIGTNLVEWIWGGFSVDKATLNRFFAFHFILP-----------------
-FTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLGLLILILLLLLLALL
SPDMLGDPDNHM----------PADPLNTPLHIKPEWYFLFAYAILRSVPNK--LGGVLA
LFLSIVILGLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQM
ASI----------LYFSIILAFLPIAGXIENY--------------

References:
Some credit of mafft_ref.py to a biostars answer found here.


Proof of Breakpoint Reversal Sorting Approximation Algorithm I'm kind of new to bioinformatics and trying to self-study. I'm reading a bioinformatics book: An Introduction to Bioinformatics Algorithms and ran into some problems about understanding the proof of the Breakpoint Reversal Sort algorithm in chapter 5.4.

Specifically, the book states a theorem:

  If a permutation $\pi$ contains a decreasing strip, then there is a reversal $\rho$ that decreases the number of breakpoints in $\pi$, that is, breakpoint ($\rho(\pi)$) < breakpoint($\pi$).

This theorem guarantees that the algorithm will terminate and lead to no breakpoint. The book later explains the proof for the theorem:

  Let $k$ be the smallest element of decreasing strip, so element $k-1$ would be the ending of an increasing strip. Therefore element $k$ and $k-1$ corresponds to a breakpoint. Then reversing the segment between $k$ and $k-1$ will decrease the number of breakpoints.

I'm quite confused about this proof and don't really understand what the authors mean by reversing the segment between $k$ and $k-1$. Can someone please explain the proof for me?

Here's the lecture slide about reversal sorting that based on the book http://csbio.unc.edu/mcmillan/Media/Lecture07Spring2015.pdf.

Denote the strips in your permutation as $\langle S_1, S_2, \cdots, S_n \rangle$. Suppose $k-1$ is in segment $S_i$, and $k$ is in segment $S_j$, for some $i, j \in [n], i \neq j$.

Recall by definition of $k$ that $k-1$ must be the maximum and therefore rightmost element in an increasing strip, and $k$ the minimum and therefore also rightmost in a decreasing strip. Consider the two cases:

$i < j$. The increasing strip comes before the decreasing. $k-1$ is on the 'right' end of $S_i$, but $k$ is on the 'right' end of $S_j$ also. In order for $k-1$ and $k$ to be adjacent, reverse the segments $S_{i+1}, \cdots, S_j$.
$j < i$. The decreasing strip comes before the increasing. $k$ is on the 'right' end of $S_j$, but $k-1$ is on the 'right' end of $S_i$ also. Again, we reverse $S_{i+1}, \cdots, S_j$.

So, "reversing the segment between $k$ and $k-1$" means to reverse all the strips in-between the two strips containing $k$ and $k-1$, excluding the 'leftmost' strip, but including the 'rightmost' strip.


How to do the reverse of bcf merge? How can one do the reverse of bcftools merge?
To get one sample back out of the merged one.
Without holding reference sites.

bcftools has the view command which can do this. Create a text file (samples) which has the names of the samples you want, one per line, and then:

bcftools view -S samples merged.vcf.gz > new.vcf

Or, if you're running a shell that supports it, you can use input redirection for a single sample name:

bcftools view -S <(echo "sample name") merged.vcf.gz > new.vcf


Distinguishing Splicings that Produce Functional Transcripts Amongst Alternates Shown This is a link to the UCSC Genome Browser showing a portion of Ch 17 containing TP53 here.

There are several representations of the gene of varying length and number of exons. I presume these are equivalent to alternate splicings.

Assuming this presumption is correct(?), how can I determine which result in RNAs that are experimentally verified?

There are two resources worth checking: GTEx portal and Ensembl.
Uniprot and NCBI will call "canonical" the longest sequence comprising valid exons. These are not necessarily the biologically valid ones.

According to Ensembl (TP53, there are a bunch that are experimentally validated and are all class 1. Note the Uniprot column being nearly homologeous: to match the transcripts one has to look at the length, in this case TP53-201 is 393aa long, which is class 1. As mentioned this is not always the case.

So which are more expressed in tissue X?
In GTEx, under the collapsed div labelled "Exon Expression" press the "isoform Expression" tab. The most expressed is ENST00000620739.4 (354aa long) join with ENST00000269305.8 (393aa).

In the Exon Expression you will find out what exons are the issue.

Now. to get the sequence of a transcript or covert a variant Ensembl drops the ball. There is no clear web tab that I know of. Nor an API route.
ftp://ftp.ensembl.org/pub/release-99/fasta/homo_sapiens/pep/Homo_sapiens.GRCh38.pep.all.fa.gz contains all the sequences however. For a project of mine (ENSTMapper in this link) to map an amino acid mutation on one transcript to another I ended up write my own short converter, although there are a few similar tools out there.


How do you normalise read coverage in a BAM file? This is a question from the Oxford Nanopore community, from user Michael Radzieta:

  I have sequenced some plasmids using the rapid barcoding kit, I have attempted to assemble the data using several programs (unicycler, spades, ra) but nothing has been able to assemble the plasmids I am expecting. I have read that if the coverage is too high, then it can negatively affect the assembly. Based on the approximate size of my plasmids and the number of filtered reads I have, the coverage is estimated at >5000.

  Can anyone recommend a tool to subsample or normalise nanopore reads to a certain depth, such as 100x coverage?

I just updated one of my SAM-normalisation tools to more accurately normalise mapped reads to a target coverage. It uses reservoir sampling on a SAM-formatted file to choose reads to bring the coverage up to 100X (if possible).

https://gitlab.com/gringer/bioinfscripts/blob/master/samNormalise.pl

Example Usage:

samtools view -h input.bam | ~/scripts/samNormalise.pl | samtools view -b > normalised.bam

The output BAM file can be converted into a fastq file using 'samtools fastq':

samtools fastq normalised.bam > normalised.fastq

Here are some example before and after shots ['samtools mpileup input.bam'; 'samtools mpileup normalised.bam']:

Notice that where the coverage is below 100, the coverage is identical. Coverage is below 100 in some areas because some reads end, and no additional reads start until later on in the reference sequence.

The options are coded as variables at the beginning of the script at the moment, but I can easily convert them to command-line options in the future if there's demand for it.


Normalize RNA seq data from multiple runs for expression analysis I have RNA samples sequenced with TruSeq Stranded Total RNA kit protocol in Illumina HiSeq (2x125bp) and NovaSeq platforms (2x150bp) - almost 100 samples altogether. I have to use the samples data for expression analysis. Question is how do we normalize samples from different runs and of read length? What is best way (with available packages) we can normalize the count (to overcome the bias) and do a DE analysis? ANy suggestion would be helpful.

There are a number of methods. If you're doing DE, you have ComBatSeq, SVAseq, RUVseq, BUSseq. You could also try Z-score normalization ($(x - \overline{x}) / \sigma$) or even quantile transformation. For the latter two, make sure you work on each batch individually, not the whole dataset at once. See more on that here.
To visually compare if transformation works, plotting PCA/UMAP/t-SNE on raw and transformed data can perhaps be of some insight.


Approximate matches with bedtools intersect I am using bedtools to test whether a set of variant calls produced by a new method is missing any variants from a set of known validated variants.

bedtools intersect -v -a newcalls.vcf -b goldstandard.vcf

It's working as intended for SNV calls, but it occasionally misses indels due to a single base pair discrepancy. It's usually an alignment artifact that doesn't affect the interpretation of the variant.

My first thought was to make a BED file with intervals that extend a few bp beyond each variant call, and comparing this BED file to the known variants. However, I was wondering whether approximate matching is possible with bedtools intersect?

You can use bedtools intersect but first use bcftools norm.  bcftools norm will left-align and normalize indels to get both vcf files consistent with each other.  This should resolve the matching if they are truly matching variants.


Find Patterns in Cluster I have a heatmap and I would like to find some rectangles.

I have already used clustermap. But here, I can not calculate these rectangles.
The order of the data should not be changed.

This Code is not working:

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
tips2 = pd.read_csv('heatmapExport.csv')
sns.set(color_codes=True)
g = sns.clustermap(tips2, cbar="false")
plt.show()

Does somebody have an idea?

Thanks a lot

Here is an example of what I want:

Based on your description I think you should have a look at a technique called 'biclustering'.

The example on this page defines the goal of this technique as 'Finding subgroups of rows and columns which are as similar as possible to each other and as different as possible to the rest.'

Since your examples are python-based, you could check out scikit-learn's implementations of biclustering.


Aligning ChIP-Seq reads to repeats for downstream peak analysis This is a brief question regarding the above. I have previously used bowtie to map reads from paired-end ChIP-Seq sequencing, then used the positions for peak-calling. However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).

For peak-calling purposes, however, I need to know reasonably confidently where precisely those repeats are located, so that MACS2 treats the coverage pileup sensibly. So my question is this: what's the best way of mapping these (short) paired-end reads from a ChIP-Seq experiment so I can do this?

EDIT: I want to be able to run a peak-caller after aligning the data, so I want a good approximation to the coverage at repeating elements after mapping for instance - I don't want to throw away all reads which multimap, as may mean that some reads could lose their mates, and so I won't be able to calculate coverage or call peaks using the data.

I want an aligner to keep both ends of the paired-end pair within about 180bp of each other (as that is the expected length of the DNA fragment that is sequenced). The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced. With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.

I have tried bowtie (and bowtie2) for this, which is where I understood that the dataset maps to repeats.

This is very good question, but as far as I know there really isn't a good answer. I will attempt to provide some comments and tips.

  However, I'm trying to do similar with a dataset which may consist of repeats (I have previously aligned the dataset using bowtie, and a lot of reads aligned to repeats).

If this dataset is also PE and the alignment is outputting either uniquely mapped reads or a single (random) position for multimappers, it is pretty much as good as it gets. A lot of repeat regions are not exact copies, so using PE there is a good change that one read of the pair mapping to a single position is helping with mapping of the other read.

  I need to know reasonably confidently where precisely those repeats are located

It is probably a matter of phrasing, but what one needs is to know the precise genomic location of the reads. The location of the repeat elements can be obtained from the Repeat Masker track of UCSC, but you probably know this already.

  With multimapped reads this step could be done by seeing which of the pair is uniquely mapped and choosing the mate which maps closest to it.

I think bowtie already does this. Did you map the read pairs together or separately?

Note that:

There will always be a number of multimapping reads (fragments) that will be discarded. The information of how many are being discarded for your samples would be helpful.
Mapping can be improved using both reads and tunning bowtie settings (see bellow)

  Finding a valid paired-end alignment where both mates align to
  repetitive regions of the reference can be very time-consuming. By
  default, Bowtie avoids much of this cost by imposing a limit on the
  number of "tries" it makes to match an alignment for one mate with a
  nearby alignment for the other. The default limit is 100. This causes
  bowtie to miss some valid paired-end alignments where both mates lie
  in repetitive regions, but the user may use the --pairtries or -y
  options to increase Bowtie's sensitivity as desired. source

  The immediate step after aligning the pairs is to flatten the pair into the DNA fragment that has been sequenced.

As indicated above, the aligner will take care of that sort of. It will gave the mapping information of each read pair, which is then used by MACS.

Other tools

I tested ShortStack which uses the coverage of uniquely mapped reads to place multimapping reads for ChIP-seq, and the results were mixed. That said, it was SE data and the code has been through some changes, so it be worth a try.

When requested, the authors of RepEnrich provided me with the scripts for an utility called RepConsensus which maps reads to a consensus reference of repeat regions. The caveats are:

mappings are better for elements whose sequences have not diverged a lot yet
it will tell the position of the protein in the repeat sequence, but no information about the surrounding genomic regions because there isn't any.

Whatever you do, do not use all alignments, that is double counting of reads, and think very carefully before using custom methods. A blog post detailing some of the potential issues can be found here (I wrote it).


Converting VCF file to PLINK bed/bim/fam files I am trying to find the best way to convert VCF files to PLINK binary bed/bim/fam files, but it seems like there are many varied ways to do this. (For example, using Plink 1.9 --vcf tag, bcftools, GATK, and vcftools).

Obviously they all probably have their uses and you could likely spend a while discussing this, but can anyone give a short, simple explanation as to whether the methods are generally the same, or whether one or another is generally best for certain situations? Thanks!

Given that Plink reads in VCF files natively, and Plink is your desired output format, I expect that plink --vcf <file> would be the best option.

There are caveats associated with the types of variants that Plink will pull in from the VCF files, mostly related to limitations of Plink's own file formats. But if your end-goal is Plink, then that should be fine.


How to subset a GRanges via an argument passed into a function? Let's say I have following example GRanges:

> library(GenomicRanges)
> gr = GRanges(
+     seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
+     ranges = IRanges(101:110, end = 111:120, names = head(letters, 10)),
+     strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
+     score = 1:10,
+     GC = seq(1, 0, length=10))
> gr
GRanges object with 10 ranges and 2 metadata columns:
    seqnames     ranges strand |     score                GC
       <Rle>  <IRanges>  <Rle> | <integer>         <numeric>
  a     chr1 [101, 111]      - |         1                 1
  b     chr2 [102, 112]      + |         2 0.888888888888889
  c     chr2 [103, 113]      + |         3 0.777777777777778
  d     chr2 [104, 114]      * |         4 0.666666666666667
  e     chr1 [105, 115]      * |         5 0.555555555555556
  f     chr1 [106, 116]      + |         6 0.444444444444444
  g     chr3 [107, 117]      + |         7 0.333333333333333
  h     chr3 [108, 118]      + |         8 0.222222222222222
  i     chr3 [109, 119]      - |         9 0.111111111111111
  j     chr3 [110, 120]      - |        10                 0
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

It is straightforward to subset this GRanges as one would subset an R data.frame, e.g.

> gr$score
 [1]  1  2  3  4  5  6  7  8  9 10

However, I am writing a function where the user would pass through a string for the metadata column to subset e.g.

subsetter = function(granges, column){
    return(granges$column)
}

I try to use the function as follows:

> subsetter(gr, score)

which outputs NULL

This obvious doesn't work, as it will subset the gr by column, not score.

> subsetter(gr, score)
NULL

If this was a data.frame, I could subset as follows:

gr[[(score)]]

but this gives the error

Error in gr[[(score)]] : this S4 class is not subsettable

Are there any other ways to subset a GRanges besides using $?

subsetter = function(gr, cname) {
    return(mcols(gr)[[cname]])
}

You can then use things like subsetter(gr, "GC") and subsetter(gr, "score").


How to construct the phylogenetic tree of the "High Recombination" genus streptomyces? I need to construct the phylogenetic tree of the whole genus of streptomyces. And my goal is to find out how my target strep species distributed in the whole strep genus and obtain the distribution pattern. Strep is a big genus, has 658 species(65% of them have the genome info), and my target species is about 25. But the traditional 16S method seems not helping when it comes to strep is a genetic recombination genus. And still lack of MLST scheme to distinguish my target species. I want to get the phylogeny for the whole genus, but I lack related genetics knowledge.

There is something I did before try to get the tree.

I ran the orthomcl program for 400 strep genomes which I got from NCBI, and the result is: I got 22 orthologues
I ran the fastani program for 200 strep genomes(type material) and got the ani result.

What can I do to achieve my goal? Thanks!

You need bootstrap support using a model-based tree building algorhithm, via maximum likelihood (a few people use Bayes). The file format is relaxed phylip format (please submit a separate question if you have difficulties here - its a bit tricky).

I use standard RAxML here, specifically raxmlHPC (easily downloadable and compiles on Linux and OSX). The codes are quite complicated and I've given them below.

A robust maximum likelihood tree is,

raxmlHPC -m GTRGAMMA -n output.tre -o outgroup-p 10000 -s alignment_file.phy

To make a bootstrap tree,

raxmlHPC -b 459 -m GTRGAMMA -n output.tre -o outgroup -p 10000 -s /path_to_alignment.phy -N 500

This tree will bootstrap for 500 replications, however to start I would use 100 replications.

Make a consensus tree of the bootstraps,

raxmlHPC -f b -m GTRGAMMA -n output.tre -o outgroup -p 10000 -t /path_to_master_tree(best_tree) -z path_to_bootfile

You require bootstrap support >80% and please repeat this with and without 5-2/5-3 (it still looks long)

The cluster you have access to is fine for the calculation, it will take aound 24 hours to complete one bootstrap calculation for one data set and obviously you need to parallelise you calculation across the 22 contigs.

Viewing the tree, FigTree (for Mac OSX) is easy.

Root

Rooting can be complicated because I don't really know your bacteria.

The recombination issue is more complicated, but I would construct 22 trees from your contigs and assess them for congruence. Panmixia is a concerned, which means too much recombination


What type of information does an .sf file contain? What is an .sf file.  What information does it contain?

You should provide more context, otherwise one can only guess what is it you are looking for! In case you are asking about Salmon's output, below is from their documentation:

  Salmon’s main output is its quantification file. This file is a
  plain-text, tab-separated file with a single header line (which names
  all of the columns). This file is named quant.sf and appears at the
  top-level of Salmon’s output directory. The columns appear in the
  following order:

Name   Length   EffectiveLength   TPM   NumReads


mapping and de novo assembly of plant without reference genome I'm studying on a plant that has no reference genome and only has one scaffold assembly and one gff3 annotation file. Can I create an index with the same assembly and gff3 in STAR and do the mapping? If I want to do de novo assembly using the Galaxy Can I import the links of 20 fastq files from the RNA-seq of the samples together to construct the Reference Genome with Trinity? If this is possible, please help me import the links of these files together into Galaxy and make the reference genome.

Thanks a lot

You have multiple questions, I will try to answer one by one:

On creating an index with STAR: As long as you provide a fasta file and its corresponding gff, you will be able to generate an index.
On creating a reference with Trinity: You have mentioned "to construct the Reference Genome with Trinity", however, with your RNA-seq data and Trinity, you can construct a de novo transcriptome and not a genome.
Uploading files to Galaxy: I suggest to do so via FTP using FileZilla or something similar. The transfer can take some time if your files are large though.

If I understand your questions correctly, you are performing RNA-seq and you are concerned about the completeness of the genome of your model organism. In that case, I would advise to bring together a de novo transcriptome from your RNA-seq data (with Trinity), use one of the quasi-mapping methods (with a tool like Salmon or kallisto).


how to merge two fasta based on full fasta header and remove the duplicate I am wondering if anyone can help to merge two fasta file based on the full header as initially gene ID may be similar in some of the sequences and to remove the duplicate entry based on full fasta header

I do have two multifasta file:

file1:

>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

file2:

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

Expected output:

>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

You can use the FastaToTbl and TblToFasta scripts I have posted here for this:

FastaToTbl file1.fa file2.fa | sort | uniq | TblToFasta > merged.fa

If you run this on your example, you will get:

$ FastaToTbl file1.fa file2.fa | sort | uniq | TblToFasta > merged.fa
$ cat merged.fa
>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATG
GTGCGGGAACGCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGT
TTGGATGTCTCTGGCGACGG
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGC
GAACGCTACGGGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGAT
GTCACCGACGACGGTGGGTG
>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAG
GTCGCGGTCGGCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTC
GACCGCCTCGAAGCCCGCCT
>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAG
GTCGCGGTCGGCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTC
GACCGCCTCGAAGCCAGCGC
>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGC
CAGTTGGTACGGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGAC
GTCGGCCTCGATGTCACGGG

This will remove any cases where both the header and the sequence are identical. If you instead want to do this based on the header only, not taking the sequences into account, you can just process both files and only print sequences from the 2nd if their ID line hasn't been seen before. For example, with awk:

$ awk '{
        if(NR==FNR){
            if(/^>/){seen[$0]++;}
            print
        }
        else{
            if(/^>/){ skip=seen[$0] }
            if(!skip){ print }
        }
       }' file1.fa file2.fa
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC
>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT


Replace lowercase characters with - I have an output from vcfutils.pl vcf2fq with specified minimal depth, and it means that nucleotides with not enough depth are lowercase.
I would like to change them to gaps. I could do it in higher language but I would like to know the solution in bash.
I have found simple solution to replace all lowercase to uppercase:

awk 'BEGIN{FS=" "}{if(!/>/){print toupper($0)}else{print $1}}' in.fna > out.fna

but I would like to replace lowercase with -. The first row need to be omitted. Could you help me? Thanks

The following sed command will do the trick.

sed -e '/^[^>]/ s/[a-z]/\-/g' in.fna > out.fna

Lines starting with > will be escaped and the lowercase characters will be replaced with - .


Biplot scores from canonical correspondence analysis I'm using the R package vegan to perform canonical correspondence analysis (CCA). As input we have two matrices, one being (sites)x(species) and the other being (sites)x(conditions).

Sample data (and source of plot) are here.

Species loadings are easily accessed with summary(cca_model)$species. What I'm trying to find is the loadings for the explanatory variables, the conditions. The only summary I can find is biplot scores. Looking through the documentation for vegan I can't find any description of how they are computed. Can I sum them across CCA components to get an idea of how much they influence the data?

This is a biplot of two CCA components. The scores are used as coordinates for the arrows.

So two questions:

  What are biplot scores in the context of CCA?

and

  Can biplot scores be used to determine how much of an effect conditions have on the response variables?

Canonical correlation analysis (CCA) not correspondence canonical analysis. Biplots are the same as in PCA: the representation in a new space of the samples and the relative position of the variables on the new space.

However, CCA seeks the values of one block that maximize the correlation with the other block. Here you don't have an independent variable that you maximize from the variables you introduced on the CCA. So the interpretation is much trickier as there isn't a direction: WaterC and FallenTwigs show the opposite effect than BarSand and CoverMoss, but there is a group of variables that explain what is the effect of water on the plants(?)


Multiple Alignment cost application Is there any biological application/case where someone would be interest in the estimated total cost of the alignment between a set of sequences (genes or amino acids) without the aligned sequence itself?
Edit: to clarify by cost I mean the cost function such as the edit distance. I have some ideas for heuristics on the multiple sequences alignment problem (which is intractable) but these would not retrieve the aligned sequence, therefore, I was wondering whether such algorithm would be useful

I'm not a fan personally, but it is context specific. The problem is that it is heading towards optimisation, particularly algorithmic optimisation. If the alignment is the final calculation - that is cool, but the situations are select. It depends how far you want to take it.
The rationale of my skepticism is simple - if an alignment was decided on a 'maximization' criteria, be it likelihood or some other metric - the subsequent calculations based on the alignment, such as phylogenetics Jukes Cantor, Kimura 2 and 3-parameter etc ... we enter circular mathematics. The critiera that defined the alignment is the criteria we are using to perform unrelated post hoc analysis. In truth there is a degree of circular mathematics in any algorithm used to produce an alignment, but we overcome this by optical verification precisely avoid this.
In this sense the alignment is simply to establish homologous amino acid/nucleotide positions for downstream analysis. Thus optical optimisation, whilst cautioned against in every other area of bioinformatics, is considered useful in alignments because of the risk that a logical error in the alignment is amplified in downstream analysis, potentially leading to an erroneous conclusion.
Algorithms are a problem where indels are concerned and aligning across an indel. It is very easy to make errors here and 'optical analysis' would simple introduce a blockwise deletion to avoid the uncertainty, i.e. if you are not sure just chop it out - better to lose data than risk artefacts on an equivocal judgement call.
I'm not saying don't try, I am saying don't look to overstate the applicability of the metric.


What are negative FSC-A values in flow cytometry data? I did flow cytometry analysis of E. coli cells 7 hours post inoculation (1% overnight culture in 3ml).
I am analysing the FCS3.0 files using the flowCore package in R. I observe that FSC-A has negative values in some events, although the corresponding FSC-H values were all positive. These points comprised just a minor fraction of the data (~250/10000 in this case). The FSC voltage was set at 788V and the threshold flow rate was 5000-10000 events/sec.
What do these negative values mean?

FSC-A vs FSC-H: plotted only negative FSC-A values

SSC-H vs FSC-A: red points indicate negative FSC-A values

This is probably a matter of calibration. If this is not the raw data (actual voltage curves), look at that.

As an area below 0 is not possible, it has to be a calibration issue, either in the machine, or in the software.

the -H plots are about the maximum height of the voltage spike in the side or forward light detector, while the -W plots are about the width (using some definition of flank-detection, so variation prone).
The area is also computed in various ways: using something like (H*W)/c (with c something like 2) or with an actual computation of the area of the spike, but again this has all kinds of predefinitions in it (what's the beginning/end of a  spike, what's the background-level, ...) Some of these predefinitions you can alter in software, some are hardcoded or even part of the functionality of the hardware, so it might be easy to impossible to get rid of the negative values.

The FSC-H values corresponding to the FSC-A values are in your plot, and positive, so the most obvious solution (a bad algorithm that calculated Area as Height times Width times Magic and did not account for negative Heights (which depending on spike definition, you might get)) is out.


Unable to read all alignments in stockholm format : UnicodeDecodeError: 'utf-8' I have a file in Stockholm format and I am trying to read it using the Biopython package. The code is as follows :

import numpy as np
from Bio import AlignIO

alignments = AlignIO.parse("Rfam.seed", "stockholm")
for alignment in alignments:
    print(alignment)

It gets read for some alignments and then throws up following error :

Traceback (most recent call last):
  File "test.py", line 5, in <module>
    for alignment in alignments:
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Bio/AlignIO/__init__.py", line 394, in parse
    for a in i:
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Bio/AlignIO/StockholmIO.py", line 408, in __next__
    line = handle.readline()
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/codecs.py", line 322, in decode
    (result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xe1 in position 2468: invalid continuation byte

How to resolve this error?

The last line indicates that this is an encoding problem. Try to see if the function you are using to read the data has an "encoding" parameter that you can set to whatever encoding your file is using. If that doesn't work, you can change the encoding of your file by:

with save as using a text editor of your choice and selecting "UTF8" as the encoding. That would be a quick work-around if you have a small number of files.

or

You can use the codecs library to check and convert file encodings within your code. Here is an example.

Warning: I haven't used Biopython in ages but encountered similar problems with Python in the past.


Recommendations for missing value imputation - DNA methylation data I'm looking for some options for imputation for a high-dimensional dataset of DNA methylation (bisulfite sequencing) data.  Dimensions on the order of 50-100 samples x ~500,000 CpG loci/features.

I've used K-nearest neighbors, but It seems that this method is not very accurate. As far as I can tell, it limits the minimum number of "genes/features" to impute at one time to something like 1500. which is kind of a lot.

I've also used missForest to impute smaller datasets with greater accuracy, but it seems computationally infeasible to do it on the full 50x500,000 dataset.

I need to impute values because I'm using the data for some statistical modeling which require complete cases.

Anyone know of better alternatives than K-nearest neighbor for large scale imputation?

I think this is still an active field of research.  I have heard of Phenix, which might be appropriate.