Proper use of BWA MEM on multiplexed GBS sample I have a multiplexed lane of GBS sequencing reads as a fastq file. I understand the first step is to demultiplex and trim the adapter sequences from the reads.
This yields many individual fastq files that correspond to individual libraries.
I would like to genotype my samples, but I need to map the reads first. My question is how should I perform the BWA MEM alignment?
Do I perform BWA MEM on each of the resulting fastq files that correspond to individual libraries...
OR
I should I perform BWA MEM on the large original fastq file that contains all my sequence reads before demux and read trimming?
Forgive my ignorance, I've only ever used TASSEL (which works great) and I'm trying to expand my understanding and flexibility.

The first thing to do is to group the FastQ files together following this post. Once you have your reads grouped, you could trim your files before mapping. Then you could perform your BWA alignment. Since you mentioned MEM, I suppose your read length should be between 70-100bp as suggested by the author:
It consists of three algorithms: BWA-backtrack, BWA-SW and BWA-MEM.
The first algorithm is designed for Illumina sequence reads up to 100bp, while the rest two for
longer sequences ranged from 70bp to 1Mbp. BWA-MEM and BWA-SW share similar features such
as long-read support and split alignment, but BWA-MEM, which is the latest,
is generally recommended for high-quality queries as it is faster and more accurate.
BWA-MEM also has better performance than BWA-backtrack for 70-100bp Illumina reads.

Then you could run your BWA MEM, supposing you have 10 CPUs available, like so:
bwa mem \                                  # bwa mem algorithm
-M \                                       # mark shorter split hits as secondary (for Picard compatibility).
-t 10 \                                    # CPUs thread
reference.fasta \                          # reference genome to be aligned to
< sample_r1.fastq \                        # fastQ read 1 (can be .gz)
< sample_r2.fastq | \                      # optional FastQ read 2 (can be .gz)
samtools view -bS -@ 10 > sample.bam       # pipe to samtools to convert SAM to BAM on the fly.

Hope this helps !


Hardware Requirements (specs) for Bioinformatics-dedicated desktop I know this is a somewhat general or vague question, but I’m interested in your opinions. I must build a desktop for general bioinformatics activities with human genomes. I will work with Python and R libraries, I will also need to install denovo sequence assembly software, alignment tools, variant viewers/browsers, sql, among other software. I don’t have much experience customizing hardware that will be dedicated to bioinformatics and genome data analyses. I’m pretty sure I work with linux and/or iOS. What is the recommended RAM, processor speed, HD capacity, graphics cards, etc. I appreciate your opinions on the matter. Thank you.

The RAM requirements for some bioinformatics analyses like assembly can be quite high (in the hundreds of Gigabytes). My recommendation is to get a fast laptop. Something with an i7 quad-core processor, 16 GB of RAM, and 1TB of storage should do. Do your analyses on this laptop. If you run out of RAM or storage, then use a cloud service to get access to bigger machines or buy an external hard drive. Once you have a feeling for what your actual needs are, start looking at a dedicated rack server solution in your local institution's server farm or again, in the cloud.

You don't mention any analyses that clearly benefit from a GPU, but as you get more familiar with your analyses, perhaps that will become useful as well.


Why do BAM files created by different tools have different file sizes? I have a BAM created by Picard. I want to filter alignments by flags with samtools view. However, I noticed that even if I apply no filters, the output BAM is different from my input BAM. Are BAMs produced by different tools also different in size? How can I check if two BAMs are the same?

It's worth bearing in mind that when outputting compressed BAM, as most tools do by default, they may well be using different levels of compression and/or different libraries, or versions of said, libraries for doing (de)compression which will result in different file sizes.  Additionally coordinate sorted BAM will compress more than unsorted BAM. The current version of Picard uses HTSJDK which in turn uses java.util.zip.Deflater/Inflater, current versions of samtools should be using the HTSlib library which in turn depends on the standard zlib library. You can see the effect of different implementations of zlib have on file size and execution time in benchmarking done by the samtools team.

However, in your case the best way to see if there is any difference between the BAM files is to rule out the effect of different levels of compression or libraries used for compression and save both BAM files as uncompressed.  Both samtools and Picard have options for disabling or changing the levels of compression, since the BAM compression standard BGZF is implemented on top of the gzip format it has inherited the ability, just like with gzip, to change the compression level from 0 to 9.

samtools view -bu will allow you to produce uncompressed BAM output (which is also handy for piping into other programs as it saves time wasted compressing decompressing what is essentially a stream).  Also note that samtools sort has a -l INT setting where INT can be set between 0 (compression off, as with -u) 1 (for fastest compression, but increased file size) or -9 (for maximal compression, with increased run time).  Some of the effects of increased runtime for higher compression settings might be ameliorated using the -@ argument which allows you to set the number of extra threads used for BAM compression, by default samtools won't use any.

Picard tools has a general setting COMPRESSION_LEVEL which is applicable to most of its tools setting this to 0, COMPRESSION_LEVEL=0 should disable compression.

So re-running whatever Picard tool you used in the first instance with COMPRESSION_LEVEL=0 will then enable you to check that the file is not further modified by samtools view -bu. The assumption here being that if both files have exactly the same content they should be the same size uncompressed, of course if they have trivial differences with regard to white space formatting things may still differ.


Calculating Z-score from logCPM values using edgeR I have the raw counts for RNA-Seq data. I converted counts data to logCPM using edgeR package.

Lets say I have a dataframe A with 15000 genes as rows and 100 samples as columns with counts data. I have used following code from edgeR tutorial.

logCPM <- cpm(A, prior.count=2, log=TRUE)
head(logCPM)[1:5,1:5]
               Sample1          Sample2          Sample3         Sample4
A1BG          -1.019429        0.8333739        -1.192960       -0.6717287
A2M            7.622280        7.8057767         7.932570        8.6409077
A2ML1          1.004541        1.0909119         5.462613        6.4149348
A4GALT         2.126502        1.8403954         5.322601        2.9469345
AAAS           5.785888        5.1639501         5.259789        5.3248026
               Sample5
A1BG          0.3121975
A2M           9.1949911
A2ML1         9.4740594
A4GALT        2.9214979
AAAS          4.3708063

After converting counts to logCPM I see there is another command for scaling which is used for heatmaps.

logCPM <- t(scale(t(logCPM)))

It is given there that scaling each row (each gene) to have mean zero and standard deviation one in page 15 edgeR

Is this nothing but z-score? or Is there any otherway to calculate z-score from logCPM data?

thanq

Yes this is nothing but a z-score. The scale function calculates column z-scores, but since you use t (transpose), it is actually performed on rows. There are many other ways to calculate z-scores, but if done correctly they will all give the same results.


Relaxed and sequential Phylip format conversion I routinely use AlignIO in BioPython for manipulating alignments and in context here moving from fasta format to phylip format. Alot of phylogeny packages accept a relaxed (sequence ID) & sequential (sequence) phylip format. This format is user-friendly, helps double-check for alignment errors and can be the essential input format.
A typical bit of code is,
'''
Converts a fasta file to phylip format
'''
from Bio import AlignIO

Class Convert (inheritClass):

    def fa2phy (self):
        fileOut, fileIn = self._inOut() # builds the path name and deletes any pre-exising outfile
        align = AlignIO.read(fileIn, "fasta") # reads the alignment in fasta
        AlignIO.write(align, fileOut, "phylip-relaxed") # writes in phylip

The output of "phylip-relaxed" is interleaved format
4 120
virusA_detailed_description ACGTGACG .... ATCCCCAA
virusB_detailed_description ACGTGACG .... ATCGGGAA
virusC_detailed_description ACGTGACG .... ATCAAAAA
virusD_detailed_description ACGTGACG .... ATTTTTAA

ACGTGACG .... ATCCCCAA
ACGTGACG .... ATCGGGAA
ACGTGACG .... ATCAAAAA
ACGTGACG .... ATTTTTAA

Output sought is relaxed-sequential
4 120
virusA_detailed_description ACGTGACG .... ATCCCCAAACGTGACG .... ATCCCCAA
virusB_detailed_description ACGTGACG .... ATCGGGAAACGTGACG .... ATCGGGAA
virusC_detailed_description ACGTGACG .... ATCAAAAAACGTGACG .... ATCAAAAA
virusD_detailed_description ACGTGACG .... ATTTTTAAACGTGACG .... ATTTTTAA

There is an option to allow a sequential format (above). HOWEVER, a "phylip-sequential" option output is as follows (it truncates the sequence ID to 10 characters,
4 120
virusA_det ACGTGACG .... ATCCCCAAACGTGACG .... ATCCCCAA
virusB_det ACGTGACG .... ATCGGGAAACGTGACG .... ATCGGGAA
virusC_det ACGTGACG .... ATCAAAAAACGTGACG .... ATCAAAAA
virusD_det ACGTGACG .... ATTTTTAAACGTGACG .... ATTTTTAA

No matter how I look through AlignIO options the above is the nearest format.  Could anyone point to a IO style solution?
I could write it, but that defeats the object (i.e. purpose not an instance of a method).

Note the "...." in the outputs above are not a typical "sequence identity" symbols, it just means that total length of a line in an interleaved format is 60 nucleotides (it might be 80).

I have nutted out a solution for you.
If this is the content of input.fa:
>virusA_detailed_description
ACGTGACGATCCCCAAACGTGACGATCCCCAA
>virusB_detailed_description
ACGTGACGATCGGGAAACGTGACGATCGGGAA
>virusC_detailed_description
ACGTGACGATCAAAAAACGTGACGATCAAAAA
>virusD_detailed_description
ACGTGACGATTTTTAAACGTGACGATTTTTAA

And this is the script:
#!/usr/bin/env python3

'''
Converts a fasta file to phylip format
'''

from Bio import AlignIO
from Bio.AlignIO.PhylipIO import PhylipWriter
import sys

infile = sys.argv[1]
outfile = sys.argv[2]
idwidth = int(sys.argv[3])
with open(infile, 'r') as input_handle:
    align = AlignIO.read(input_handle, "fasta")
    with open(outfile, 'w') as output_handle:
        PhylipWriter(output_handle).write_alignment(align, id_width=idwidth)

Then, save this script as phylpwrtr.py and do:
chmod +x phylpwrtr.py
./phylpwrtr.py input.fa output_longids.phy 30 # 30 is the width of the header, adjust the value to adjust header length

The contents of output_longids.phy looks like:
 4 32
virusA_detailed_description    ACGTGACGAT CCCCAAACGT GACGATCCCC AA
virusB_detailed_description    ACGTGACGAT CGGGAAACGT GACGATCGGG AA
virusC_detailed_description    ACGTGACGAT CAAAAAACGT GACGATCAAA AA
virusD_detailed_description    ACGTGACGAT TTTTAAACGT GACGATTTTT AA


GTF from consensus sequence I am new to using bioinformatic tools and I was hoping this community could help clear something up for me. I need to generate a gtf file. My data are a set of complete HA genes for influenza B viruses in fasta format. From reading through forums and other internet searches, my understanding is that fasta cannot be transformed to gtf or gff because these last two provide location information which fasta does not (unless available in the header). I saw someone suggest mapping sequences to a reference sequence and generating the gtf from that mapping.

My question is, would it make sense to use the consensus sequence selected for alignment as this reference sequence? The alignment algorithm is a simple majority rule. Because all my data are complete genes, I would think the mapping of locations (start codon and stop codon) would be the same no matter which sequence in my data set I choose to use. Is this not so?

I am looking at viral evolution in the HA segment and my end goal is to calculate the dS/dN ratio for a set influenza B viruses collected during a particular season. The program I want to use requires a gtf, but likely other options would only require the reference genome. Would it make more sense to use the strain used as the vaccine candidate as the reference genome?

I appreciate any help you can provide me!

"Transforming" Fasta files to GTF or GFF3 files is a common request, as are similar tasks such as "converting" BAM files to VCF files. But these are underdefined questions, and "transform" or "convert" isn't really the right way to think about the task.

As far as your particular request is concerned, Fasta files and GTF files contain fundamentally different types data.

Fasta files contain sequences, which could be chromosomes, assembled contigs, polypeptides, transcripts, or other biological sequences.
GTF (and GFF3 and BED) files contain intervals, indicating that some feature of interest occurs at such-and-such position on such-and-such sequence. A hypothetical example: a transcription factor binding site occurs on chromosome 2 from position 40011 to 40017.

So you are correct: it's not really helpful to think of your task in terms of "transforming" Fasta to GTF. The process of creating a GTF/GFF3/BED file is often called annotation. There are many different methods to annotate many different types of genomic features: transposable elements, segmental duplications, protein-coding genes, non-coding RNAs, and so on.

As a result, it's meaningless to say "I need a GTF file" without clearly specifying what that GTF file will annotate. (It's kind of like saying "I need an XML file" or "I need an Excel spreadsheet" without indicating what the contents of those files would be.) The things you need to clarify are:

What will the intervals in my GTF file represent?
What precise sequences will serve as the reference and provide interval coordinates for the GTF file?
How are those intervals to be determined?

It sounds like 3. is a big part of your question. If I'm reading things correctly, you want to annotate a particular class of genes on an influenza genome. Is that correct? Did you perform de novo sequencing of influenza viruses, and now hope to characterize variation in their HA genes? Stating what you hope to get out of your analysis would be helpful in finding a good way forward.


Clustering RNAseq fold-change data I have a dataset that looks something like this

Treatment1
Treatment 2
Control

Sample1
3.23.
0.87.
1

Sample2
1.71.
1.
1

Sample3
2.88.
5.65
1

Sample4
0.77.
1.34
1

The numbers describe a fold change with respect to the control (where the control is always set to 1).
Does it make sense to compute values such as averages, sums, and differences (e.g. the average fold change for treatment1)? And yes, wary that doing things on log fold-change is not a thing, but how about just fold change?
By extension, does it also make sense to compute something like Euclidean distances (e.g. fold change vector of Sample1 vs. Sample2) for clustering fold change data? I think clustering on counts/TPM makes more sense, but wasn't sure about doing it on fold changes.
Something tells me that it's not possible because the fold changes are relative values as opposed to measured numbers. I've opted for looking at correlations (e.g. Spearman) for now though; are there alternative metrics/strategies to consider?
If there are some good papers talking about clustering and manipulating fold change data I'd be equally grateful, thanks!

Does it make sense to compute values such as averages, sums, and differences (e.g. the average fold change for treatment1)? And yes, wary that doing things on log fold-change is not a thing, but how about just fold change?

This is problematic, if you want to compare treatments then compute the actual fold change and its significance. You can't really base anything off of the fold-change (log or otherwise) alone.

does it also make sense to compute something like Euclidean distances for clustering fold change data

It's not unheard of to do this. The problem is you lose all information about how reliable those fold-change estimates are. This is something maintained when you cluster sample-level information (e.g., TPM or robust log transformed data). Not that using variance-stabilized data for clustering will work better, since it tamps down on how much highly-expressed genes influence your clustering (assuming you don't want that to occur).
Have a look at the DESeq2 vignette for some nice examples of what good practices are.


Looper: Duplicate sample name error I am using Looper. While running looper run project_config.yaml I am getting duplicate sample name error. I have two fastq files for each sample, so in the annotation.csv file I have:

sample_name,library,file
abeta-24-1,RNA-seq,data/merged/abeta_24h_1_r1.fq.gz
abeta-24-1,RNA-seq,data/merged/abeta_24h_1_r2.fq.gz
abeta-24-2,RNA-seq,data/merged/abeta_24h_2_r1.fq.gz
abeta-24-2,RNA-seq,data/merged/abeta_24h_2_r2.fq.gz

My pipeline_interface.yaml file looks the following way:

protocol_mapping:
  RNA-seq: sandro_rna_seq

pipelines:
  sandro_rna_seq:
    name: sandro_rna_seq
    path: sandro_rna_seq.py  # relative to this
pipeline_interface.yaml file
    looper_args: True
    arguments:
      "-L": "IU"
      "-T": "gencode_mouse_m13"
      "-E": "/scratch/nv4e/binf-
pipelines/bulk_rna_seq/salmon/pypiper/jyvo_experiment-metadata.yaml"
      "-C": "/scratch/nv4e/binf-
pipelines/bulk_rna_seq/salmon/pypiper/salmon_rna_seq.yaml"
      "-X": "True"
      "-P": "6"
      "-I": sample_name
    resources:
      default:
        file_size: "0"
        cores: "20"
        mem: "120000"

Why am I getting this error? How could I fix it?

As explained in the pep docs, the sample_name column should be considered a unique identifier for a sample. The PEP requires that One row corresponds to one sample You have multiple rows with the same name; this is not allowed.

You simply need to recast your annotation sheet such that it has file1 and file2 columns, but only 2 rows.


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

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

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

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


TagReadWithGene missing when using latest version of Drop-seq_tools I am using Drop-seq_tools to analyze scRNA-seq in a similar way of this paper. I am using the same data (GSE97930) to get used to Drop-seq_tools.

On the GEO repo the authors provide the relevant code (GSE97930_README.txt). However, one tool used in the script is missing: TagReadWithGene.

The latest version of Drop-seq_tools does not have such method (as it should be provided by this library), and there is no version history for the whole package, so I cannot download previous versions.

The best matches for the missing method are TagReadWithGeneExon and TagReadWithInterval.
In the methods of the cited paper they say:

  Reads that mapped to intronic or exonic regions of genes as per the GENCODE gene annotation were recorded.

I think that in this case the tool to use should be TagReadWithInterval since the other one should consider only exons.

I cannot find any other documentation about the tools. Do you know/already encountered this problem?

As per this response from Drop-seq_tools team, at the moment intronic regions cannot be annotated and TagReadWithGeneExon should be used instead.

However, in the version of the public repository:

  TagReadWithGeneExonFunction is the old code, for people who want the more simple tagging methodology and rely on it for 3rd party tools, and TagReadWithGeneFunction, which for your purposes replaces TagReadWithGeneExon. [...]  You’ll see additional options to extract expression for CODING, UTR, INTRONIC reads [...]


How to find the position-specific weight matrix (PSWM) for a list of genes having different sizes? Wiki has the algorithm for finding position-specific weight matrix (PSWM) for a list of nucleotide sequences having an equal number of nucleotide. How do we find the PSSM for a list of nucleotide sequences having different sizes?

If your sequences are different lengths, this implies a gapped alignment with insertions/deletions. This is not possible with a regular PSSM, so you have to create a separate PSSM for each group of homologous sequences of the same sequence length.

The alternative is to use a different approach that allows gaps like profile HMMs (as implemented in HMMER3 or HH-Suite3).


Obtaining identifier from plot of sequence GC% I was looking at this script from https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc297 and wanted to use it to determine outliers for tentative Horizontal gene transfer.
The data set is provided in the example, under the section 20.2.2 Plot of sequence GC%.
I want to yield a identifier for the x axis instead of just an arbitrary number of genes or cds' so I can actually determine which sequences have the outlying values, as I can not with the current graph.

I am not great at python or understanding list comprehensions, so if anyone could augment the code provided in the example to do so, it would be appreciated.
I provided the code below also.
#from : https://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc297, data set under
#section 20.2.2 Plot of sequence GC%, ls_orchid.fasta

from Bio import SeqIO
from Bio.SeqUtils import GC

gc_values = sorted(GC(rec.seq) for rec in SeqIO.parse("ls_orchid.fasta", "fasta"))
import pylab
pylab.plot(gc_values)
pylab.title("%i orchid sequences\nGC%% %0.1f to %0.1f" \
            % (len(gc_values),min(gc_values),max(gc_values)))
pylab.xlabel("Genes")
pylab.ylabel("GC%")
pylab.show()

The tricky bit here is that you need a list of sequence records sorted in the same order as the CG values. You can accomplish this by making tuples of the seq records and GC values, then sorting the tuples specifying the GC values as the sorting key:

gc_w_records = sorted(((GC(rec.seq), rec) for rec in SeqIO.parse("ls_orchid.fasta", "fasta")), key=lambda tup: tup[0])
gc_values = [x[0] for x in gc_w_records]
records = [x[1] for x in gc_w_records]

I don't know much about matplotlib and it's going to be challenging to display 100 IDs, but you'll probably be interested in the extremes of the GC distribution, which you can obtain directly from the sorted list of records:

records[0:2]
# [SeqRecord(seq=Seq('CGTCACGAGGTCTCCGGATGTGACCCTGCGGAAGGATCATTGTTGAGATCACAT...CAT', SingleLetterAlphabet()), id='gi|2765587|emb|Z78462.1|PSZ78462', name='gi|2765587|emb|Z78462.1|PSZ78462', description='gi|2765587|emb|Z78462.1|PSZ78462 P.sukhakulii 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]),
# SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACAT...AGG', SingleLetterAlphabet()), id='gi|2765568|emb|Z78443.1|PLZ78443', name='gi|2765568|emb|Z78443.1|PLZ78443', description='gi|2765568|emb|Z78443.1|PLZ78443 P.lawrenceanum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])]

or:

records[-3:-1]
# [SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGACAGCAG...TAA', SingleLetterAlphabet()), id='gi|2765656|emb|Z78531.1|CFZ78531', name='gi|2765656|emb|Z78531.1|CFZ78531', description='gi|2765656|emb|Z78531.1|CFZ78531 C.fasciculatum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]),
# SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGCTCATTGTTGAGACCGCAA...AAG', SingleLetterAlphabet()), id='gi|2765625|emb|Z78500.1|PWZ78500', name='gi|2765625|emb|Z78500.1|PWZ78500', description='gi|2765625|emb|Z78500.1|PWZ78500 P.warszewiczianum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[])]


DMPfold error while running protein structure prediction I am using DMPfold to predict protein structure of a 382 residue protein sequence. All goes wel, and even the files are produced but during the "seq2maps.csh", I get this error in between:

Running PSIPRED & SOLVPRED
Running PSI-BLAST with sequence c600m2-tail-before-lysin-sequence.temp.fasta ...
Segmentation fault (core dumped)
[blastpgp] WARNING: Unable to open BLOSUM62
[blastpgp] WARNING: BlastScoreBlkMatFill returned non-zero status
[blastpgp] WARNING: SetUpBlastSearch failed.
BLASTP 2.2.17 [Aug-26-2007]

Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.

Reference for composition-based statistics:
Schaffer, Alejandro A., L. Aravind, Thomas L. Madden,
Sergei Shavirin, John L. Spouge, Yuri I. Wolf,
Eugene V. Koonin, and Stephen F. Altschul (2001),
"Improving the accuracy of PSI-BLAST protein database searches with
composition-based statistics and other refinements",  Nucleic Acids Res. 29:2994-3005.

Query= C600M2_00100 tail protein [Salmonella phage BPS17S6]
         (377 letters)

Database: c600m2-tail-before-lysin-sequence.a3m
           46,978 sequences; 11,323,713 total letters

SearchingFATAL: Error whilst running blastpgp - script terminated!
Running PSICOV
Running FreeContact
Running CCMpred
 _____ _____ _____               _
|     |     |     |___ ___ ___ _| |
|   --|   --| | | | . |  _| -_| . |
|_____|_____|_|_|_|  _|_| |___|___|
                  |_|

using CPU (10 thread(s))

Reweighted 30404 sequences with threshold 0.8 to Beff=22917.9 weight mean=0.753779, min=0.0232558, max=1

Will optimize 62686429 32-bit variables

iter    eval    f(x)        ║x║         ║g║         step
1       1       1.50338e+07 152357      3.021629e+11    1.55e-06
2       1       1.47369e+07 152356      2.1447226e+11   1.18e-06
3       1       1.44292e+07 152355      1.4646483e+11   1.22e-06
4       1       1.41268e+07 152353      1.0388224e+11   1.5e-06
5       1       1.38295e+07 152358      7.337667e+10    1.92e-06

Can you please let me know if I need to rerun seq2maps.csh because of this error? And, how this error can be solved? Thank you!

The comment by Devon Ryan is correct - looks like your blastpgp command is not working properly. For a protein of that length it's more likely to be a problem with the BLAST executable running on your system than a memory issue.

Try running the blastpgp command (https://github.com/psipred/DMPfold/blob/master/bin/runpsipredandsolvwithdb#L26) on its own and see what happens:

~/path/to/blastpgp -a 12 -b 0 -j 3 -h 0.001 -d yourfile.a3m -i yourfile.fasta


Understanding Single Cell RNAseq Plots I have just started some single cell RNAseq analysis with 500 cells. I have used Seurat for analysis. Following the guideline I have generated few plots. Now I am trying to understand how to read those plots. Can the community provide their input in the following plots. To give a brief background of the experiment. This data was processed in MinION. Its from human brain sample from our collaborator. In their shortseq data they identified 15 clusters. However, both the libraries were processed independently. The main question I am trying to answer is celltype specific isoforms. but in my dataset I can find only 6 Cell clusters. So, I am trying to understand why I am getting such a small cluster? Is it because I stared with only 500 Cells? Also I want to know what these two PCA plot is indicating. Thanks

Identifying 15 clusters from 500 cells is challenging. It will depend on the relative abundance of each cell type within the tissue. You can use the cell type frequencies found in the shortSeq data set to see what to expect from a 500 cell dataset using the howmanycells tool

About the pc plots:
The first is an elbowplot which explains how many variation is explained by each PC. The knee basically defines when components become less informative.

EDIT: I noticed that the y-axis of the elbowplot starts at 4 suggesting that the dimensionality of your data is most likely ~10 dimensions. But I expect many of the clusters that don't show up on your Umap plot to be present at very low frequencies.

The second plot is a jackstraw plot and provides a visualization of a statistical test that calculates the significance of your principal components. This one shows that all your components are significant. You can combine this with the elbowplot to estimate the dimensionality (10-15 dimensions) of your data and decide which components to use for downstream analysis. In general the more they dropoff to the dotted line the less significant they become.

While these plots are useful as a first impression, always rely on biological information and common sense to decide what information (here: components) to include in downstream analysis such as Umap, clustering, etc. For example: usually you don't want to include components mainly defined by ribosomal, mitochondrial or hsp genes.


It means -inf. It has to do with how the PSSM score is calculated.


truncate fastq.gz file without error in Snakemake All,

As the first step in my Snakemake pipeline, I generate small fastq files from all input fastq.gz files. This is to run a quick round of "infer_experiment.py" (from RSeQC, to get the kit strandedness). However, Snakemake says this step fails. I think I know why, it is because a command like this

zcat fastq_raw/20190011-043_S43_R1_001.fastq.gz | head -n 400000 > infer_strandedness/20190011-043_S43_subset_R1.fastq

Ends with an error: gzip: stdout: Broken pipe and Snakemake uses Bash strict mode. (Snakemake reports: (one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)). So my question is:

How to avoid this, or, are there other truncating methods that do finish cleanly? Preferebly they are also efficient (seqtk offers subsampling but this parses the entire fastq file!).

Snakemake is using bash strict mode to ensure best practice error reporting in shell commands. Unfortunately some tools violate bash strict mode.

More details of the error, here. Alternatively, you can use if you don't want disable the strict mode.

zcat file.fastq.gz  | awk '(NR<=400000)' > subset.fastq


What causes the difference in total length of assembled contigs and scaffolds in SOAPdenovo2? I use SOAPdenovo2 to assemble a large genome (4.8G) using ~20X paired-end reads. The  total length of contig sizes is 6.3G while total length of scaffolds is 2.7G. Note that this is a haploid genome, so there is no issue of heterozygosity for scaffolding.

I am wondering what happened during the scaffolding procedure? How comes that sum of scaffolds is so much smaller than sum of contigs?

Something else that could be happening is that contigs that are being collapsed into "heterozygous" groups. This would be a particular problem when a genome has a substantial amount of repeated sequence. Digging deep into the supplementary information of the SOAPdenovo2 paper, I've found the following information:

  In SOAPdenovo2, heterozygous contig pairs are recognized by utilizing the information of contig depth and the locality of contig. The recognized heterozygous contig pairs should obey the following rules: 1) the similarity between contigs should be high enough, for example, ≥ 95%; 2) the depth of both contigs should be near half of the average depth or all contigs, complying Poisson distribution; 3) the two contigs should be located adjacently in a scaffold and have no relationship to each other inferred by paired-end reads information. The normal contigs neighboring the heterozygous regions, if they exist, could be connected to both of the heterozygous contig pairs (H1 and H2). Only the contig with relatively higher depth in a heterozygous contig pair were kept for scaffolding. The method reduces the influence of genome heterozygosity on final scaffold length. All heterozygous contig pairs were outputted to a file to facilitate further analysis. However, the trade-off of this method is that it might incorrectly remove paralogous contigs. This problem could be relieved by a gap-filling procedure while the removed copy of paralogous contigs would be represented by gaps during scaffolding.


How to select the most representative pathways from a gene enrichment analysis? I have perform an enrichment analysis to a cluster of genes. The output is a list of pathways and their p-value (the pathways are selected because p-value < 0.05). The list is still quite long, so I want to reduce it. For that purpose I have a calculated the Dice coefficient of the pathways in a matrix $p$x$p$ where $p$ is the number of pathways in the list. I want both the ones that are more different (they overlap less, their Dice coefficient is lower) and the pathways more representative of the most similar pathways (So if a there is a group of 5 pathways that overlap over 0.8 take just one).

How can I select the most representatives pathways?

There is a similar tool for GO but it relays on discarding not significant GO, while here all the initial pathways are already significant.

If I do a clustering of the pathways using the Dice coefficient matrix I don't know where (or how) to cut.

I tried using the height to select the pathways. But I am unsure of the interpretation of height.

Some other tools I have seen use a multidimensional scaling plot, but I am not sure if performing it and cutting at certain point of the first dimension would help.

This sounds like something that might be amenable to a clustered heatmap plot, or a correlation matrix plot, or something similar. Have you looked at a correlation matrix of the dice coefficient matrix (or maybe just a heatmap plot of that matrix without the correlation matrix)?

The corrplot package looks like it might be useful, in particular the hclust / drawing rectangles presentation.

I can't vouch for this package though; it's just something I found by a search for "R plot correlation matrix".


Is it possible to create a DESeqDataSet with a user-provided design matrix? I'm trying to run a differential gene expression analysis using DESeq2, with counts coming from kallisto. I have imported them using tximport and I'm creating the DESeqDataSet (dds) using the DESeqDataSetFromMatrix function.

> dds <- DESeqDataSetFromMatrix(counts,
                                samples,
                                ~ strain + batch)

And I get the following error, expected given my experimental design:

Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

Now, I know that I can just remove one column from the design matrix to make it work, but is there a way to supply my own design matrix to DESeq2? The following code raises an error:

> design <- model.matrix(~strain+batch, s2c)
> design = design[, -9] # removing a column to avoid full-rank
> dds <- DESeqDataSetFromMatrix(counts, s2c, design=design)
converting counts to integer mode
Error: $ operator is invalid for atomic vectors

Is there a way to provide my own model.matrix?

p.s. the modified model works in sleuth, but I would like to use DESeq2 for this particular analysis.

Provide rank sufficient design to DESeqDataSetFromMatrix and then use your custom model matrix in DESeq. In essence:

dds = DESeqDataSetFromMatric(counts, s2c, design=~batch)
design <- model.matrix(~strain+batch, s2c)
design = design[, -9]
DESeq(dds, full=design)

See this thread on the bioconductor site for details.


Loading local FASTA file in igv.js I'm trying to use igv.js for a very simple visualization of some genomic alignments, but I can't get past the first step.  I have a FASTA file that I want to use as the reference, but igv.js consistently gives an error when I try to load that file.  The file is a valid FASTA (have tested others from other sources to be sure - same error), and I can specify a remote https:// file, but not a local file.

Here's the HTML page I'm working on:

<!DOCTYPE html>

    <meta charset="UTF-8">
    <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no">

    <title>Local FASTA Test</title>

Local File Test

and here's the FASTA file

>HSBGPG Human gene for bone gla protein (BGP)
GGCAGATTCCCCCTAGACCCGCCCGCACCATGGTCAGGCATGCCCCTCCTCATCGCTGGGCACAGCCCAGAGGGT
ATAAACAGTGCTGGAGGCTGGCGGGGCAGGCCAGCTGAGTCCTGAGCAGCAGCCCAGCGCAGCCACCGAGACACC
ATGAGAGCCCTCACACTCCTCGCCCTATTGGCCCTGGCCGCACTTTGCATCGCTGGCCAGGCAGGTGAGTGCCCC
CACCTCCCCTCAGGCCGCATTGCAGTGGGGGCTGAGAGGAGGAAGCACCATGGCCCACCTCTTCTCACCCCTTTG
GCTGGCAGTCCCTTTGCAGTCTAACCACCTTGTTGCAGGCTCAATCCATTTGCCCCAGCTCTGCCCTTGCAGAGG
GAGAGGAGGGAAGAGCAAGCTGCCCGAGACGCAGGGGAAGGAGGATGAGGGCCCTGGGGATGAGCTGGGGTGAAC
CAGGCTCCCTTTCCTTTGCAGGTGCGAAGCCCAGCGGTGCAGAGTCCAGCAAAGGTGCAGGTATGAGGATGGACC
TGATGGGTTCCTGGACCCTCCCCTCTCACCCTGGTCCCTCAGTCTCATTCCCCCACTCCTGCCACCTCCTGTCTG
GCCATCAGGAAGGCCAGCCTGCTCCCCACCTGATCCTCCCAAACCCAGAGCCACCTGATGCCTGCCCCTCTGCTC
CACAGCCTTTGTGTCCAAGCAGGAGGGCAGCGAGGTAGTGAAGAGACCCAGGCGCTACCTGTATCAATGGCTGGG
GTGAGAGAAAAGGCAGAGCTGGGCCAAGGCCCTGCCTCTCCGGGATGGTCTGTGGGGGAGCTGCAGCAGGGAGTG
GCCTCTCTGGGTTGTGGTGGGGGTACAGGCAGCCTGCCCTGGTGGGCACCCTGGAGCCCCATGTGTAGGGAGAGG
AGGGATGGGCATTTTGCACGGGGGCTGATGCCACCACGTCGGGTGTCTCAGAGCCCCAGTCCCCTACCCGGATCC
CCTGGAGCCCAGGAGGGAGGTGTGTGAGCTCAATCCGGACTGTGACGAGTTGGCTGACCACATCGGCTTTCAGGA
GGCCTATCGGCGCTTCTACGGCCCGGTCTAGGGTGTCGCTCTGCTGGCCTGGCCGGCAACCCCAGTTCTGCTCCT
CTCCAGGCACCCTTCTTTCCTCTTCCCCTTGCCCTTGCCCTGACCTCCCAGCCCTATGGATGTGGGGTCCCCATC
ATCCCAGCTGCTCCCAAATAAACTCCAGAAG

It looks as though the script is looking for a URL (fastaURL) so you could try specifying the file as file://directory/test.fasta

By default though, javascript cannot read local files from disk for security reasons. Are you sure it is even possible to specify a local file in your script in this way?  If yes, then the problem is likely that you need to specify the path correctly.


Given a Genomic Ranges of SNPs, how to inject these SNPs in genome via BSGenome? Let's say I have the genome hg19 loaded into R via BSGenome

library("BSgenome")
hg19genome = getBSgenome('BSgenome.Hsapiens.UCSC.hg19', masked=FALSE)

I then have a list of SNPs loaded as a GRanges object, gr

library(GenomicRanges)

 > gr
 GRanges object with 212 ranges and 3 metadata columns:
        seqnames               ranges strand |     width     REF    ALT
           <Rle>          v <IRanges>  <Rle> | <numeric>
    [1]        1 [86099032, 86099032]      * |         1     C      T
    [2]        1 [86099033, 86099033]      * |         1     C      A
    [3]        1 [86099199, 86099199]      * |         1     G      A
  ....

Is there a straightforward way to inject these SNPs into hg19genome?

This seems relatively complicated given the structure of a BSGenome object.

The creator of the package answered this question previously on the Bioconductor support forums:
https://support.bioconductor.org/p/86665/#86757

  We don't provide an easy way to inject arbitrary SNPs in an arbitrary
  BSgenome at the moment. However, it should not be too hard to forge
  the BSgenome ... package. First you would need to
  compute the sequences of the mutated chromosomes (you can
  use replaceLetterAt for this), then write them to a 2bit file (put
  them in a DNAStringSet object and call rtracklayer::export on it),
  then use that 2bit file to forge the BSgenome ....
  package (see the BSgenomeForge vignette in the BSgenome package for
  how to do this).
  - Herve Pages

Here is another similar question with more info about replaceLetterAt: https://support.bioconductor.org/p/26199/


Location of Reads in Sequencing I have questions regarding sequencing:

Are paired-end reads from different stands?
What about single reads, are they coming from both strands?

Mates in a pair are giving signal for a single strand, though this is only apparent in bisulfite-treated samples. The reason for this is that you are sequencing the ends of fragments of material and it's random which strand of this (or both in different clusters) will get sequenced. Different pairs of reads will then randomly sequence one strand or the other.

Single-end reads are the same as paired-end reads, the second read in the pair is just not made.


How to interpret methylation calls from Bismark on opposite strands? I'm looking at Reduced representation bisulfite sequencing (RRBS) data from ENCODE, and to align the FASTQ files I've used Bismark with Bowtie 1. When I load the resulting BAM file into R with Rsamtools and GenomicRanges, I get something like this:

> reads
GRanges object with 2 ranges and 5 metadata columns:
      seqnames         ranges strand |                     seq                                   XM          XR          XG        NM
         <Rle>      <IRanges>  <Rle> |          <DNAStringSet>                          <character> <character> <character> <integer>
  [1]     chrM [16523, 16558]      - | ATAAAACCTA...CCCTTAAATA .....h........h.......z.............          CT          GA         3
  [2]     chrM [16524, 16559]      + | TAAAGTTTAA...TTTTAAATAA .....hh.......hhh.h.z...hhhh........          CT          CT        11
  -------
  seqinfo: 25 sequences from an unspecified genome

For the read on the - strand (which is aligned to the "G->A"-converted reference), how should I compare methylation calls to those on the + strand, since they don't line up?

For example, the z in the methylation string of the - strand is one position ahead of the z in the + strand (which makes sense because of symmetric CpG methylation). But how should I determine whether these two methylation calls are "essentially the same" or not?

The strand the bismark reports is related to the strand from which the read originated, not necessarily how it's aligned. So, alignments on the + strand shouldn't have calls overlapping those on the - strand, since you can't have a C in the same place on the same strand. One should often see Z/z next to each other on opposite strand, like in your example, since these are CpG (so it should be a Z/z on the + strand and then a Z/z on the following base on the - strand). Thus, in your example you have two reads supporting unmethylation for the CpG as a whole (1 read for each of the Cs).

The most confusing thing about BSseq is that reads (or paired-end reads) are only ever informative for a single strand. This is actually true for all sequencing with Illumina instruments (single-stranded fragments are loaded, after all), but since strands are almost always complementary we can usually ignore this fact.

As an aside, you might find MethylDackel useful (full disclosure, I wrote it). It'll be much faster at extracting methylation calls than bismark and supports nice things like excluding regions of bias methylation and likely variant positions.


using blastn to search [error] so I'm fairly new to using linux and blast - have encountered an error in my work.

I have created a fasta database with the results in this directory path:
~/class/hw/db

I'm currently trying to run a search against my created database against query sequences in a file called seqs.fas.

My code is as follows:
blastn -db $db -query seqs.fas -out q1results.txt

However I keep getting this error:
Error: Too many positional arguments (1), the offending value: seqs.fas

Error:  (CArgException::eSynopsis) Too many positional arguments (1), the offending value: seqs.fas

I'm not sure what I'm doing wrong here, and would appreciate any advice to fix this, thanks!

I remember having this problem in the past, it was encoding related. Rewriting the -s manually (the dash sign) should fix the issue, please see this.

  Problem is a copy & past from the BLAST documentation, all dashes
  ('–') need to be replaced by hand at the command line.


Total read count of shotgun affected by few highly abundant sequences I am building statistical models to analyse output from Illumina shotgun sequencing (HiSeq 4000) on stool samples (but RNA-seq data should behave similarily). The raw counts are a statistical sample of the total collection of DNA molecules used as input for the sequencing machine, and I think we can assume that the sequencing platform has limited capacity, such that detected reads are always a small fraction of the the original sequencing library.

Consider two different samples: one control sample and one treatment sample, which is identical apart from a single species (or gene) that has much higher absolute abundance relative what it in the control. As a result, the reads of this species will be much higher in the treatment sample, and also the 'real estate' of the remaining genes in that sample will be decreased. This is simply the nature of relative abundances.

However, size factor estimation that uses the "median ratio method" (http://dx.doi.org/10.1186/gb-2010-11-10-r106) for instance, are based on "experience with real data [...] [that] shows [that] a few highly and differentially expressed genes [or abundant species] may have strong influence on the total read count." (http://dx.doi.org/10.1186/gb-2010-11-10-r106). I am curious to see experimental evidence for this last statement. Or does someone have at least an explanation of how total read count can be affected by the absolute abundance of a species/gene?

That sentence in the Anders and Huber paper is a bit misleadingly worded (I don't think that was their intent). What they mean is that highly expressed and highly differentially expressed genes can out-compete the majority of the other genes, which thereby leads to them appearing to have lower expression if one uses total read counts for normalization (basically, what you said in the second paragraph).

An alternate interpretation of what they wrote (likely the one you understandable understood to have been their meaning) is that some highly expressed genes will inflate the total number of reads a sample actually gets. I have yet to see this in practice, but it certainly COULD happen. There are technical biases in sequencing such that some sequences (often due to length and GC content) are easier to sequence on common machines. A good example of this would be PhiX, which forms clusters much more easily on Illumina machines than real samples. I expect that this is relatively rare in practice, but using a more robust normalization would still help in protecting against it.


Expression of a gene in different groups I would like to check the expression of a gene in different groups like Disease vs Normal samples. I want to make a plot out of that to check whether it is significant or not.

From this paper lncRNA I see that Figure 1C they used RPKM value. But, I'm using read counts data. My question is can I use log cpm to make a plot like that or should I need to use only RPKM?

You should never use RPKM. It’s simply obsolete in the age of paired-end sequencing, and has been replaced by FPKM (which is, strictly speaking, a synonym).

The linked blog post explains more generally the problems that measures such as FPKM and CPM suffer from. A more robust measure is the TPM (transcripts per million), which scales CPM by the lengths of the transcripts.

Lastly, none of these measures perform very well for comparing gene expression across replicates. Since you have read count data, you should run DESeq2 or edgeR to perform differential expression analysis across all genes, and go from there.

As the comment noted, the paper you cite does not represent scientific best practice. In fact, a competent reviewer would have rejected this analysis.


Smallest group size for differential expression in limma (bulk RNA-Seq) I am reading Smyth et al. (ref. 1). I want to run differential expression analysis on a bulk RNA-Seq dataset in which each group is composed by 2 samples. In the paper previously cited it is written that:

  Genes must be expressed in at least one group (or in at least three samples across the entire experiment, where three was chosen as this is the smallest group size) to be kept for downstream analysis.

Is it possible to use limma DE analysis also with groups composed by only 2 samples? NB. It is possible that in this particular dataset the smallest group size is 3.

If not, which alternative should I use?

Update

I have 17 samples. My idea is to test one group vs the rest of groups (all together). So at least the reference group would be composed by > 10 samples. In this case, what can be a feasible analysis for DE?

Smyth, G. K., Law, C. W., Alhamdoosh, M., Su, S. & Ritchie, M. E. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. F1000Research 5, 1408 (2016).

The problem with most of the methods is that use the gene's variance for each group, which can't be calculated (reliable) when the sample is <= 2. Also, statistically, it would have extremely low power, so the conclusions couldn't be trusted much.

You could also calculate the "raw" fold change by yourself (ie without the variance estimation and adjusting). If you had one more sample per group you could use DESeq2 which was thought for such cases. But from the paper:

  However, if there are two or fewer replicates for a condition, these samples do not contribute to outlier detection, as there are insufficient replicates to determine outlier status.

Which makes me unsure if it will work well for just two samples as it won't help to detect outliers genes.


How do I create Protein Ligand images I was wondering how I can create images showing the bonds and the residue codes involved in the bound Protein Ligand complex obtained after Molecular Docking as depicted in the image below.

Pymol, in either its commercial or open source edition, will be sufficient to make figures such as your examples. Pymol recognises "polar contacts" for the intermolecular interactions and has labels that can be used for residue codes.
It is also scriptable and is widely used.


.cif .pdb conversion with python can you please give me an advice how to convert .cif files into .pdb preferably using python?
Thanks in advance!
Best,
Balint Biro

As @Michael said PyMOL can do this. PyMOL is not only an app, but also a python package —installed via conda not via the regular download (you can have both).
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load(infile,'myprotein')
     pymol.cmd.save(infile.replace('.cif', '.pdb'), selection='myprotein')

This runs the parallelisable instance, so you can have that context manager block running on how many threads or processes you fancy.


Where are .motif files from homer knownResults? I have been using homer's findMotifsGenome.pl, but with my new version (v4.9.1) of homer I don't get .motif files in the knownResults folder. I do get them in the homerResults (de novo) folder, though. With my previous version of homer I did get them and used them to find which peaks are in which motif. Which is what I would like to do

Does anyone know why there are no .motif files? I have only .svg files. I have seen that for the homerResults the .motif files are concatenated already into a file called homerMotifs.all.motifs, but no such file for knownResults.

PS. I tried to call Chuck for help, but he was not answering his phone...

PS2. I will ask this question on Biostars to hopefully get more response than here.

I have provided a solution in Biostars:

If you edit the file findKnownMotifs.pl and move the call to printMotif to come before testing $svgFlag you will find that .motif files are printed again.

In the following I have commented out the call from where it need not be, and introduced it as the first line, right before the if block:

   printMotif($_, $directory . "/$knownDir/$matrixFile"); # this should be done regardless of $svgFlag
    if ($svgFlag) {
        my $str = HomerSVGLogo::motif2svg($_,$bitFlag,$rnaFlag2);
            open SVG, ">$directory/$knownDir/$logoFile.svg";
            print SVG $str;
            close SVG;
            print MAIN "$str\n";
    } else {
  #printMotif($_, $directory . "/$knownDir/$matrixFile"); # this should be done regardless of $svgFlag
            `profile2seq.pl "$directory/$knownDir/$matrixFile" 100 > "$tmpFile"`;
            my $width = $widthFactor * $_->{'len'};
            `seqlogo -a $bits -f "$tmpFile" -F PNG -S -c -o "$directory/$knownDir/$logoFile" -h $height -w $width`;
  `seqlogo -a $bits -f "$tmpFile" -F PDF -S -c -o "$directory/$knownDir/$logoFile" -h $height -w $width 2> /dev/null`;
  print MAIN "<IMG src=\"$knownDir/$logoFile.png\"/>";
}

I expect this change was mistakenly introduced while refactoring the code to support svg as alternate output format for logos.

I have reported this to the author of this code.


remaining human genome variation that hasn't been sequenced? Given all the genome variation information we have accumulated from resequencing human genomes up until now (early 2018), how much variation is left that hasn't been sequenced?

EDIT: given one of the answers, and the fact that 99% of the SNPs at 1% or more are found, is there a quantitative measure of remaining variation?

In theory, almost any base in the human genome may mutate, so you have billions of variants to go. Ok, this is not so useful.

A related and potentially useful question is: given a human, what is the fraction of his/her variants seen in the the previously sequenced samples? For a Wright-Fisher population, there is an analytical answer: x% of variants on a haplotype has frequency below x%. Replace "x" with "1": 1% of variants on a haplotype have frequency below 1%. The 1000 genomes project concluded that 99% of SNPs at frequency 1% or higher have been called by the project. If we all came from a Wright-Fisher population, 99% of of SNPs have been called in 1000g. Given that the project also calls SNPs at lower frequency, nearly all SNPs in a newly sequenced sample are present in 1000g.

The above is a theoretical analysis. In practice, there are two major complications. First, we are not a Wright-Fisher population. Due to recent population expansion, each haplotype harbors more new mutations than the Wright-Fisher model would predict. Given a known population history, it is actually possible to numerically compute $f(x)$, the fraction of variants with frequency below x (for Wright-Fisher, $f(x)=x$). I don't have this result, though. Second, 1000g and many other genome projects didn't call SNPs in repetitive regions and missed a significant portion of high-frequency indels and structural variations. With Illumina, you can only call ~30% of long deletions callable with PacBio data, but we only have a dozen of PacBio human genomes in public. There are still a lot we don't know.


WGCNA module preservation error code with plotting summary statistics I am trying module preservation with WGCNA using my own data. I am plotting the summary statistics including Zsummary, medianRank, observed IGP and IGP permutation p-value for two comparison data sets.

When plotting the scatterplots for IGP & p-value vs. module size I get following error code. I have re-run code and done practice data with no problems but unsure how to correct for this error? Does anyone have any experience with this and know how to resolve? Thank you kindly

 plotData[[2]][, 3] = cr[[ref]][[test]]$Acutal.IGP
  plotData[[2]],[, 4] = -log(cr[[ref]][[test]]$p.value + 1e -4)

errror in plotData[[2]][, 3] = cr[[ref]][[test]]$Acutal.IGP
number of items to replace is not a multiple of replacement length`````

My first guess is that there is a typo in your code, and Acutal.IGP should be Actual.IGP. I checked my tutorials and the word Acutal is in none of them, so I would assume you mistyped the variable name.


two aligners combined results I would like to ask if someone uses two different aligners to produce count matrices and then run the same alogrithm for DEG analysis, would it make sense to find the intersection of the DEGs in order to conclude to the "real" DEGS?

This sort of procedure is done in tool papers to benchmark results compared to different aligners/mappers. See the salmon or kallisto papers for examples. Using the intersect of DE results from different mappers will end up yielding somewhat conservative results, not unlike using different differential expression packages and taking the intersect of their outputs. In practice you could just take a more conservative p-value cut-off and save the hassle of using multiple tools, the results would end up being similar.


Remove Redundant Sequences from FASTA file in Python I'm attempting to remove redundant sequences from a fasta file (from NCBI). When I execute this code, it returns the number of spots, not the number of sequences. (Number of spots: 408,293, Number of sequences: 40,830)
Please let me know if you have run a code that can remove the duplicate ideas or have any ideas/suggestions.
Thank you SO much.
from Bio import SeqIO
import time

start = time.time()

seen = []
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
    if str(record.seq) not in seen:
        seen.append(str(record.seq))
        records.append(record)

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

print(f"Run time is {(end- start)/60}")

Your code (as currently formatted) should generate the correct output, but it is not efficient. Make seen a set to improve these O(n) in calls to O(1) time complexity- this will speed it up a lot
seen = set()
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
    if record.seq not in seen:
        seen.add(record.seq)
        records.append(record)


ATAC seq : Plotting insert sizes My tutorial link
library(GenomicAlignments)

atacReads <- readGAlignmentPairs("MACS2_peak_call/ALL_ATAC_data/HSC1.bam", param = ScanBamParam(mapqFilter = 1,
            flag = scanBamFlag(isPaired = TRUE, isProperPair = TRUE), what = c("qname",
        "mapq", "isize"), which = GRanges("chr20", IRanges(1, 63025520))))
# length(atacReads)
atacReads

atacReads_read1 <- GenomicAlignments::first(atacReads)
insertSizes <- abs(elementMetadata(atacReads_read1)$isize)
head(insertSizes)

library(magrittr)
library(dplyr)
library(ggplot2)
fragLenPlot <- table(insertSizes) %>% data.frame %>% dplyr::rename(InsertSize = insertSizes,
    Count = Freq) %>% mutate(InsertSize = as.numeric(as.vector(InsertSize)),
 Count = as.numeric(as.vector(Count))) %>% ggplot(aes(x = InsertSize, y = Count)) +
  geom_line()

fragLenPlot + theme_bw()

fragLenPlot + scale_y_continuous(trans = "log2") + theme_bw()

How to i specify the Grange to span the total chromosome. Now I use only "chr20" as shown in the example.

You can exclude the which statement and get the entirety of all of the chromosomes. Note however that R is a very bad platform to use for this, since you end up reading everything into memory (I hope the BAM file isn't huge).
Note that you're seriously reinventing the wheel. There are a number of packages that will directly produce the plot you want, such as bamPEFragmentSize from deepTools.


How to detect a mutation and predict its consequence? I read lecture notes about mutations, and wondered what kind of algorithms are there to detect mutations?
How do we know if the gene is mutated or whether it is a sequencing error?
I saw this thread which is related, but I am not sure how to use the CIGAR information, e.g., is it 100% accurate? What is the underline mechanism to detect mutation? Is there a way to predict what the mutation will cause?

I'll follow up to the great answer from Kamil S Jaron:

Regarding predicting what the variant ("mutation" is a very loaded term) will do, there are a variety of tools. Chief among these are annovar and VEP. The general idea behind these is to classify the variants according to their overlap with genes, which codons they change (if any), how big that change is (e.g., changes in charge are more likely detrimental) and so on. One could also consider conservation, since if a position is highly conserved then changes in it are more likely to be detrimental.

If you really want to predict how a variant will change a protein's function then that usually requires prior knowledge about the proteins in question. Eventually someone will use machine learning to cull the literature and provide good predictions, but I haven't seen that yet.


Which data is being used for violin plot? Sorry I just got totally confused conceptually.

If this is my raw count data

Likely, Seurat divides each value by sum of the column afterward times by 10000.

which gives so

                   s1.1        s1.2 s1.3       s1.4 s1.5
DDB_G0267178 0.00000000 0.009263254    0 0.01286397    0
DDB_G0267180 0.00000000 0.000000000    0 0.00000000    0
DDB_G0267182 0.00000000 0.000000000    0 0.03810585    0
DDB_G0267184 0.00000000 0.000000000    0 0.00000000    0
DDB_G0267188 0.02640801 0.000000000    0 0.01286397    0

Literally Seurat claims that final step would be taking natural log of the above matrix as I tested with

library(seurat)
mat <- matrix(data = rbinom(n = 25, size = 5, prob = 0.2), nrow = 5)
mat

    [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    1    1    1
[2,]    1    0    2    2    3
[3,]    1    1    1    1    2
[4,]    1    1    0    1    0
[5,]    0    0    0    1    3
mat_norm <- LogNormalize(data = mat)
mat_norm
[1,] 7.824446 .        7.824446 7.419181 7.014015
[2,] 7.824446 .        8.517393 8.112028 8.112028
[3,] 7.824446 8.517393 7.824446 7.419181 7.706713
[4,] 7.824446 8.517393 .        7.419181 .
[5,] .        .        .        7.419181 8.112028

My confusion is: seurat@data dose not give me natural log transformed data rather returns only data divided by column sum and scalled by 10000.
So when plotting a violin plot if data are log normalised plot should not look as below picture

When I used

VlnPlot(object = y, features.plot = "DDB_G0277853", x.lab.rot = FALSE, y.log = TRUE)

plot changed totally. Could someone please help me in getting idea about seurat@data that is not really log transformed and by which data violin plot is being produced??? So if a gene has 10000 read counts, without log transformation violin plot should be so

This is my seurat object

Those two graphs don't look radically different.   The y axis is definitely not log scaled on the top one, and it looks like it is scaled in the second one (y.log = TRUE)


Finding the nearest gene to a specified gene region I've got some data on genomic regions of interest in a bedgraph-like format from the Mus musculus [mm10] genome, and would like to find out the nearest gene to these regions.

Some of the regions will reside within a gene (in which case I'd like to know that they're within a gene), some regions may reside within multiple genes (in which case I'd like to know all gene names), and some regions are not within a gene (in which case I'd like to know the closest gene: the shortest sequence gap between either end of the gene and either end of the tagged genomic region. In this case I want to do an unstranded search, because it's possible that I might have regions describing suppressing transcripts. Here's some example data:

Ref         Start             End      fwd     rev
chrX     73716152        73716152        0      -2
chrX    167207094       167209162        0      -3
chrY     30770844        30772724        3       0
chr1     24613189        24613641        0       3
chr1     24613971        24614812        0       4
chr1     35788729        35789387       -5       0
chr1     93151586        93160225        0      -2
chr1    107511454       107525597       -2       0
chr1    130882067       130887413        0      -3
chr1    134182466       134189942        2       0
chr1    135818954       135832956       -2       0
chr1    149520498       149521157        0      -3
chr1    166166500       166166500       -2       0
chr2      5379059         5379374        0      -2
chr2     32379110        32381837        0      -2
chr2     32381854        32381858        0      -2
chr2     32384633        32384633        0      -2
chr2     32384636        32387773        0      -3
chr2    120539390       120563799        0      -2
chr2    125946785       125975527        2       0
chr2    127633311       127633684        0      -2
chr2    127634841       127656210        0      -3
chr2    164562717       164562717       -2       0
chr2    164562719       164562720       -2       0
chr2    164568498       164568499       -2       0
chr2    181188007       181188388       -2       0

I'd like to add another 'Gene' column to this table that gives the official names of the nearest gene(s) for each region defined by <Ref>:<Start>-<End>.

You're looking for bedtools closest, which thankfully defaults both considering overlapping genes and outputting all hits in case of ties.

  bedtools closest -a foo.txt -b genes.bed | awk '{foo = sprintf("%s:%i-%i", $6, $7, $8); OFS="\t"; print $1, $2, $3, $4, $5, foo}'

You'll need to strip the header and ensure the file is sorted. I've piped the output to awk to reformat things. The output then is something like:

chr1    24613189    24613641    0   3   chr1:24613188-24613971
chr1    24613971    24614812    0   4   chr1:24613973-24614651
chr1    35788729    35789387    -5  0   chr1:35806973-35810462
chr1    93151586    93160225    0   -2  chr1:93151348-93160948
chr1    93151586    93160225    0   -2  chr1:93151582-93160870
chr1    107511454   107525597   -2  0   chr1:107511422-107525598

Those results are for the Gencode M15 release (the BED file I used has transcripts rather than genes, so you'll get a bit different coordinates). As an aside, the gene name would normally be output as well, just add column 9 to the output.


Mathematical notation for formulating a rank score This is a cross-post from Mathematics forum. As no one has answered it yet, I fIgured post here as well.

I would like to describe a transform I used to rank my data points. I have recoded my variable columns from $10$ measurements with different range to $n$ (for each measurement) where $n \in \{0,0.5,1\}$.

For columns where measurements were expressed as categorical values I simply recoded the categorical value to either $0$, $0.5$ or $1$ and for continuous variables I recoded the lower quartile as $0$, interquartile as $0.5$ and upper quartile as $1$ and finally I summed the recoded values to produce a single score $\sum n_{i1,..i10} $ for each row.

I am trying to write my method and I'd like to have it in mathematical notation. I was wondering if anyone could help me with this!

Here is my solution to the problem. I am posting it here in case anyone else come up with the same idea but did not know how to formulate it in mathematical notation!

  $$ \psi_i = \sum_{i=1}^9 S_i $$

  where $S_i$ is the score function and defined as:

  $$
S_i = \cases{
  0 & \text{if $\theta(x_i) < Q_1$ or $\theta(x_i) = $ benign/neutral} \\
  0.5 & \text{if $Q_1 \leq \theta(x_i) < Q_3$ or $\theta(x_i) = $ possibly damaging/uncertain} \\
  1 & \text{if $\theta(x_i) \geq Q_3$ or $\theta(x_i) = $ damaging}
}
$$

  The $\theta(x_i)$ is the pathogenicity or conservation score for variant $x$ as defined by model $i$ and $Q$ denotes the quartile range for scores from M-CAP, CADD, GERP and Phylop models.


How does Cellminer's "Cross-correlations of transcripts, drugs, and microRNAs" work I have seen one of Cellminer tools. I am not sure how do they calculate the cross correlation of the genes, what does it actually mean? Based on what databases?

For example if I take their example input (HUGO format):

abcb1
BRCA2
CNBP

I get the following cross correlation matrix (based on what signals).

Identifier      abcb1   BRCA2   CNBP
abcb1           1      -0.142   0.069
BRCA2          -0.142   1       0.176
CNBP           0.069    0.176   1

In the article describing Cellminer says about the cross-correlations:

  ... a cross-correlation table of the resultant z-scores can be
  generated ...

On page 3503

And the z-scores are defined previously as:

  The tool output includes relative transcript intensity presented
  as z-scores...


Visualize similar interactions in two homologous networks Background
I have two networks with homologous genes from two different species.
To identify homologs I performed a functional comparison. One of the species is not well studied while the other one is the house fly (D. melanogaster). I used ClueGO to predict the functions of the genes in the D. melanogaster network.
I found a function of interest influenced by a set of genes and predicted that the same genes in the other network (thus far uncharacterized) should also perform a similar function.
Questions
Q1: How do I check if there are overlapping patterns between both the networks on Cytoscape?
Q2: There are more than 100 homologs in this list and I would like to minimize the list by only considering those homologous genes that are arranged the same as well.

An easy way to visualize this is to make a tab-delimited table that contains edge information for both species and both networks, color the edges depending on the species, and make the thickness and transparency dependent on the interaction strength.

For example, say you have a list of orthologs for the two species and just assign them a number or a name. Each edge has a "source" and a "target". The edge weight is the interaction strength.

(Note do not copy-paste the below table as the mixed tabs/spaces may make the graph topology incorrect). In my toy dataset the data look like this:

source  target  variable    value
gene1   gene2   speciesA    3
gene1   gene3   speciesA    42
gene2   gene3   speciesA    3
gene1   gene2   speciesB    9
gene2   gene3   speciesB    5
gene2   gene4   speciesB    41
gene3   gene4   speciesB    32

Import your network to cytoscape as a network from a file:

Then go to the style tab:

... and make sure edge is selected.

Then you can adjust the transparency, edge color, and edge stroke settings, like so using edge stroke color:

:

... and transparency panels as examples.

Don't forget to save your style when it is complete!

When you have adjusted these settings you will have something that looks like this:


Why is there no GT-entry in my .vcf file from bacteria I mapped raw reads obtained from an E. coli genome (.fastq file) to my reference genome. Next, I exported the .vcf (variant calling file), which looks something like the header example shown below (total >50 lines):
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Variants:_LGC19-XL01_S63_L001_R_001_(trimmed)
CFC381_K12_Bw25113  360287  .   GGGAAT  AATTCC  1952.76 .   NS=1;VF=0.380;SB=0.533;SB50=0.60;SB65=1.0;TYPE=Substitution;AVQUAL=30   DP:AO   231:89

How can I back-fill or quickly assign the missing GTs to the .vcf ?
For my downstream process (plink --pca, via bfctools merge), I need that there are GT entries in every .vcf, so I can merge them and the respective GT is shown in the merged .vcf. I only found solutions that need the .bam file, but that is laborious and I don't think I need the .bam in the case of haploid data.

This method needs .bam: https://www.biostars.org/p/319304/
This method needs .fastq: https://training.galaxyproject.org/training-material/topics/variant-analysis/tutorials/non-dip/tutorial.html#examining-the-results

However, the.fastq is actually from a mixed culture, so I wonder do I really need to follow the .bam-strategy? Suggestion for a tool? But in that case, why did it not export the GTs in the fist place (I exported from "geneious" sequence analysis program)

The representation for a variant present in a haploid genome should be 1. A heterozygote in a diploid genome would be 0/1 (or 0|1 or 1|0  if phased). So, to add a 1 as the genotype, you would need to first insert a GT to the beginning of each INFO line and then the 1 at the beginning of the last field. But, this should happen only on lines that are not headers (do not start with #). This should do it:
$ awk -F'\t' -v OFS="\t" '{ if(/^#/){print}else{$9="GT:"$9;$10="1:"$10;print}}' file.vcf
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Variants:_LGC19-XL01_S63_L001_R_001_(trimmed)
CFC381_K12_Bw25113  360287  .   GGGAAT  AATTCC  1952.76 .   NS=1;VF=0.380;SB=0.533;SB50=0.60;SB65=1.0;TYPE=Substitution;AVQUAL=30   GT:DP:AO    1:231:89

So just redirect the output to a fixed file:
$ awk -F'\t' -v OFS="\t" '{ if(/^#/){print}else{$9="GT:"$9;$10="1:"$10;print}}' file.vcf > fixed.vcf

Now, in this specific case, you could also do something like this:
$ sed -i.bak '/^[^#]/{s/\<DP:AO\t/GT:DP:AO\t1:/}' file.vcf

But that assumes that all lines will have the exact same INFO field of DP:AO which is not a safe assumption since each variant can have different INFO fields.


How hard is it to clean and QC gene expression microarray data? I am a PhD student working on developing ideas for my dissertation papers. One of my planned papers will be working with some (human) gene expression data.

I have had a class on working with gene expression data, but no practical experience. From what I remember from my project in that class, you do a bunch of pre-checks and QC, then you remove batch effects, then you summarize expression across probes and normalize. I thank that part of my class project took me about three days.

But when I talk to other people about it, they say that cleaning/QCing/normalizing gene expression data is really hard, that it will take forever to get it working, and that if I do this dissertation paper I am going to end up taking extra years to graduate.

Am I missing something here? Why do the people I talk to seem to think this is so hard and time consuming -- is it really that bad? Should I avoid having a dissertation paper that uses gene expression data if I want to graduate on time?

(The data will be from Affymetrix microarrays. There are enough samples to use two 96-well plates, so batch effects will be more of a concern than they were in my class project.)

I am not sure how much you know about bioinformatics already, can you use R? For a bioinformatician looking at QC for microarrays should not be a big deal, at least for me it would take maybe a day (or two) to get this done. However, if you never used R and want to start from scratch, it depends on how quickly you learn how to deal with arrays and QC. There are many tutorials and review papers about normalization and QC of arrays. Try there and see if you can easily understand what these papers are telling you.

Read the limma user's guide, and try to follow the tutorials therein. If you understand that quickly you can try a real experiment of your own. If you don't understand what's described in there, you might want to look for another project for your thesis.


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

Yes, if you have a ranked list, you can use GSEA to calculate the enrichment of the strains along the ranked list.

Note that hypergeometric test is for count data, so you would test if there are more strains of one type that fit on certain conditions.

I would recommend this article comparing different enrichment methods to understand well what they do.


local BLAST error: BLAST Database error: Error: Not a valid version 4 database I want to use a local blastn command to BLAST a multifasta file ~700 sequences using following command: blastn -db nt -query fasta_all.fasta -num_alignments 2 -out fasta_blasted.txt but I receive error :BLAST Database error: Error: Not a valid version 4 database.

I use a local nt database which I downloaded today (11th of May 2020) by running update_blastdb.pl --decompress nt [*]. I downloaded this but I did not do any post-processing and just the "raw" download files are in the folder (see image below which shows only a part of the files)

The blast+ I downloaded and installed is ncbi-blast-2.10.0+-x64-linux.tar.gz from https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

So, both database and blastn version should be the latest. It probably has something to do with either the version of the database or the software, but as it both are the latest version, I find it strange that these would not be compatible. Any ideas what causes this error?

first I recommend you to check your downloaded data with md5 to be sure there was no problem with that step. How many files did you download?
Also I was curious if you had another BLAST+ version installed on your computer. Be sure that blastn is called is a version corresponding to the database you are using.


How do I write tests for a snakemake pipeline? I've been playing with snakemake, which is really powerful, however, I'd like to create a 'unit test' for the whole pipeline. e.g. a way to call the pipeline that tests that all the rules work together, all the steps complete, etc.
I assumed I could find some guides about how to do this, but I don't find any clear information about how to tests pipelines.
Any advice or suggestions would be welcome.

You can run CI on a snakemake workflow when provoding a minimal example dataset. Maybe an example from the snakemake-workflows repo will be helpful: https://github.com/snakemake-workflows/dna-seq-varlociraptor


Output of Seurat FindAllMarkers parameters I compared two manually defined clusters using Seurat package function FindAllMarkers and got the output:

Now, I am confused about three things:

What are pct.1 and pct.2?
How come p-adjusted values equal to 1? What does it mean?
If we take first row, what does avg_logFC value of -1.35264 mean when we have cluster 0 in the cluster column? Is it that in cluster 0 the Cttnbp2 gene downregulated by a factor of 2^1.35264?

pct.1– The percentage of cells where the gene is detected in the first group
p_val_adj– Adjusted p-value, based on bonferroni correction using all genes in the dataset. This is also known as a false discovery rate (FDR) adjusted p-value. An adjusted p-value of 1.00 means that after correcting for multiple testing, there is a 100% chance that the result (the logFC here) is due to chance.
You are correct. The expression of Cttnbp2 donwn-regulated by that factor. If you were to look in the table at the cluster 1 DE genes you would see Cttnbp2 with a logFC value of 1.35264. When you are looking at genes differentailly expressed in only two clusters, the sign is the only thing that changes.

Vignette posted above by @heathobrien.


awk for motif splicing from an amino acid alignment I'm looking to splice out specific amino acid alignment positions as a continuous motif from an alignment (~15 000 sequences). The code is written in Perl (I agree not many developers now) but splicing method is awk which a lot of you know. Please note the hash is "tied" (not shown), this just means the order of the sequences is preserved.  I can get the awk code below to return single amino acid positions.

I am struggling to get a loop working in awk.

I agree the Perl wrapper makes it more complicated, but its part of a huge data pipeline (all written in Perl).

#!/usr/local/ActivePerl-5.20/bin/env perl
#use strict;
use warnings;

# I loaded the alignment into a hash, where key = sequence and value = protein
my %dna = (
    virus1 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVAHLALIAAFKVRPALLVSFIF",
    virus2 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVXXXXXXAAFKVRPALLVSFIF",
    virus3 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGZZZZZZLIAAFKVRPALLVSFIF"
);

# I define the start and end of the alignment positions I am splicing
my $n = 1;
my $z = 5;

# I turn the alignment into a string with separated by carriage returns
my @array = map {"$dna{$_}"} keys %dna;
my $string = join "\n", @array;

# This is the awk loop, note in Perl the "``" means issue a shell command

my $splice = `echo "$string" | awk 'BEGIN {FS=""} {x = $n; while (x < $z); print \$ x; ++x}' 2>1`;
print $splice;
__DATA__

Expected output

LGVLV
LGVLV
LGVLV

Observed output
I get the original sequences returned

I am familiar with Python/Perl solutions of line by line processing but awk is more efficient IMO.

I think you are trying to use awk to get characters 3 through 7 (inclusive) of each input sequence. If so, a working version of your approach would be this:

my $splice = `echo "$string" | awk 'BEGIN {FS=\"\"} {x = $n; while (x < $z){ printf \"%s\",\$x; ++x} print ""}' 2>1`;

That, when added to your perl script, produces:

$ foo.pl
VLVIL
VLVIL
VLVIL

That isn't the output you show, but the output you show is just the first 5 letters of each sequence, which doesn't seem to correspond to what your script is doing.

Now, the reason your version didn't work, is that the while was like so:

while (x < $z); print \$ x; ++x}

Assuming the \$ x was a typo, the problem is that the ++x isn't part of the while loop. So that should never even have exited. Seeing as we're in the same time zone, and given the time, I am guessing you're tired and that you miscopied something when trying to give us a reproducible example. In any case, to get the loop right, you need {}:

while (x < $z){print \$x; ++x}

However, if all you want is a substring, a loop is a very inefficient way of doing that. You can instead use awk's substr command:

  substr(s, m[, n ])
             Return  the at most n-character substring of s that begins at position
             m, numbering from 1. If n is omitted, or if n specifies  more  charac‐
             ters than are left in the string, the length of the substring shall be
             limited by the length of the string s.

So, in your script, that would be:

my $splice = `echo "$string" | awk '{print substr(\$0,$n,5)}'

Since you're already using a language that excels in string manipulation, I really don't understand why you would try to use awk instead. It will not be faster or more efficient than the same thing in perl:

print substr($_,$n,5) . "\n" for @array;

By calling awk you are also calling the shell which will be slower by definition, especially if you're echoing thousands of lines. In fact, I tried this with 15000 sequences (just repetitions of your first one) and using either the awk loop I gave above, or the awk+substr, or the native perl substr:

awk loop, average time taken to run the same thing 100 times: 0.0525
awk substr, average time taken to run the same thing 100 times: 0.0522s
perl substr, average time taken to run the same thing 100 times: 0.0464

And you can make it even faster by skipping the array and iterating over the keys directly:

print substr($dna{$_},$n,5) . "\n" for keys(%dna);

That averaged at 0.0376s over 100 runs.


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

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

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

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


Verify a predicted protein in one genome in a different genome of the same species I have two genome assemblies of the same non-model species, call them Assembly 1 (generated from Illumina data) and Assembly 2 (generated from PacBio data).

For Assembly 1, I also have predicted proteome data, generated with EVM. Say there is a protein, call it Protein X, for which I have the fasta sequence from predictions, but suspect an error in Assembly 1 (because the Protein X has a long segment of auto-identity). For this reason, I would like to verify the assembly of the locus coding for Protein X in Assembly 2 and/or verify the presence of Protein X in the Assembly 2.

Data in my possession:

both genome assemblies (fasta)
raw rna-seq data
reference-based transcriptome assembly for each genome assembly (fasta/gtf)
set of predicted protein sequence based on Assembly 1 (fasta)

Things I was thinking of:

build the EVM predictions on the Assembly 2 and verify the identity of Protein X by simple fasta identity check
look for the protein sequence in the Assembly 2. I could use something similar to Reverse Translate to get the most likely nucleotide sequence from the input amino-acid sequence of Protein X, then align it to the Assembly 2 to check the coordinates of the locus it originates from - but I am not very sure of the relevance and/or correctness of this procedure
something allowing 'blast' of the protein sequence in the Assembly 2 (but I am not aware of any software that could do that)

Any ideas on how I could do that?

This should be very easy to do. Here are some options:

Use a tool like Exonerate or GeneWise that can align protein sequences to genomic DNA while attempting to model splice sites etc.
As you said, blast. There's nothing special in what you describe, really, this is what tBLASTn is for. Just use your protein as the query sequence and tBLASTn as the blast flavor which will take a protein query and match it against all 6 possible reading frames of a DNA database. Set up your assembly as the blast DB, and then blast as you would normally.

Of the two, using the first approach is much better since it will try to build a valid gene model for you instead of simply finding the regions of high homology as blast will do.

Finally, just for completeness, you can go old school and run a full de-novo gene prediction tool on your assembly. Something like GeneID for instance. But that's certainly overkill if you're just looking for a single gene.


How can I bypass 'out of memory' problem using a HISAT2 for human genome indexing? I am trying to perform a HISAT2 indexing using a GCF_000001405.39_GRCh38.p13_genomic.gtf and GCF_000001405.39_GRCh38.p13_genomic.fna files downloaded from https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.39 using the procedure described in  Pertea et al. (https://www.ncbi.nlm.nih.gov/pubmed/27560171), e.g. first I create extract_splice_sites.py and extract_exons.py and then my code looks like this

#!/bin/bash
#SBATCH --time=02:30:00
#SBATCH --account=def-myusername
#SBATCH --mem=160000M
#SBATCH --mail-user=mymail@gmail.com
#SBATCH --mail-type=All

module load nixpkgs/16.09  intel/2018.3
module load hisat2/2.1.0

hisat2-build --ss imam.ss --exon imam.exon GCF_000001405.39_GRCh38.p13_genomic.fna imam002

I am trying to receive a set of imam002*.ht2 files but every time the procedure ends with a Ran out of memory message. From here https://github.com/griffithlab/rnaseq_tutorial/wiki/Indexing I learned that I must allocate 160Gb of RAM if I want to index the entire genome. Even when I do that and submit my script via sbatch on a HPC I still got this

slurmstepd: error: Detected 1 oom-kill event(s) in step 30827017.batch cgroup.

Just in case here's the window that I have when I don't sbatch my script, I still have an error

Inactive Modules:
  1) openmpi/2.1.1

The following have been reloaded with a version change:
  1) gcccore/.5.4.0 => gcccore/.7.3.0     2) icc/.2016.4.258 => icc/.2018.3.222     3) ifort/.2016.4.258 => ifort/.2018.3.222     4) imkl/11.3.4.258 => imkl/2018.3.222     5) intel/2016.4 => intel/2018.3

Settings:
  Output files: "imam002.*.ht2"
  Line rate: 7 (line is 128 bytes)
  Lines per side: 1 (side is 128 bytes)
  Offset rate: 4 (one in 16)
  FTable chars: 10
  Strings: unpacked
  Local offset rate: 3 (one in 8)
  Local fTable chars: 6
  Local sequence length: 57344
  Local sequence overlap between two consecutive indexes: 1024
  Endianness: little
  Actual local endianness: little
  Sanity checking: disabled
  Assertions: disabled
  Random seed: 0
  Sizeofs: void*:8, int:4, long:8, size_t:8
Input files DNA, FASTA:
  GCF_000001405.39_GRCh38.p13_genomic.fna
Reading reference sizes
  Time reading reference sizes: 00:00:23
Calculating joined length
Writing header
Reserving space for joined string
Joining reference sequences
  Time to join reference sequences: 00:00:13
  Time to read SNPs and splice sites: 00:00:15
  Ran out of memory; automatically trying more memory-economical parameters.
./hisat2_index.sh: line 24:   943 Killed                  hisat2-build --ss imam.ss --exon imam.exon GCF_000001405.39_GRCh38.p13_genomic.fna imam002

When I use the same script for another genome (not human) - I don't have this problem.

Is it really a problem of RAM and, if it is, how is usually human genome gets indexed by HISAT2? Is it even possible?

You haven't allocated the requisite 160GB of RAM required. Or rather, you allocated 160GB of memory, which is what you were told to do, but they most likely meant 160GiB.

Computers use binary, not decimal, so a gigabyte is 1024 megabytes. This is the same annoying mix-up people run into when their hard drives don't have the space the marketing department says it does. By requesting 160000M in your job script, you actually requested a few gigabytes less than what your source says you need.

$$
160000 \text{ megabytes} \times \frac{1 \text{ gigabyte}}{1024 \text{ megabytes}} = 156.25 \text{ gigabytes}
$$

$$
160 \text{ gigabytes} \times \frac{1024 \text{ megabytes}}{1 \text{ gigabyte}} = 163840 \text{ megabytes}
$$

If "160GB" is what you need, then 163840M is the correct value for the job script.

I would honestly suggest cutting the Gordian knot and not dealing with this "gigabytes vs. gibibytes" ordeal by simply using the G suffix.1

#SBATCH --mem=160G


Scaling by linear regression against the number of reads I am trying to build the preprocessing pipeline presented in The Tabula Muris Consortium et al. (pp).

It is a pipeline to preprocess single-cell sequencing data. There is one step that is not clear:

  Counts were log-normalized (log(1 + counts per N)), then scaled by linear regression against the number of reads (or UMIs), the percent of reads mapping to Rn45s, and the percent of reads to ribosomal genes.

I understand the first part (I assume that log in this context is log2), but I need help on understanding how to scale by linear regression against the number of reads, the percent of reads mapping to Rn45s, and the percent of reads to ribosomal genes.

I don't know if this question has been solved already, but what they try to do is equalize the depth of sequencing for each cell. Therefore, they scale for the total number of reads. If you regress out (via linear or negative binomial regression) the differences in the number of reads per cell, you end up with cells that have been sequenced with the same depth.

In my opinion they use the ribosomal genes in the same way. They are considered to be a kind of housekeeping genes, which you can use to equalize the sequencing depth.

And yes, the log in this context is log2, used to get fold-change values instead of counts.


How to check if indels in VCF files are left or right aligned? I downloaded a VCF file from dbSNP, and I'm curious if the indels in the file are left-aligned for GRCH37 genome. The documentation doesn't say anything.

How can we tell if a VCF file has left or right aligned indels?

One possibility is to run bcftools norm on the file. At the end it will print out a statistic about how many variants were realigned. I did this for dbSNPv151 (GRCh38) only for chromosome 22 like this:

$ bcftools view -i 'TYPE="indel"' dbSNP_GRCh38.p7_b151.vcf.gz 22|bcftools norm -Ou -f GRCh38.fa > /dev/null
Lines   total/split/realigned/skipped:  901503/0/490/0


Prediction model for spread of coronavirus A previous thread on Bioinformatics StackExchange citing an article by the Economist which references a cartoon of predictive outcomes by the CDC (Economist article on coronavirus) suggests a model in epidemiology using the "basic reproductive rate", commonly known as R0. Although I'm no bioinformatician, I imagine that one could tweak a model based on similar infections like SARS and HDN1. Yet, I found no models predicting or giving confidence intervals for the spread of coronavirus. I wonder if the the data is still too sparse, e.g. we don't know if it's seasonal (mentioned in the above thread) and we cannot infer the virulence.

Would such a model be reliable? If not, how many months of data would we need to reliably use a model predicting the spread?

The original R0 calculation was developed in the early 1900s to model malaria transmission and is known as the Ross-MacDonald model so your idea is over 100 years old. The model was very famous for identifying the optimal strategy for malarial control by targeting the adult mosquito rather than the mosquito larvae. What it said was reducing the mosquito population and reducing malaria transmission were two different issues. So the bednet strategy is based in this idea and particularly spraying the inside of houses with insecticide.

R0 was famously used in the 1980s to predict the HIV pandemic when WHO thought it would remain localised in Africa.

The problem with R0 and Coronaviruses is we don't know the asymptomatic population who could also transmit the virus.

This population is potentially very large. Epidemiologists will attempt extrapolation from cruise ship quarrentine, eg the Diamond Princess where everyone is tested to understand this transmission parameter and project this onto the rest of the planet. Clearly this could have limitations or inaccuracies. The other issue I am baffled by is seasonal transmission and whether this will occur with COVID-19.

One thing R0 of COVID-19 will identify is the latent period transmission which occurs before the symptoms, i.e. the incubation period. Essentially, people with cough/mild symptoms before they self isolate and become really ill with a lung infection will transmit the virus. In this context R0 will be/was very influential for understanding when to move between quarrentine and containment phases of the pandemic.