What methods exist to calculate RNA expression profile similarity Some of the work in our lab requires a comparison of a strain across several experimental conditions. We are looking to identify most similar experimental conditions based on the gene transcription response similarity from the cell.

While we could easily invent and create home-grown methods to do it, their implementation and testing are a laborious project in themselves and are outside the scope of our current work.

Are there any methods for RNA expression profile similarity computation that has been already published? If yes, what is your experience using them?

There's a new paper that's just appeared on the subject of "proportionality", including a method by which RNA expression might be compared.

This is a new concept for me, and the article is not easy enough for me to read that I can write up a quick summary; the authors don't seem to devote a section in the paper to defining "proportionality". However, here's an interesting chunk from the article:

  We graphed the network of relationships between these mRNAs (S5 Fig.), an approach similar to gene co-expression network [12] or weighted gene co-expression analysis [13] but founded on proportionality and therefore valid for relative data. The network revealed one cluster of 96, and many other smaller clusters of mRNAs behaving proportionally across conditions.
  ...
  We are also keen to raise awareness that correlation (and other statistical methods that assume measurements come from real coordinate space) should not be applied to relative abundances. This is highly relevant to gene coexpression networks [12]. Correlation is at the heart of methods like Weighted Gene Co-expression Network Analysis [13] and heatmap visualization [14]. These methods are potentially misleading if applied to relative data.


Are there open-sourced graph-based variant callers? I would like to compare GATK with graph-based variant callers. I have seen the Fast and accurate genomic analyses using genome graphs paper from SevenBridges and the paragraph tool by Illumina, though it does not look like folks use their platform or it is not open-sourced (in the case of SevenBridges's GRAF).
The idea is to see if graph-based approaches detect better both SNPs and structural variants compared with existing methods. Additionally, I believe that if there exists some novel graph-based aligner/assembler one could use the metadata downstream to identify variants.

Yes, there are. There are some suggestions in the comments (VG, WhatsHap, GraphAligner, Minigraph).
However, to be clear, the current default variant calling algorithm used in the GATK (the haplotype caller) is a graph-based algorithm. It constructs a De Bruijn graph from mapped reads for local assembly around variants. Other tools that also follow this strategy include Platypus and freebayes.
You might also consider multi-color genome assemblers. For example, Cortex and Corticall use colored De Bruijn graphs to aggregate reads from multiple genomes and call variants.


Significance and timing of "mux scans" I'm using MinIONQC to do quality control on some ONT data. The software plots several read characteristics over time (hours passed during the sequencing process). These plots contain several vertical red lines. From the documentation:

  Muxes, which occur every 8 hours, are shown as red dashed lines

Reading this paper (doi:10.12688/f1000research.7201.1) I was able to get a rough idea of what a "mux" or a "mux scan" is. If I understand correctly, it's a calibration step that scans each well containing 4 nanopores to determine which nanopores will be used for data collection in which order.

Do I understand the basic idea of the mux scan correctly? Does it occur at the same frequency on all ONT instruments? Does it always occur in 8hr intervals, even if the sequencing reaction is run for longer than the standard run time?

Yes, your understanding is largely correct. This originates from the situation that for each detector on a nanopore array there are 4 pores. I'll explain mux scans and groups, but this is outdated information as now another system is used.

So a MinION FC has 2048 pores with 512 sensors. Not all of these pores will be equally suitable for sequencing. At the beginning of the run originally a mux scan was performed, to test for each sensor which the best pore would be (if any) out of the 4 available. The best pores will be used for sequencing first, the rest will be turned off. That is group 1. After 8 hours the expectation is that many of these pores have degraded/accumulated damage or are not functional for other reasons. Then it is time to swith to group 2. So you change to another set of pores, which in these cases often lead to a temporary increase in yield. The default script used 8 hours, although some used adaptive software to select new pores when the first group started going bad. Also note that during the run no new mux scan was performed: the software did not check if the selection it made in the beginning still made sense. Pores in group 1 got turned off after 8 hours and switched to group 2, without taking their current 'health' into account. Cumulative yield plots get a characteristic bumpy pattern from this. Each bump is less than the previous one, eventually flatting off.

Now that's history, as the software now uses dynamic muxing and will continuously check for the best pore per channel - and as such maximize yield by having as many healthy pores as possible available for sequencing. Cumulative yield plots are now smoother.


Using BLAT command line tool to blat split sequences I have the nucleotide sequence:

AATTGAGGCACATTTTTTTTTAGACAGTCTTGCTCTGTTGCCCAGGCTGGAGTGCAGTGGTGTGATCATAGCTCACTGCAGCCTCGACCTCCTGGGCTCAACAAAGCACACAGTGGGCGGATCCCCACCAG

When I blat this on UCSC Genome Browser, the first hit is a match which spans 6572 basepairs on chromosome 19, and essentially matches the first portion of the sequence to one spot on the genome, and the second portion to another part of the genome approximately 6.5kb downstream.

However when I use the BLAT command line tool, it doesn't find the second half of this match, it only matches the first portion.

Here are the flags I'm using:

-t=dna

-q=dna

-minIdentity=90

-out=blast8

How do I get the command line tool to mimic the behavior of BLAT on UCSC Genome Browser?

The FAQ offers an answer:

  I'm setting up my own Blat server and would like to use the same parameter values that the UCSC web-based Blat server uses.

  We almost always expect there to be some small differences between the hgBlat/gfServer and the stand-alone command-line blat. The best matches can be found using pslReps and pslCDnaFilter utilities. The web-based blat is tuned permissively with a minimum cut-off score of 20, which will display most of the alignments. Other than to confirm that your command-line blat is working, there is little use in perfectly replicating the web-based blat results. We advise deciding which filtering parameters make the most sense for the experiment or analysis. Often these settings will be different and more stringent than those of the web-based blat. With that in mind, use the following settings to replicate the search results of the web-based blat:

  [ . . . ]

  standalone blat:

  blat search:

blat -stepSize=5 -repMatch=2253 -minScore=0 -minIdentity=0 database.2bit query.fa output.psl


Why doesn't Biopython AlignIO.read() recognise the 'mauve' format? On this Biopython tutorial, they describe how to import a multiple sequence alignment in the Mauve (XMFA: extensible multi fasta format). So I imported the AlignIO module:

from Bio import AlignIO
alignment = AlignIO.read(open("alignment.xmfa"), "mauve")

but I get the following error:

ValueError
Traceback (most recent call last) <ipython-input-30-396e54791d7a> in <module>()
----> 1 alignment = AlignIO.read(open("/Users/cr517/Documents/alignment.xmfa"),
"mauve")

/Users/cr517/anaconda/lib/python3.6/site-packages/Bio/AlignIO/__init__.py
in read(handle, format, seq_count, alphabet)
    422     iterator = parse(handle, format, seq_count, alphabet)
    423     try:
--> 424         first = next(iterator)
    425     except StopIteration:
    426         first = None

/Users/cr517/anaconda/lib/python3.6/site-packages/Bio/AlignIO/__init__.py
in parse(handle, format, seq_count, alphabet)
    366                                                 seq_count=seq_count)
    367         else:
--> 368             raise ValueError("Unknown format '%s'" % format)
    369
    370         # This imposes some overhead... wait until we drop Python 2.4 to fix it

ValueError: Unknown format 'mauve'

You're using version 1.68 or older. Mauve support was added in 1.70.


Parse out exon coordinates from bed file for each gene This is my bed file for all the exon coordinates I want to take out all the exons of a given gene. Let say I have gene in chr 1 which starts from chr1 11868   12227 so I want to parse out all the exons that comes in between 11868    12227.

This is my small subset:

cat exon.bed | head -10
chr1    11868   12227   +   exon
chr1    11871   12227   +   exon
chr1    11873   12227   +   exon
chr1    12009   12057   +   exon
chr1    12178   12227   +   exon
chr1    12594   12721   +   exon
chr1    12612   12697   +   exon
chr1    12612   12721   +   exon
chr1    12612   12721   +   exon
chr1    12974   13052   +   exon

How do I parse out?

I use mostly R and bit of shell script but I'm not sure if I can use R, may a few lines of perl or shell script can help me solve my problem.

There are multiple ways to go about it. On the command line you can make a 1 line BED file:

chr1   11868   12227

And then bedtools intersect with it.

In R, you could load your original BED file and use GenomicRanges:

library("GenomicRanges")
bed = read.delim("foo.bed", header=F)  # Rename this
# N.B., BED files use 0-based coordinates, I've switched to 1-based here
exons = GRanges(seqnames=levels(bed$V1),
                ranges=IRanges(start=bed$V2+1, end=bed$V3),
                strand=bed$V4)
# Again, I switched to 1-based coordinates...you don't have to do that
region = GRanges(seqnames="chr1",
                 ranges=IRanges(start=11869, end=12227),
                 strand="+")
theExonsYouWant = subsetByOverlaps(exons, region, type="within")

You'll need to either specify a strand in your region or use ignore.strand=T in subsetByOverlaps(). I hope you don't have any genes that overlap on the same strand, since they'll prove problematic.


BioMart does not return a result I am supplying to BioMart 81 ensemblIDs, but getting 80 results where 1 ensembleID: ENSMUSG00000083840 is absent. What does it mean if BioMart can not find the gene? Does it mean that it is a non-coding one or that there is just no symbol attached to it? Here is my BioMart query:

library('biomaRt')
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
genes <- rownames(mat)
G_list <- getBM(filters= "ensembl_gene_id", attributes=
c("ensembl_gene_id","mgi_symbol","description"),values=genes,mart= mart)

Whenever you are wondering about things like this, just look up the identifier on the Ensembl web page. If you look up ENSMUSG00000083840, you will see:

  This identifier is not in the current EnsEMBL database

In this case, it doesn't seem to have been replaced by anything and it wasn't a protein coding gene, so it looks like it was simply a wrong annotation which has been removed.


To end up with a single table made by merging all your tables, you can give a try with reduce() of purrr:  The first argument .x would be the list of tables to merge (in your case the annotation file plus all the other files) and the .y argument would be your function, i.e. inner_join(). reduce() would merge all the tables within the given list, left to right or right to left. ... would enable you to specify column(s) to use as "anchors" while merging.

If you would prefer to keep your tables separately then lapply() is the answer: lapply(list_of_files, function(x) inner_join(annotation_file, x)). This would give a list of merged tables.


How to get PSSM of protein by protein sequence with python? I didn't find way that convert protein sequences in  FASTA to PSSM by python,
and how to write program in python to find PSSM for our data set ?
Maybe better to question how to find PSSM from NCBI alignments ?

Thank you

We solve our problem and create gitHub repository (bio-protein).


Wrong format .ped in PLINK I know that I don't have the correct format in my .ped file for PLINK but I don't know how to convert it.
My actual file is something like:

BSCAN BSCAN20016 0 0 2 0 0 1 1 1....
BSCAN BSCAN20017 0 0 2 0 0 1 1 1....
BSCAN BSCAN20018 0 0 2 0 0 1 1 1....

Where column 1 is FamilyID, column 2 is ID, column 3 is FatherID, column 4 MotherID, column 5 Sex and then the genotype of each individual.
I only need to create binary files, so I tried to run:

plink-1.9-rc --file admixed --make-bed --no-pheno --cow --out admixed

but the program shows: Error: half-missing call in .ped file at variant 1, line 1.

Anyone know how can I change the format or recode it to make a correct input?

Thanks in advance!

So, columns 7 and 8 should be the genotype of variant 1 (each column is one of the alleles). You have them as 0 1, which means absent and major allele.

I think you're just missing a column. From the doc:

Column 6 is Phenotype value ('1' = control, '2' = case, '-9'/'0'/non-numeric = missing data if case/control)

Ideally you would know which SNPs are heterozygous based on your SNP calling and have such SNPs as 1 2, but if you have a homozygous organism or if you're working with data that's been flattened you can simply generate the .ped file with the same number twice for each SNP.


How can I easily get the read size distribution of reads mapping on a certain set of regions? Suppose I have a BAM file indicating where reads in a library have mapped, and a bed file describing a set of genomic regions.

Is there a way to easily get the size distribution of the reads mapping on this set of regions?

This can be done rather simply by combining this samtools based answer and bioawk:

samtools view -L regions.bed aligned_reads.bam \
    | bioawk -c sam '{hist[length($seq)]++} END {for (l in hist) print l, hist[l]}' \
    | sort -n -k1

The samtools command extracts the alignements in the region given with the -L option, in sam format. The -c sam option of bioawk makes the $seq variable (corresponding to the read sequence field) available in the awk command.

The final sort ensures the size classes are ordered numerically.

We obtain a two-column size    count table.

Examples with timings:

$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | bioawk -c sam '{histo[length($seq)]++} END {for (l in histo) print l, histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.592s
user    0m0.816s
sys     0m0.036s

Since sam format is actually simple to parse, one can try to use mawk instead for slightly higher speed (the sequence being in the 10th field):

samtools view -L regions.bed aligned_reads.bam \
    | mawk '{hist[length($10)]++} END {for (l in hist) print l"\t"hist[l]}' \
    | sort -n -k1

(Contrary to bioawk, the default output field separator is not a tabulation.)

Example:

$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | mawk '{histo[length($10)]++} END {for (l in histo) print l"\t"histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.362s
user    0m0.500s
sys     0m0.004s


Rosetta ab initio run gives 0.00 RMSD for generated 2000 models? (run still in progress) I am running rosetta ab initio prediction using the command line:

"mpiexec -np 16 /main/source/bin/AbinitioRelax.mpi.linuxgccrelease -database /main/database/ @options -mpi_tracer_to_file log > log1.log &"

"Options" file looks like this: https://ibb.co/V9VrQHV

RMSD vs Energy plot is:
https://ibb.co/1mjzwDj

There is no native structure as this is an ab initio prediction. Could you please let me know about your comments on the RMSD being 0.00 for all 2000 models generated till now? Is this normal and will produce results after more runs?

RMSD is a measure of similarity between two structures. The relevant RMSD in protein structure prediction is between a model and the native structure. If you don't have a native structure you cannot have an RMSD in this sense. Ignore the fact it is 0.0, there is possibly some error that means you are calculating the RMSD between a model and itself.

And that's the difficulty of protein structure prediction... you don't know whether you are right or not! The structures you have generated are your results, and you'll need to find ways to filter or use them that are not reliant on RMSD. You could use estimation of model accuracy methods or clustering, for example.


Are fgsea and Broad Institute GSEA equivalent? Several gene set enrichment methods are available, the most famous/popular is the Broad Institute tool. Many other tools are available (See for example the biocView of GSE which list 82 different packages). There are several parameters in consideration :

the statistic used to order the genes,
if it competitive or self-contained,
if it is supervised or not,
and how is the enrichment score calculated.

I am using the fgsea - Fast Gene Set Enrichment Analysis package to calculate the enrichment scores and someone told me that the numbers are different from the ones on the Broad Institute despite all the other parameters being equivalent.

Are these two methods (fgsea and Broad Institute GSEA) equivalent to calculate the enrichment score?

I looked to the algorithms of both papers, and they seem fairly similar, but I don't know if in real datasets they are equivalent or not.

Is there any article reviewing and comparing how does the enrichment score method affect to the result?

According to the FGSEA preprint:

  We ran reference GSEA with default parameters. The permutation number
  was set to 1000, which means that for each input gene set 1000
  independent samples were generated. The run took 100 seconds and
  resulted in 79 gene sets with GSEA-adjusted FDR q-value of less than
  10−2. All significant gene sets were in a positive mode. First, to get
  a similar nominal p-values accuracy we ran FGSEA algorithm on 1000
  permutations. This took 2 seconds, but resulted in no significant hits
  due after multiple testing correction (with FRD ≤ 1%).

Thus, FGSEA and GSEA are not identical.

And again in the conclusion:

  Consequently, gene sets can be ranked more precisely in the results
  and, which is even more important, standard multiple testing
  correction methods can be applied instead of approximate ones as in
  [GSEA].

The author argues that FGSEA is more accurate, so it can't be equivalent.

If you are interested specifically in the enrichment score, that was addressed by the author in the preprint comments:

  Values of enrichment scores and normalized enrichment scores are the
  same for both broad version and fgsea.

So that part seems to be the same.


What is the difference between FASTA, FASTQ, and SAM file formats? I'd like to learn the differences between 3 common formats such as FASTA, FASTQ and SAM. How they are different? Are there any benefits of using one over another?

Based on Wikipedia pages, I can't tell the differences between them.

Let’s start with what they have in common: All three formats store

sequence data, and
sequence metadata.

Furthermore, all three formats are text-based.

However, beyond that all three formats are different and serve different purposes.

Let’s start with the simplest format:

FASTA

FASTA stores a variable number of sequence records, and for each record it stores the sequence itself, and a sequence ID. Each record starts with a header line whose first character is >, followed by the sequence ID. The next lines of a record contain the actual sequence.

The Wikipedia artice gives several examples for peptide sequences, but since FASTQ and SAM are used exclusively (?) for nucleotide sequences, here’s a nucleotide example:

>Mus_musculus_tRNA-Ala-AGC-1-1 (chr13.trna34-AlaAGC)
GGGGGTGTAGCTCAGTGGTAGAGCGCGTGCTTAGCATGCACGAGGcCCTGGGTTCGATCC
CCAGCACCTCCA
>Mus_musculus_tRNA-Ala-AGC-10-1 (chr13.trna457-AlaAGC)
GGGGGATTAGCTCAAATGGTAGAGCGCTCGCTTAGCATGCAAGAGGtAGTGGGATCGATG
CCCACATCCTCCA

The ID can be in any arbitrary format, although several conventions exist.

In the context of nucleotide sequences, FASTA is mostly used to store reference data; that is, data extracted from a curated database; the above is adapted from GtRNAdb (a database of tRNA sequences).

FASTQ

FASTQ was conceived to solve a specific problem arising during sequencing: Due to how different sequencing technologies work, the confidence in each base call (that is, the estimated probability of having correctly identified a given nucleotide) varies. This is expressed in the Phred quality score. FASTA had no standardised way of encoding this. By contrast, a FASTQ record contains a sequence of quality scores for each nucleotide.

A FASTQ record has the following format:

A line starting with @, containing the sequence ID.
One or more lines that contain the sequence.
A new line starting with the character +, and being either empty or repeating the sequence ID.
One or more lines that contain the quality scores.

Here’s an example of a FASTQ file with two records:

@071112_SLXA-EAS1_s_7:5:1:817:345
GGGTGATGGCCGCTGCCGATGGCGTC
AAATCCCACC
+
IIIIIIIIIIIIIIIIIIIIIIIIII
IIII9IG9IC
@071112_SLXA-EAS1_s_7:5:1:801:338
GTTCAGGGATACGACGTTTGTATTTTAAGAATCTGA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII6IBI

FASTQ files are mostly used to store short-read data from high-throughput sequencing experiments. The sequence and quality scores are usually put into a single line each, and indeed many tools assume that each record in a FASTQ file is exactly four lines long, even though this isn’t guaranteed.

As with FASTA, the format of the sequence ID isn’t standardised, but different producers of FASTQ use fixed notations that follow strict conventions.

SAM

SAM files are so complex that a complete description [PDF] takes 15 pages. So here’s the short version.

The original purpose of SAM files is to store mapping information for sequences from high-throughput sequencing. As a consequence, a SAM record needs to store more than just the sequence and its quality, it also needs to store information about where and how a sequence maps into the reference.

Unlike the previous formats, SAM is tab-based, and each record, consisting of either 11 or 12 fields, fills exactly one line. Here’s an example (tabs replaced by fixed-width spacing):

r001  99  chr1  7 30  17M         =  37  39  TTAGATAAAGGATACTG   IIIIIIIIIIIIIIIII
r002  0   chrX  9 30  3S6M1P1I4M  *  0   0   AAAAGATAAGGATA      IIIIIIIIII6IBI    NM:i:1

For a description of the individual fields, refer to the documentation. The relevant bit is this: SAM can express exactly the same information as FASTQ, plus, as mentioned, the mapping information. However, SAM is also used to store read data without mapping information.

In addition to sequence records, SAM files can also contain a header, which stores information about the reference that the sequences were mapped to, and the tool used to create the SAM file. Header information precede the sequence records, and consist of lines starting with @.

SAM itself is almost never used as a storage format; instead, files are stored in BAM format, which is a compact binary representation of SAM. It stores the same information, just more efficiently and, in conjunction with a search index, allows fast retrieval of individual records from the middle of the file (= fast random access). BAM files are also much more compact than compressed FASTQ or FASTA files.

The above implies a hierarchy in what the formats can store: FASTA ⊂ FASTQ ⊂ SAM.

In a typical high-throughput analysis workflow, you will encounter all three file types:

FASTA to store the reference genome/transcriptome that the sequence fragments will be mapped to.
FASTQ to store the sequence fragments before mapping.
SAM/BAM to store the sequence fragments after mapping.


Improve a reference genome with sequencing data I have a DNA sample which I know doesn't quite match my reference genome - my culture comes from a subpopulation which has undergone significant mutation since the reference was created.

From visual inspection with IGV, a significant number of both SNPs and SVs appear to be present, but an assembly built entirely from my own sequencing data is not high enough quality for my purposes.

How can I modify this reference genome to match my sample with new sequencing data (preferably with Oxford Nanopore Technologies long reads, but I can also use these to scaffold short reads if necessary), taking advantage of my knowledge that the existing reference is mostly very good, without having to access the reads which were originally used to construct the reference genome?

One approach to this is to use whatever data you have to iteratively update the reference genome. You can keep chain files along the way so you can convert coordinates (e.g. in gff files) from the original reference to your new pseudoreference.

A simple approach might be:

Align new data to existing reference
Call variants (e.g. samtools mpileup, GATK, or whatever is best for you)
Create new reference incorporating variants from 2
Rinse and repeat (i.e. go to 1)

You can track some simple stats as you do this - e.g. the number of new variants should decrease, the number of reads mapped should increase, and the mismatch rate should decrease, with every iteration of the above loop. Once the pseudoreference stabilises, you know you can't do much more.


creating population AF from set of individuals in 1000 genomes population? What's the recommended way to create population-level Allele Frequency variant entries from a set of individuals (vcfs?) in the 1000 genomes (or similar) project?

For example, given the IBS individuals in the 1000 Genomes project:

http://www.internationalgenome.org/data-portal/population/IBS

I would like to create an IBS .vcf.gz that contains an IBS population-level AF entry for each of the variants present.

Any ideas on how to achieve this?

EDIT: this is related to my earlier question on how to create the most complete file for population allele frequencies, here:
what is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets?

I have seen that there is an Allele Frequency Calculator tool in EnsEMBL, which is at least available with GRCh37:

http://www.internationalgenome.org/allele-frequency-calculator-documentation

http://grch37.ensembl.org/Homo_sapiens/Tools/AlleleFrequency?db=core

There are also instructions on how to do this genome-wide using vcftools:

vcf-subset -c ~/1000g/IBS.samples.list ~/1000g/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz | fill-an-ac | bgzip -c > ~/1000g/IBS.chr1.20130502.vcf.gz

This seems to be generating correct files, and the AF can be calculated with a tool that combines AC and AN values, such as it is described here:

https://www.biostars.org/p/123468/

Any other ideas / warnings / pitfalls of the approach?

bcftools has a plugin called fill-tags which should do the job for you easily, without having to mess around with lots of complicated bash scripting.

You could do something along the lines of this, presuming you have a .vcf file with all the IBS individuals in it:

bcftools +fill-tags in.vcf.gz -Oz -o out.vcf.gz

This will give each position in the .vcf file new INFO tags, one of which will be allele-frequency.

You can then pull off the frequency for a particular variant (CHR:POS) by:

bcftools query -t CHR:POS -f'%CHROM\t%POS\t%AF\n' out.vcf.gz


Why does the FASTA sequence for coronavirus look like DNA, not RNA? I'm looking at a genome sequence for 2019-nCoV on NCBI. The FASTA sequence looks like this:

>MN988713.1 Wuhan seafood market pneumonia virus isolate 2019-nCoV/USA-IL1/2020, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC
TAATTACTGTCGTTGACAGGACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTG
...
...
TTAATCAGTGTGTAACATTAGGGAGGACTTGAAAGAGCCACCACATTTTCACCGAGGCCACGCGGAGTAC
GATCGAGTGTACAGTGAACAATGCTAGGGAGAGCTGCCTATATGGAAGAGCCCTAATGTGTAAAATTAAT
TTTAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGGAGAATGACAAAAAAAAAAAA

Coronavirus is an RNA virus, so I was expecting the sequence to consist of AUGC characters. But the letters here are ATGC, which looks like DNA!

I found a possible answer, that this is the sequence of a "complementary DNA". I read that

  The term cDNA is also used, typically in a bioinformatics context, to refer to an mRNA transcript's sequence, expressed as DNA bases (GCAT) rather than RNA bases (GCAU).

However, I don't believe this theory that I'm looking at a cDNA. If this were true, the end of the true mRNA sequence would be ...UCUUACUGUUUUUUUUUUUU, or a "poly(U)" tail. But I believe the coronavirus has a poly(A) tail.

I also found that the start of all highlighted genes begin with the sequence ATG. This is the DNA equivalent of the RNA start codon AUG.

So, I believe what I'm looking at is the true mRNA, in 5'→3' direction, but with all U converted to T.

So, is this really what I'm looking at? Is this some formatting/representation issue? Or does 2019-nCoV really contain DNA, rather than RNA?

That is the correct sequence for 2019-nCov. Coronavirus is of course an RNA virus and in fact, to my knowledge, every RNA virus in Genbank is present as cDNA (AGCT, i.e. thydmine) and not RNA (AGCU, i.e. uracil).

The reason is simple, we never sequence directly from RNA because RNA is too unstable and easily degraded by RNase. Instead the genome is reverse transcribed, either by targeted reverse transcription or random amplification and thus converted to cDNA. cDNA is stable and is essentially reverse transcribed RNA.

The cDNA is either sequenced directly or further amplified by PCR and then sequenced. Hence the sequence we observe is the cDNA rather than RNA, thus we observe thymine rather than uracil and that is how it is reported.


How to calculate the memory usage of storing kmers in RAM I want to write a program in C++ that stores kmers in a hash or in a trie. How can I calculate how much RAM I would need for each type of data structure? For this application my kmers are strand-specific, so I cannot reduce the size complexity by considering canonical kmers.

For a hash, for example, I know that there are $4^k$ possible kmers of length $k$. I also am assuming that DNA character with an alphabet $\in \{G,A,T,C\}$ can be stored in two bits ($2^2 = 4$). Using this logic the kmer hash should occupy $(4^k) \cdot 2 \cdot k$ bits. Plugging this into Wolfram Alpha ((4^13) * 2 * 13 bits to gigabytes), I calculated the amount of RAM that some small kmer sizes would require to store the entire hash in memory in GB.

k = 13 requires    0.2181 GB
k = 15 requires    4.027  GB
k = 17 requires   73.01  GB
k = 19 requires 1306   GB

This escalated rather quickly, and I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM?

The kmers that I want to store are sometimes sparse compared to the entire space (for k=13, only 29 million or so out of 67 million possibilities).

Constraints

I do not need to store the number of occurences of a kmer, just its existence
I cannot risk false positives. Bloom filters won't work in this case. :^/

To follow up on Devon Ryan's answer, I thought it would be a little fun to write a Python script that demonstrates using a bit array to maintain a presence/absence table.

Note: I wrote a C++ port that includes a custom bitset implementation that can be sized at runtime. This and the Python script are available on Github: https://github.com/alexpreynolds/kmer-boolean

#!/usr/bin/env python

import sys
import os
import bitarray

# read FASTA
def read_sequences():
    global seqs
    seqs = []
    seq = ""
    for line in sys.stdin:
        if line.startswith('>'):
            if len(seq) > 0:
                seqs.append(seq)
                seq = ""
        else:
            seq += line.strip()
    seqs.append(seq)

# build and initialize bit array
def initialize_bitarray():
    global ba
    ba = bitarray.bitarray(4**k)
    ba.setall(False)
    sys.stderr.write("Memory usage of bitarray.bitarray instance is [%ld] bytes\n" % (ba.buffer_info()[1]))

# process sequences
def process_sequences():
    global observed_kmers
    observed_kmers = {}
    for seq in seqs:
        for i in range(0, len(seq)):
            kmer = seq[i:i+k]
            if len(kmer) == k:
                observed_kmers[kmer] = None
                idx = 0
                for j in range(k-1, -1, -1):
                    idx += 4**(k-j-1) * bm[kmer[j]]
                ba[idx] = True

def test_bitarray():
    test_idx = 0
    for j in range(k-1, -1, -1):
        test_idx += 4**(k-j-1) * bm[test_kmer[j]]
    test_result = ba[test_idx]
    if test_result:
        sys.stdout.write("%s found\n" % (test_kmer))
        sys.exit(os.EX_OK)
    else:
        sys.stdout.write("%s not found\n" % (test_kmer))
        sys.exit(os.EX_DATAERR)

def main():
    global k
    k = int(sys.argv[1])

    global bm
    bm = { 'A' : 0, 'C' : 1, 'T' : 2, 'G' : 3 }

    read_sequences()
    initialize_bitarray()
    process_sequences()

    try:
        global test_kmer
        test_kmer = sys.argv[2]
        if len(test_kmer) == k:
            test_bitarray()
        else:
            raise ValueError("test kmer (%s) should be of length k (%d)" % (test_kmer, k))
    except IndexError as err:
        keys = list(observed_kmers.keys())
        for i in range(0, len(keys)):
            sys.stdout.write("%s found\n" % (keys[i]))

    sys.exit(os.EX_OK)

if __name__== "__main__":
    main()

Note that this doesn't look at canonical kmers, e.g., AG is considered a distinct 2mer from its reverse complement CT.

To use this script, you pipe in your FASTA, specify the k, and an optional kmer that you want to test for presence/absence, e.g.:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 2 CG
Memory usage of bitarray.bitarray instance is [2] bytes
CG found

Or:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 3 AAA
Memory usage of bitarray.bitarray instance is [8] bytes
AAA not found

Or if the optional test kmer is left out:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 5
Memory usage of bitarray.bitarray instance is [128] bytes
CATTC found
ATTCT found
TTCTC found
TCTCG found
CTCGG found
TCGGG found
CGGGA found
GGGAC found
TTATA found
TATAT found
TTTAT found
TTATT found
TATTA found
ATTAG found
TTAGA found
TAGAC found
AGACC found
GACCT found
ACCTC found
CCTCT found

Or for the ~67M kmers in a 13-mer set, for which a roughly 8.4MB bit array is reserved:

$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 13
Memory usage of bitarray.bitarray instance is [8388608] bytes
TTTATTAGACCTC found
TTATTAGACCTCT found


TRAL does not find "Phobos result file" I want to use TRAL to annotate tandem repeats in the reference genome of Caenorhabditis elegans. For this, I need to install some external software, such as Phobos. I've downloaded Phobos and I am using it by typing the path to it's executable like this:

input/software/phobos-v3.3.12-mac/bin/phobos-mac-intel-64bit-start-from-command-line genome.fa output.txt

This works quite well.

In order for TRAL to find Phobos, I changed TRAL's configuration file as explained here.

I've replaced PHOBOS = phobos with PHOBOS = /Users/user/Documents/.../project/input/software/phobos-v3.3.12-mac/bin/phobos-mac-intel-64bit-start-from-command-line

Then, I followed the tutorial Run and parse de novo repeat detection software and typed the following commands in Python 3 in the Terminal:

import os
from tral.sequence import sequence
from tral.paths import PACKAGE_DIRECTORY

proteome_HIV = os.path.join(PACKAGE_DIRECTORY, "examples", "data", "HIV-1_388796.faa")
sequences_HIV = sequence.Sequence.create(file = proteome_HIV, input_format = 'fasta')

So far so good. But the following command creates an error:

tandem_repeats = sequences_HIV[0].detect(denovo = True, detection = {"detectors": ["Phobos"]})

Did not find Phobos result file in /var/folders/z2/rt6qxntd3tg5bjmqlqnm0mdm0000gn/T/tmpvc20y3i4/001/PHOBOS/phobos.o

Note that if I use PHOBOS instead of Phobos in the command, I get the following error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/sequence.py", line 162, in detect
    **kwargs['detection'])[0]
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/repeat_detection_run.py", line 962, in run_detector
    Detectors(detectors, sequence_type)
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/tral/sequence/repeat_detection_run.py", line 872, in Detectors
    FINDER_LIST.keys()))
Exception: Unknown TR detector supplied (Supplied: ['PHOBOS']. Known TR detectors: dict_keys(['HHrepID', 'Phobos', 'TRED', 'T-REKS', 'TRF', 'TRUST', 'XSTREAM']))

I am working on Mac OS X. I installed TRAL with pip install tral.

TRAL 0.3.5 has bugs running and parsing Phobos. You can follow progress towards a fix on this issue: https://github.com/acg-team/tral/issues/16

Until the issue is fixed, you should remove Phobos from the repeat detector list in your tral config file (~/.tral/config.ini) and use a different repeat detector.

[sequence]
    [[repeat_detection]]
        DNA = TRED, T-REKS, TRF, XSTREAM


Heteroplasmy Calculation For calculating heteroplasmy in mtDNA I was using MitoSeek, the tool provides position wise heteroplasmy percentage. So, will it be wise to add up all sites to get a position-wise percentage in a single sample, or calculate it by adding up all counts (variants) from an aligned file and divide it by total read count?

Heteroplasmy at any site (as long as it's a reliable variant) suggests that there is mitochondrial heteroplasmy, it doesn't matter where it is. In that case, the "mitochondrial heteroplasmy" statistic would simply be the maximum value at any [reliable] location.

Due to the low amount of variation in the mitochondrial genome, and because it's a haploid genome, I don't think that a calculation like "average heteroplasmy" makes much sense for mitochondria.


RNAseq DE comparison with samples of different read length I have RNA samples but with different read lengths (Eg, HiSeq 2x125 and NovaSeq 2x150bp data). I would like to do DE analysis on these samples. What do you recommend? Do you recommend to trim the reads to 125bp for the NovaSeq to compare with the other batches from HiSeq? Or any other suggestions.

You don't need to do anything, this is taken care of by the linear model - if you don't have any confounding effects.
If you do have confounding effects, e.g. all control samples on HiSeq and all treated on NovaSeq there's not a lot you can do anyway, other than be honest when you present your results.
Taken care of by the model means adding an extra 2 level factor for sequencer by sample, you need to do this anyway to take care of the batch effect between sequencer runs.


Multiple input & output files? I'm trying to print out the measured distances for each .pdbqt file, and then write a new output file for each .pdbqt, using the original name of the .pdbqt file as the new .txt file name. Below is the script.
I'm using python3 & I'm trying to do this via Linux CLI.
from pymol import stored
stored.a, stored.b, stored.c, stored.d = [],[],[],[]
cmd.iterate_state(1, "r55_01///UNL`1/N*", "stored.a.append(ID)")
cmd.iterate_state(1, "r55_01///UNL`1/O*", "stored.b.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASP`113/OD2", "stored.c.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASN`312/ND2", "stored.d.append(ID)")

f=open("tt.txt", 'w')

for a in stored.a:
  for c in stored.c:
      d1=cmd.distance("id %s" % a, "id %s" % c,cutoff=5)
      f.write("%s D113 %s\n" % (a,d1))
      print(d1)

for b in stored.b:
  for d in stored.d:
      d2 = cmd.distance("id %s" % b, "id %s" % d, cutoff=5)
      f.write("%s N312 %s\n" % (b,d2))
      print(d2)

f.close()

My folder looks something like this
.
├── dst6.txt
├── sp
│   ├── r55_01.pdbqt
│   ├── r55_02.pdbqt
│   ├── r55_03.pdbqt
│   ├── r55_04.pdbqt
│   ├── r55_05.pdbqt
│   ├── r55_06.pdbqt
│   ├── r55_07.pdbqt
│   ├── r55_08.pdbqt
│   ├── r55_09.pdbqt
│   ├── r55_10.pdbqt
│   ├── r56_01.pdbqt
│   ├── r56_02.pdbqt
│   ├── r56_03.pdbqt
│   ├── r56_04.pdbqt
│   ├── r56_05.pdbqt
│   ├── r56_06.pdbqt
│   ├── r56_07.pdbqt
│   ├── r56_08.pdbqt
│   ├── r56_09.pdbqt
│   └── r56_10.pdbqt
└── t6.py

THANKS IN ADVANCE!!

You made the hardest work, you were not far to achieve what you wanted ;)
This script.py should work after you made it executable (in a terminal type: chmod u+x script.py), just type ./script.py in your terminal in your working directory (so here you should be in your "sp" directory). Note that I got inspired by: https://www.biostars.org/p/42474/ to know how to launch pymol within the script and also that I would personally not use pymol for doing this kind of analysis/data extraction from a pdb/pdbqt file but rather use a dedicated analysis package (there are plenty like MDanalysis). I personally use pymol (pymol script) for visualization purpose.
#!/usr/bin/env python

import __main__
__main__.pymol_argv = ['pymol', '-qc']  # Quiet and no GUI

from pymol import stored
import pymol, os

pymol.finish_launching()

#put what you've done into a function
def myfunc(file):
    spath = os.path.abspath(file)
    sname = spath.split('/')[-1].split('.')[0]

    pymol.cmd.load(spath, sname)

    stored.a, stored.b, stored.c, stored.d = [], [], [], []
    pymol.cmd.iterate_state(1, "%s///UNL`1/N*" % sname, "stored.a.append(ID)")
    pymol.cmd.iterate_state(1, "%s///UNL`1/O*" % sname, "stored.b.append(ID)")
    pymol.cmd.iterate_state(1, "%s//R/ASP`113/OD2" % sname, "stored.c.append(ID)")
    pymol.cmd.iterate_state(1, "%s//R/ASN`312/ND2" % sname, "stored.d.append(ID)")

    f = open('%s.txt' % sname, 'w')

    for a in stored.a:
        for c in stored.c:
            d1 = pymol.cmd.distance("id %s" % a, "id %s" % c, cutoff=5)
            f.write("%s D113 %s\n" % (a, d1))
            print(d1)

    for b in stored.b:
        for d in stored.d:
            d2 = pymol.cmd.distance("id %s" % b, "id %s" % d, cutoff=5)
            f.write("%s N312 %s\n" % (b, d2))
            print(d2)

    f.close()

mypdbqt = os.popen('ls *.pdbqt').read().split()

for i in mypdbqt:
    myfunc(i)

pymol.cmd.quit()


How to subset a BAM by a list of QNAMEs? I have a text file 'qnames.txt' with QNAMEs in the following format:
EXAMPLE:QNAME1
EXAMPLE:QNAME2
EXAMPLE:QNAME3
EXAMPLE:QNAME4
EXAMPLE:QNAME5

I would like to subset my BAM file.bam via all of these QNAMEs into a new SAM.
Naturally, I can do this individually, e.g.
samtools view file.bam | grep 'EXAMPLE:QNAME1' > subset.bam

But I'm unsure how to do this for a list of QNAMES:

How do I write a for loop that will do all of these queries, outputting the correct SAM needed?
I could write a for loop that creates n SAM files, and then cats them…

Is there a way to specifically grep by QNAME? The above may grep reads that may not be associated with the correct QNAME?

How do I keep the BAM header?

use picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads

  READ_LIST_FILE (File) Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null.


How to concatenate "by chromosome"-VCFs? I have a several VCFs which are VCF which only contain information by chromosome. That is, there's a chromosome 1 VCF (with only chr1), a chromosome 2 VCF (with only chr2), etc.

I checked to make sure that these VCFs were valid via VCFtools, i.e.

$ vcf-validator chr1.vcf

which works---these are valid VCFs I was given.

Now, I would like to combine these VCFs into one VCF.

I naïvely tried the following cat operation:

$ cat chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf > total_chroms.vcf

This doesn't work properly though

$ vcf-validator total_chroms.vcf
The header tag 'contig' not present for CHROM=chr1. (Not required but highly recommended.)
Could not parse the line, wrong number of columns: [##fileformat=VCFv4.2\n]
 at /path/vcftools-0.1.14/perl/Vcf.pm line 172, <__ANONIO__> line 191016.
     Vcf::throw('Vcf4_2=HASH(0x1ae7208)', 'Could not parse the line, wrong number of columns: [##filefor...') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 335
     Vcf::next_data_array('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 3457
     Vcf4_1::next_data_array('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 2574
     VcfReader::run_validation('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14//bin/vcf-validator line 60
     main::do_validation('HASH(0x16ada68)') called at  /path/vcftools-0.1.14//bin/vcf-validator line 14
$

What tools are available to merge these VCF's together into a total VCF?

I would recommend bcftools concat.  You can't just cat them together because each file has a header section.  The bcftools command will handle all that for you.  Each vcf file must be sorted prior to calling concat

bcftools concat -o total_chroms.vcf chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf


Why is 42 the seed used in RunPCA() in Seurat? The RunPCA() in Seurat 3.0 function has seed.use=42 by default. How does one choose a seed number? Is this chosen simply because it is the Answer to the Ultimate Question of Life, the Universe, and Everything? Or maybe Lewis Carroll's curious obsession with the number?

Given that it is virtually impossible for a human to predict what random numbers will be generated given a certain seed, and what effects they will have for a given application, the choice of the seed can be considered arbitrary.

What matters is that fixing the seed should fix the behaviour of the application. What seed has been chosen does not really matter, as long as it is documented.

The choice of 42 could therefore well come from one of the reasons you propose. Only the authors of the programs could confirm.


How to use variables while calling pybedtools I am trying to use bedtools and fetch sequences from a whole genome fasta file inside the script I get region coordinates as variables.
For example:

chr_ = "chr1"
start = 3000
end = 3402

I am not sure how to wrap my function, because the example provided by BedTool.sequence only uses coordinates in a string form.

Can you please help me to have something like this?

a = pybedtools.BedTool(chr_, start, end, from_variables=True)

You have to take your variables and make them into a single string variable.

for example:

chr_ = "chr1"
start = 3000
end = 3402
as_str = ' '.join([chr_, str(start), str(end)])
# or
as_str2 = "{} {} {}".format(chr_, start, end)

Once you have the string built you can use the example in the bedTool docs.

a = BedTool(as_str, from_string=True)


Easy way to create gene plot in gviz? I would like to create a simple plot of a given human gene containing the main features, e.g. exons, introns, ORFs, transcripts etc. I understand that gviz  is a popular choice for this kind of thing, but most tutorials are massive overkill for what I want to do (I don't need fancy features such as importing BAM files and plotting coverage).
I was wondering whether someone can point me towards some sample code for creating a basic gene plot. It probably would be good to start by creating a dataframe containing all the features of my gene - is there a handy package that can do that?

If Python is an option, you could check BioPython and GenomDiagram. This should get the job done.


Should one be especially conscious of removing low counts, when applying batch correction? When (as a mistake) I did not remove low counts at all (beside those that equal zero for all samples), I got the following ma plot (using Glimma) :

On  the right you see the individual counts. The gene that is displayed is signficant. However, the reality seems that the gene is not expressed both in treatment and in control. The batch I have corrected for made it a significantly DE however.
The batch correction algorithm of DESeq2 doesn't care whether the counts of the gene are almost all zeroes, and does not take this as a special case.
In reality, it seems to me that it is a special case - if it is questionable whether the gene is expressed in the first place, it makes no sense to apply batch correction to zero counts.
The solution to this seems to be to make sure to remove low counts when applying batch correction. Usually it's only a recommendation; but when applying batch, it seems a necessity.
Would you agree?
I even wonder if one should remove counts that are low per-condition. Suppose that a gene is not expressed in treatment (and only in treatment), does it really make sense to apply batch correction to this gene?

I don't think this has anything to do with batch correction.
The averages of that gene are distinctly different between groups, so it's correct that the gene is statistically significant.  But with one sample driving all the difference, it's likely to be an artifact, or otherwise not fruitful to follow up on.
As geek_y suggested, a filter like "omit genes unless 2 samples have > 10 counts" would drop that gene.


How do I obtain a vcf read from a pyvcf.Reader? I wondered how would I be able to obtain variant based on chrom and location range.
Optimally, this function would be like:
import vcf
vcf_reader = vcf.Reader(filename=some_file_name)
vcf_read = vcf.variants_between(CHROM, START, END)

Does anyone how to do so?

I think vcf_reader.fetch is what you want.  You have to index the file first though.
If starting from regular text VCF you can run:
bgzip example.vcf
tabix example.vcf.gz

Then the Python code is almost the same as yours:
import vcf
vcf_reader = vcf.Reader(filename="example.vcf.gz")
vcf_read = vcf_reader.fetch("NC_023663.1", 7803980, 7804000)
for record in vcf_read:
    print(record)

For this example VCF you get:
Record(CHROM=NC_023663.1, POS=7803981, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803982, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803983, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803984, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803985, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803986, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803987, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803988, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803989, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803990, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803991, REF=A, ALT=[G, <*>])
Record(CHROM=NC_023663.1, POS=7803992, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803993, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803994, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803995, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803996, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803997, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803998, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803999, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7804000, REF=T, ALT=[<*>])

API reference:
https://pyvcf.readthedocs.io/en/latest/API.html#vcf.Reader.fetch
Note the "zero-based, half-open coordinate system" it describes (which is why 7803980 in the second Python argument becomes POS=7803981 from the VCF but then it stops at POS=7804000).


How can I export a full alignment from IGV as an image? The IGV browser lets you export an alignment as an image (File => Save Image). However, this image only contains those reads that fit in the viewing window:

As you can see in the image above by looking at the scroll bar on the right, this is only a small subset of aligned reads. Is there any way to export an image of the aligned reads that includes all reads and not only those currently shown in the window?

Running an IGV batch script helps increase the number of reads per image, compared with File | Save image, but it does not solve the issue by itself.

To run a script, first make a script file like so with 1 command per line:

snapshotDirectory ./
goto chr2:28582094-28582790
maxPanelHeight 10000
snapshot chr2_28582094_28582790.png

Run the script using Tools | Run batch script menu. The importat trick is the maxPanelHeight 10000 which is what makes this include more reads than the default snapshot.

You can also convert this into a little shell script:

#!/bin/bash
mkdir -p /tmp/igv
sleep 2
echo "snapshotDirectory /tmp/igv"
sleep 2
goto chr2:28582094-28582790
sleep 2
echo "maxPanelHeight 10000"
echo "snapshot chr2_28582094_28582790.png"

Then, just open IGV, load the relevant BAM file and run the script like this to connect to IGV via telnet (credit to Pierre Lindenbaum):

sh foo.sh |  telnet 127.0.0.1 60151


Can we use GTEx data as control data for TCGA data? I am using Recount2 TCGA data and was wondering is it right to use GTEx data as control data for this. I would really appreciate your views on this?

Definately not directly. Probably not even if used indirectly. Different projects, done at different times, with different protocols, different machines, different library preps etc.

We have tested combining GTEx and TCGA data for several different purposes where one might expect the analyses' internal controls would have allowed compatibilty, but even then we have always found a persistent batch effect which we were unable to remove.


R not being able to find the file despite being in correct directory I have question regarding a project I have to do about a health student study. I cannot import the csv file via read.csv or read.csv2 despite having transformed all files to csv and being in the correct folder.
Here are the command lines:
setwd ("C:/Users/fek/Desktop/Analyse i-Share")
getwd()
TA<-read.csv2(file=Biobanque.csv, header=TRUE, sep=",", dec=",", fill=TRUE)

And here are the results in the console window:
setwd ("C:/Users/fek/Desktop/Analyse i-Share")

getwd()

[1] "C:/Users/fek/Desktop/Analyse i-Share"
TA<-read.csv2(file=Biobanque.csv, header=TRUE, sep=",", dec=",", fill=TRUE)

[1] Error in read.table(file = file, header = header, sep = sep, quote = quote,  :
objet 'Biobanque.csv' introuvable
To translate from french, it says it can't find the file!

The problem is just what it says.  By omitting the quotes around "Biobanque.csv", R thinks you are talking about an object whose name is Biobanque.csv. And R is confused because of course you have not created an object by that name.


Number of edges in Graphical Analysis using Jupyter with Go Programming language I'm trying to learn graphical analysis in a biological setting.
The online course suggested using the Go Programming language, so I'm new to it and help to understand why the number of edges that I read from the following code is different than from the site I'm downloading it from.
import (
    "fmt"
    "log"

    "github.com/kortschak/graphprac"
)

infile := "YeastL.dot"
g, err := graphprac.NewGraph(infile)
if err != nil {
    log.Fatalf("failed to read graph: %v", err)
}
{
    fmt.Printf("G has %d nodes and %d edges.\n", g.Nodes().Len(), g.Edges().Len())
}

The output I'm getting is: G has 2361 nodes and 6646 edges.
But from the website http://vlado.fmf.uni-lj.si/pub/networks/data/bio/Yeast/Yeast.htm it's suppose to read 7182 edges.
Thanks!

the dot file in the Git repo has 6646 edges, the pajek file in your link has 7182.
Chances are the latter has been updated since.


How to use Python to count k-mers? I have some FASTQ sequence files and a FASTA file for some regions I'm interested in.

I would like:

Build an index for the FASTA file
Use the index to count number of k-mers occurred in my sequence files

I know how to do this in many k-mer counting tools, such as Kallisto, Jellyfish, DSK, Salmon. However, I'm more interested in Python (khmer is a Python package, all the other k-mer counting tools are C++).

Unfortunately, I can't figure out how to do that in khmer from the documentation:

  http://khmer.readthedocs.io/en/v2.1.1/introduction.html

It talks about RNA-Seq, metagenomics, de Bruijn graph etc... All I need are just counts. Khmer is a popular k-mer tool, so it should be able to do it, but how?

Q: How to build an index and count k-mers using khmer?

EDIT:

I don't need number of unique k-mers, I'm looking for a solution like Jellyfish where we build an index and count k-mers to the index.

Something like:

MY_SEQ_1  100 counts
MY_SEQ_2  134 counts

I wrote a command-line k-mer counter called kmer-counter that will output results in a form that your Python script can consume: https://github.com/alexpreynolds/kmer-counter

You can grab, build and install it like so:

$ git clone https://github.com/alexpreynolds/kmer-counter.git
$ cd kmer-counter
$ make
$ cp kmer-counter /usr/local/bin

Once the binary is in your path, you might use it in Python like so:

k = 6
fastaFile = '/path/to/some/seqs.fa'
kmerCmd = 'kmer-counter --fasta --k=%d %s' % (k, fastaFile)
try:
    output = subprocess.check_output(kmerCmd, shell=True)
    result = {}
    for line in output.splitlines():
        (header, counts) = line.strip().split('\t')
        header = header[1:]
        kmers = dict((k,int(v)) for (k,v) in [d.split(':') for d in counts.split(' ')])
        result[header] = kmers
    sys.stdout.write("%s" % (str(result)))
except subprocess.CalledProcessError as error:
    sys.stderr.write("%s" % (str(error)))

Given example FASTA like this:

>foo
TTAACG
>bar
GTGGAAGTTCTTAGGGCATGGCAAAGAGTCAGAATTTGAC

For k=6, you would get an iterable Python dictionary like this:

{'foo': {'TTAACG': 1, 'CGTTAA': 1}, 'bar': {'GTTCTT': 1, 'AGAACT': 1, 'GAGTCA': 1, 'ATGGCA': 1, 'GAACTT': 1, 'ATTCTG': 1, 'CTAAGA': 1, 'CTTCCA': 1, 'ATTTGA': 1, 'GGAAGT': 1, 'AGGGCA': 1, 'CCTAAG': 1, 'CTCTTT': 1, 'AATTTG': 1, 'TCTGAC': 1, 'TTTGCC': 1, 'CTTAGG': 1, 'TTTGAC': 1, 'GAAGTT': 1, 'CCCTAA': 1, 'AGAATT': 1, 'AGTCAG': 1, 'CTGACT': 1, 'TCTTAG': 1, 'CGTTAA': 1, 'GTGGAA': 1, 'TGCCAT': 1, 'ACTCTT': 1, 'GGGCAT': 1, 'TTAGGG': 1, 'CTTTGC': 1, 'TGGAAG': 1, 'GACTCT': 1, 'CATGCC': 1, 'GCAAAG': 1, 'AAATTC': 1, 'GTCAAA': 1, 'TGACTC': 1, 'TAGGGC': 1, 'AAGTTC': 1, 'ATGCCC': 1, 'TCAAAT': 1, 'CAAAGA': 1, 'AACTTC': 1, 'GTCAGA': 1, 'CAAATT': 1, 'TAAGAA': 1, 'CATGGC': 1, 'AAGAAC': 1, 'AAGAGT': 1, 'TCTTTG': 1, 'TTCCAC': 1, 'TGGCAA': 1, 'GGCAAA': 1, 'AGTTCT': 1, 'AGAGTC': 1, 'TCAGAA': 1, 'GAATTT': 1, 'AAAGAG': 1, 'TGCCCT': 1, 'CCATGC': 1, 'GGCATG': 1, 'TTGCCA': 1, 'CAGAAT': 1, 'AATTCT': 1, 'GCATGG': 1, 'ACTTCC': 1, 'TTCTTA': 1, 'GCCATG': 1, 'GCCCTA': 1, 'TTCTGA': 1}}

You can use standard Python calls to manipulate this dictionary object and get sums of counts per record, for sequence, etc. which seems to answer your question. Please feel free to clarify what you're looking for if this object representation is not clear.

For a fully-fleshed out demonstration, see: https://github.com/alexpreynolds/kmer-counter/blob/master/test/kmer-test.py


Mapping bisulfite reads to a reduced size genome Is it worth mapping bisulfite reads (WGBS) to a reduced size genome?

I'm interested only in modifications in CpG context, thus instead of mapping to a whole genome (human genome) I would:

Extract DNA sequence around each CpG (-1000/+1000 bp)
Map reads (using bismark) to this subset

My idea is to reduce mapping time, however afraid that there might be some mapping bias or other issues.

Don't restrict the portion of the genome you're mapping to, you'll just increase false-positive alignments by doing so. If you're worried about mapping time, then (A) use more cores or (B) use a different mapper (BSMAP was the quickest out there last I looked, though in my opinion it's fine to wait a bit longer for better quality results (I normally suggest bwameth to people)).

As pointed out by user 719016 in the comments, the above does not apply to amplicon sequencing. With that, you'll get VASTLY better performance if align to a reduced-size genome...and without an appreciably increased false-positive rate (unless your target-enrichment for sequencing is poor).


Joining scaffolds to get full genome? I am using Ray to assemble a about 7.2 kb picornavirus genome. I am using k-mer of 55. Following is the exact code that I used.

mpiexec -n 1 Ray -k 55 -s 2-CA-SVA-2_S2_L001_R1_001.fastq_q30_trim_LR_15.fasta  -o rayd55tr

This worked quite well as I got a scaffolds.fasta file with two scaffolds one with 5730 bp and another with 1548 bp.

Now, I am confused how I can join those scaffolds to get one linear genome. I joined them manually and blasted just to check and I got something like below :.

Is there any tool to join these two contigs together to have one complete linear genome? I have reference genome too if needed.

Ragout worked well for me. It is a reference-based scaffolder.

ragout.py --outdir output/ input/ragout.rcp > message.txt 2> error.txt --threads 8 --overwrite

The ragout.rcp file contains a 'recipe' with the reference sequence and the contigs in FASTA format.:

.references = ensembl
.target = assembly

ensembl.fasta = ref.fa
assembly.fasta = contigs.fa

In my case, it reduced the number of contigs 10-fold and increased the N50 10-fold. It added N's between the contigs.


Error in as.vector(x) : no method for coercing this S4 class to a vector I tried to run the following code in R studio. Everything worked fine, except at the last step [write.table(mdat, "recount_mdat.csv")] when I tried to export the 'mdat', I got the following error:

Error in as.vector(x) : no method for coercing this S4 class to a vector

My code:

 library('dplyr')
 library('recount')
 library('magrittr')
 library('ffpe')
 library('RSkittleBrewer')
 library('SummarizedExperiment')
 library('devtools')
 trop <- RSkittleBrewer::RSkittleBrewer('tropical')
 colon_proj <- c('SRP025982', 'SRP42161')
 if(any(!file.exists(file.path(colon_proj, 'rse_gene.Rdata')))) {
 sapply(colon_proj, download_study)
 }

 blood_proj <- c('SRP066834', 'SRP041736', 'SRP060416')
 if(any(!file.exists(file.path(blood_proj, 'rse_gene.Rdata')))) {
 sapply(blood_proj, download_study)
 }
 proj <- c(colon_proj,blood_proj)

 dat <- lapply(proj, function(x) {
 load(file.path(x, 'rse_gene.Rdata'))
 return(rse_gene)
 })
 proj

 sapply(dat, dim)

 metadata <- all_metadata('sra')

 write.table(metadata, "recount_metadata.csv")

 mdat <- do.call(cbind, dat)

 write.table(mdat, "recount_mdat.csv")

is(mdat)
## "RangedSummarizedExperiment" "SummarizedExperiment" "Vector" "Annotated"

You can't cbind a bunch of obscure object types.

If you want merged count tables you should do this:

mdat <- do.call(cbind,lapply(dat,assay))

Where row.names are Ensembl gene IDs and col.names are the SRR accessions.

Then run your table writing command.

-

If you want the coordinates of your genes then do this to make a bed with the genomic locations:

library(rtracklayer)
export.bed(rowRanges(dat[[1]]),'gene_coords.bed')


You can use the same approach as for two fields, just add the third:

awk 'NR==FNR{a[$1,$2,$3]=$4;next} ($1,$2,$3) in a{print $0, a[$1,$2,$3]}' ret lcls


Testing for viral/bacterial infection FASTQ files What's the best way to test if a pair of Illumina FASTQ files from blood DNA extraction contain any reads from viral/bacterial infection from DNA in human blood? Thanks in advance.

EDIT:

Looking at the available options, I started testing Kraken for samples where I have seen an indication of viral reads, I will update with results / further questions.

If you're talking RNA, I have a blog post from a while back about Diamond/MEGAN:

http://achri.blogspot.com/2017/04/data-for-free-viral-metagenomics-in.html

In principle this will work for DNA if you have a non-trivial number of pathogen reads (since they tend to be gene dense in their genomes). Otherwise I use Kraken on a database for all the NCBI nt database, but this is slower and doesn't let you sanity check the false positive taxonomy assignments as well as I describe in the post. Of course, ideally, you run multiple methods if you have the resources to do so.


Database of RNA-Seq and SNP genotyping data I need to obtain publicly available RNA-seq + SNP genotyping data for a cohort - ideally one where the RNA-seq data has been sequenced from brain tissue.

I have found CMC and GTEx to be sufficient for my needs, however the procedure to obtaining data from both these sources are extremely long and they request us to share our data with them as well; which we can not do because it is not ours to share.

I was wondering if anyone knows of another resource that can be used for obtaining such data?

Thanks

Human RNA-seq data, like genetic data is very sensitive data as it can be used to identify people (and their relatives), so people will likely not distribute it without entering a formal agreement that ensures that people can be held accountable if the data is misused, stolen or made public. Depending on the consent form signed by the patients (or their relatives), they may not be allowed to share it at all. One possibility would be to start a collaboration with a group that has access to such data and ask them to run your analyses on their data. You can then have a shared publication and it's a win for everybody.


Dealing with duplicate keys I got the following error:

attempting to use to_dict to index a multiple sequence FASTA file with duplicate keys.

Are there methods for working around duplicate keys in biopython? I was thinking to simply remove the duplicate entries, but I would need to delete the entire entry associated with the keys.

Yes, Biopython deliberately does not allow duplicate keys with .to_dict(), e.g.
Input file:
$ cat test.fa
>1
ATG
>2
AAA
>1
CG

Code:
from Bio import SeqIO

print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta")))

Traceback:
  File "script.py", line 3, in <module>
    print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta")))
  File "/path/lib/python3.7/site-packages/Bio/SeqIO/__init__.py", line 789, in to_dict
    raise ValueError("Duplicate key '%s'" % key)
ValueError: Duplicate key '1'

However, you can easily write your own function, e.g.:
from Bio import SeqIO

def to_dict_remove_dups(sequences):
    return {record.id: record for record in sequences}

print(to_dict_remove_dups(SeqIO.parse("test.fa", "fasta")))

Output:
{'1': SeqRecord(seq=Seq('CG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[]), '2': SeqRecord(seq=Seq('AAA', SingleLetterAlphabet()), id='2', name='2', description='2', dbxrefs=[])}

Alternatively, you Biopython's .to_dict() accepts a key_function argument, e.g. if you wish to make all the headers unique
from Bio import SeqIO
from itertools import count

c = count()

print(SeqIO.to_dict(SeqIO.parse("test.fa", "fasta"),
                    key_function = lambda x: f"{x.id}_seqNum{next(c)}"))

Output:
{'1_seqNum0': SeqRecord(seq=Seq('ATG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[]), '2_seqNum1': SeqRecord(seq=Seq('AAA', SingleLetterAlphabet()), id='2', name='2', description='2', dbxrefs=[]), '1_seqNum2': SeqRecord(seq=Seq('CG', SingleLetterAlphabet()), id='1', name='1', description='1', dbxrefs=[])}


Selection of differential expressed genes I'm working with RNA-seq data. I have 40 tumor samples and 5 Normal samples. Differential analysis with Deseq2 based on Fold change > 1.2 and alpha < 0.05 gave very low number of differentially expressed genes. Only 2 upregulated genes.

res <- results(dds, lfcThreshold = log2(1.2), alpha = 0.05)

To get more number of differential expressed genes I have few questions now

Instead of FDR < 0.05 can I use FDR < 0.1 (or) FDR < 0.5. Will there be any low confidence with this?
Can I select differential expressed genes only based on FDR < 0.05 without any fold change cutoff?
As I get very low number of DEGs with FDR < 0.05 & Foldchange > 1.2, Can I select DEGs based on Foldchange > 1.2 and p.value < 0.01 or 0.05 ?

PCA of the samples looks like this:

You can usually get away with FDR < 0.1, but that's as high as you can go. This all presumes you're doing follow-up experiments of some sort, of course. I guess you could increase the FDR more, but you're then really increasing the odds that your follow-up experiments will fail. Obviously increasing the FDR will decrease the confidence in the results, it's called "false discovery rate" for a reason.
Yes, fold-change cut-offs are mostly useful when you need to prioritize hits further. Small changes are almost impossible to validate with less-sensitive methods (qPCR) and less likely to be biologically relevant.
You're just going to waste your time if you ignore FDR.

Remember that all thresholds are intended to increase the odds of success in follow-up experiments. If you have infinite time and resources then you can set these to whatever you want, but you have neither so you better prioritize reasonably.

It's likely that you have a large batch effect in your data, effect sizes in cancer tend to be much larger than this.


Why are the clusters in the left side of a tSNE plot and not in the centre? I have no idea why I used the subsetdata in Seurat after using CCA the TSNEPlot shows the cluster all in the left side and merge， did I miss something? Here are the code and figure I got.

cluster <-SubsetData(object = immune.combined, ident.use= "3")

cluster <- NormalizeData(object = cluster, normalization.method = "logNormalize", scale.factor = 10000)

cluster <- FindVariableGenes(object = cluster, mean.function = ExpMean, dispersion.function = LogVMR, x.low.cutoff = 0, x.high.cutoff = 6, y.cutoff = 0, do.plot = F)

cluster <- ScaleData(object = cluster, display.progress = FALSE)
cluster <- RunPCA(object = cluster, pc.genes = cluster@var.genes, pcs.print = 1:10, genes.print =5)

PCAPlot(object = cluster, group.by= "orig.ident", dim.1 = 1, dim.2 = 2)

cluster <- FindClusters (object = cluster, reduction.type = "pca", dims.use = 1:10, resolution =0.3, print.output = 0, save.SNN = TRUE, force.recalc = T)

TSNEPlot(object = cluster, do.label = TRUE,do.return = TRUE,pt.size = 2,label.size = 5, no.legend=FALSE)

I'm going to stick my neck out here. Overall PCA followed by tSNE is a very powerful method of unsupervised learning however, I agree with @GWW,

cluster <- ScaleData(object = cluster, display.progress = FALSE)

If this is a max-min rescale transformation (common in supervised learning) then the outliers will badly affect the clustering. You can remove outliers that are over 2 sd from a formal analysis and re-run the calculation. That is just the risk of a max-min rescaling.

The one issue I'm concerned about is the initial quasi-standardisation followed by log-transformation,

log1p(value/colSums[cell-idx] *scale_factor)

In this case you are multiplying by 10000, this is to remove the ratio (used as a quasi-standardisation) and then taking the log ??? That is not the way to transform a ratio, it would be better to take the reciprocal of the ratio. There are other (better) methods, a formal standardisation would be useful, before a transformation is considered. Possibly if that was done correctly you wouldn't have the problem of outliers potentially crunching your data.

If you fiddle with the parameters you might get something - but its a symptom not the cause of the problem. However, if removing the outliers works (I have not checked the ScaleData function), then you are sorted and all else is forgotten (at least for this data set). The information I read was here.


Delete chains from a protein We have a docking target protein with multiple chains (A,B and C). However, I just need chain A. How can I delete the 2 other chains from the docking target protein.

If the tag r is incorrect, the easiest way is to download and install PyMOL: https://pymol.org/2/
Drap and drop your pdb file and type remove not chain A.
Then save molecule via the menu as pdb or type save some_file_name.pdb (will be in your home root folder).


Remove repetitive region vcf file using repeatmasker bed file I have a 1000 genomes vcf file for chromosome 14. And I want to remove variants in repetitive regions flagged by repeatMasker. I have a bed file for those repeats. I downloaded that file from UCSC. How can I eliminate/remove those sites in vcf file using repeat masker's bed file? Insights will be appreciated.

A couple fast alternatives:
$ bedops --chrom chr14 -n 1 <(vcf2bed < variants.vcf) <(sort-bed repeats.bed) > answer.bed
$ bedops --chrom chr14 -n 1 <(vcf2bed < variants.vcf) <(rmsk2bed < repeats.rmsk) > answer.bed

Using --chrom chr14 focuses set operations on that chromosome.


Is there a way to save a full pdb file in PYMOL after uploading biological unit? pdb code: 2obt

The pdb file has a coordinates of a monomer.

I have used following commands to generate an assembly, or in other words to upload one more biological unit so that the protein would be a dimer:

set assembly, ""
fetch 2obt, asu, async=0

set assembly, "1"
fetch 2obt, assembly1, async=0

After generating the functional dimer in  pymol under name assembly1, I have an important issue that I cannot save a dimer structure in pdb file with dimer coordinates. Does anybody has an idea how to go around it or is it even possible to do that?

So on the GUI form of PyMOL on the right hand side you have a list of what is ambiguously called objects in PyMOL (models elsewhere), which can be enabled or disabled by clicking on their name.
The name also has a 1/1 in green or similar. This indicates the state. In your case you have 2 states.

To split a multistate object use

split_state prot_A

To combine two objects one uses the command

create new_name, prot_A or prot_B

Now, this has two caveats which if violated result in a nasty piece of modern art:

if the chains, segment ids and residues are the same
if you had two overlapping macromolecules

To deal with the first, change the chain id on one of them.

alter prot_B, chain='B'
   sort

Segment identifier is a weird thing that most comp-biochem tools dislike, so I like to remove that

alter prot_B, segi=''
   sort

Now when you go to save your new_name object, in the popup of the GUI you get a dropdown menu defaulted to enabled (visible objects that have not be clicked to be disabled). There you can choose new_name.

Aligning two

If someone that is in not the OP author has a template protein and needs to make a dimer of another based on that, the process is the same.

create twin, prot_B
align twin, prot_A and chain B
alter twin, chain='B'
sort

SDF files

If someone else has a multicompound sdf where the state is the header name, I strongly advise not to use split_state. Yes, the python pymol module (not the app) handles states in a really weird way, but states can be selected with state 1 selection algebra and state is an argument for most commands that do not accept state selector.


Running MaxQuant on Linux I'm trying to run MaxQuant on a linux laptop, hower I'm constantly running into problems with MaxQuant crashing at different stages

I tried to use

mono 6.8
mono 5.4.1

on

centos
ubuntu

each with MaxQuant

1.6.3
1.6.7
1.6.11

(everything using docker), but I had no success so far.

I also tried to use the ahead of time compilation of mono, also without sucess.
MQ mostly stalls at the "preparing searches" stage.

Has anyone managed to get MQ running properly on linux and can tell me, which linux distro/version, mono version and MQ version you are using?

Any other help is of course also welcome

There is a conda package for maxquant that should work on most platforms (I've used it on CentOS and I think it's been used on Ubuntu as well). There's also a docker container version of that that we (Bioconda) have posted to the biocontainers project that should also work.

Please note that large parts of a maxquant analysis are single threaded, regardless of what you specify. We've also had single analysis runs lasting over a month (even when using 90 threads), so the program can be very slow if your search space is large.


Differential gene expression bias due to effect of an individual sample I am analysing a human single cell RNA seq experiment, where we have 4 groups, four samples each. Data has been analysed using Seurat, with the canonical workflow. I have tried DE using various methods (Wilcox, MAST, LR), and for each of them the top genes (biggest fold change, lowest adjusted p) are genes over expressed by specific samples, e.g.:

Looking at the average expression for these genes per each subject, and plotting on a Boxplot, it does not seem to be a genuine difference (e.g. for the gene above):

Besides the subject, there is no other parameter (e.g. age, sex) that is biased across groups and that shows up as differentially represented in the PCA plot.

Is there a DE method that takes into account this behaviour and corrects for it?
How can these results be presented?
Besides increasing number of subjects (not really feasible), is there a solution to analyse these data?

Thanks a lot!

Leave-one-out cross validation is an iterative method that leaves out one sample until each sample has been left out once. You can then identify the intersection of the differentially expressed genes from the sample space of repeated analyses to produce a set of genes that are consistently differentially expressed independent of individual sample biases.
more info here:
http://efavdb.com/leave-one-out-cross-validation/


How is prior.count used by edgeR's cpm edgeR's cpm function has an argument called prior.count. Based on my understanding of the documentation, it is supposed to be adding a fixed number per sample which is proportional to the library size of said sample. Average of all numbers added to all samples would be equal to prior.count.

However, looking at actual data this does not seem to be the case. Given an imaginary data frame of

df = data.frame(a = c(1,2,3,0),c = c(1,2,3,0)*3, b = c(1,2,3,0)*2)

We can try to calculate log cpms by doing

logCPM = cpm(df,log = TRUE,prior.count = .5)

We can also calculate regular cpms by doing

CPM = cpm(df)

If we wanted to see, what is the number that is added to the CPMs before getting logged, we could do

difference = 2^logCPM - CPM

To make it look pretty lets use tibble

tibble::as.tibble(difference)

# A tibble: 4 x 3
                     a                    c                    b
                 <dbl>                <dbl>                <dbl>
1  25641                25641                25641
2  12821                12821                12821
3 -    0.0000000000582 -    0.0000000000582 -    0.0000000000582
4  38462                38462                38462

Here we see that the number that is added has no relation to the library size of the sample. What I want to learn is, how is this number that is added to each cpm is calculated based on prior.count.

I have been digging through the code but that part goes into C++ territory which hinders my understanding of what is going on

The prior count ends up getting scaled by the ratio of a library size to the average library size and then multiplied by 2 before getting added to each library size (I'm sure there's a good reason for that, but I don't know what it is). Using the example df data frame in your post, let's walk step-by-step through what cpm() is doing:

df = data.frame(a = c(1,2,3,0),c = c(1,2,3,0)*3, b = c(1,2,3,0)*2)
prior.count = 0.5
# First, we need to calculate a library size
lib.size = colSums(df)
# Calculate the average library size and the adjusted priors
ave.lib = mean(lib.size)
adjusted.prior = prior.count * lib.size / ave.lib
# Update the library sizes
adjusted.lib.size = lib.size + 2*adjusted.prior
# Now we can compute the CPM
customCPM = t((log(t(df) + adjusted.prior) - log(adjusted.lib.size) + log(1000000))/log(2))

The matrix transposition looks really ugly in the last step, but it's needed so R adds things like adjusted.prior to each row of df.

Note that the prior count is only used if you compute log(CPM), it's unused otherwise.


Do I need to fill in missing side chains and loops before protein-ligand docking for VLS? I am going to screen a set of small molecules against a CYP450 2C9 protein target. The structure I am using exactly is 5A5I. Do I need to fill in missing side chains and loops before performing VLS? I assume while preparing the protein I should leave water molecules in the structure as there is a bond between $Fe^{2+}$ and water that is involved in the binding of a ligand.

In general
Most if not all docking applications will tolerate a jump.
Furthermore, you really don't want to dock anything to the flexible regions even if you are using a program that does properly protein flexibility.
Furthermore, if your protein adopts a closed conformation upon native substrate binding (not really the case here), you should start with that for better docking results.
Cytochrome P450
So most Cytochrome P450s are strange though and in many those loops play a large part in their catalysis. I am on a paper under review (bioarxiv) that deals exactly with this in a different cytochrome (P450-BM3, the biotech one). Namely even if the active site is fairly large for an active site in the core of an enzyme, the ligand (testosterone) did not really rotate around and instead the specificity of oxidation was controlled by how the ligand got in. So those crystallographically unresolved loops act as gatekeepers for what can bind. In terms of VLS, one aims for blocking off the cofactor, i.e. the heme. So you can do docking without worrying about them in a first instance, but there is a possibility that a lot of active drugs instead interact with the loops.
Certainly, dynamic undocking and other MD methods may be required a few iterations in, to see why some candidates failed and why others worked.
An interesting unofficial further possible spin is the fact that solvent-exposed sequential residues Asp-Pro can hydrolyse under acidic conditions, in a mechanism similar to isoaspartate formation, and several cytochromes have these on the entrance loops and have weird kinetics and smeary SDS-PAGE gels.
Lastly, you are totally correct in your observation. The water binding the heme in the crystal is very special. It is taking the place of oxygen. That could be modelled as is, changed to a pre-encounter O2, or event changing the heme to a heme bound to a radical oxygen —note that some forcefields (CHARMM, Rosetta Ref2015) do not like this weird valency iron.


ERROR: (subscript) logical subscript too long I am trying to remove a row of a dataframe in R however because the dataframe is so large when I use:

rawCounts <- rawCounts[rawCounts != geneID, ]

I get the error:

ERROR in rawCounts[rawCounts != geneID, ]: (subscript) logical subscript too long

Is there another way of removing a specific row in a large dataframe that is not so memory intensive?

Thanks in advance.

I have quite a bit of memory, memory.limit() returns:

> memory.limit()
[1] 16318

I am answering my own question here.

I realised I should have been using

rawCounts <- rawCounts[rownames(rawCounts) != geneID, ]

Problem solved.


Error pasing the following line in pdb I converted an .sdf file containing around 50,000 anti-viral agents to individual .pdb files using OpenBabel. However, when I tried opening the .pdb files with Autodock I kept getting the error below. What do you think is causing the error.
Thank you for your time and patience.
Error Message:

“ Error pasing the following line in pdb:
HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000 ."

A PDB file has strict spacing rules as documented here.
The coordinates should be:
31 - 38       Real(8.3)      x             Orthogonal coordinates for X.
39 - 46       Real(8.3)      y             Orthogonal coordinates for Y.
47 - 54       Real(8.3)      z             Orthogonal coordinates for Z.

So trying this on your line:
'HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000'[30:39]

I get 9.051 0., which is not a number.
The whole line is formatted incorrectly however.
So as you say you are using Autodock why not pass a mol2 to prepare_ligand4.py instead which is an acceptable format?
Alternatively, you could convert your multi-entry sdf to pdb files via RDKit in Python:
suppl = Chem.SDMolSupplier(sdfile)
for mol in suppl:
      # mol = Chem.AddHs(mol) # has hydrogens?
      # AllChem.EmbedMolecule(mol) # is 3D?
      AllChem.ComputeGasteigerCharges(m) # sdfs do not have partial charges
      Chem.MolToPDBFile(m, f"{mol.GetProp('_Name')}.pdb")


how is the DNA Integrity Number (DIN) calculated in Bioanalyzer/TapeStation? For DNA/RNA quantification machines like the Bioanalyzer or TapeStation, the DNA Integrity Number (DIN) or RNA Integrity Number (RIN) numbers are quoted as a measure of the fragmentation of the material.

How is the DNA Integrity Number (DIN) calculated in Bioanalyzer/TapeStation? Is there any open-source software to parse/collate the data from the machines?

Note that DIN (DNA Integrity Number) and RIN (RNA Integrity Number) are different scores. I've had trouble in the past finding code or formulas for calculating both of these scores, which is very frustrating considering how frequently they are used in research papers and for NGS QC. The closest that I have been able to get to this are demonstrative graphs for various RIN numbers and DIN numbers, with a description of the variables that are included in the model.

The code that @719016 provided appears to attempt to calculate the RIN, not the DIN.

Just in case it's helpful, here is an extract of that code, only including lines that seem to be directly involved in calculating the actual B_RIN values. This is not runnable, because it depends on particular input files, and it's still pretty unintelligible:

myfiles<-list.files()
results<-myfiles[which(substr(myfiles,69,79)=="Results.csv")]
n <- which(substr(myfiles,69,79)=="Results.csv")
targets <- read.delim(myfiles[n[1]], sep=",",header=F,skip=14,as.is=T)
targets2<- read.delim(myfiles[n[1]], sep=",",header=F,skip=0,as.is=T)
str3<-paste0(substr(targets2[1,2],1,67),"_S")
chip_num<-length(n)
samp_num_chip<-length(grep(str3,myfiles))
samp_num<-length(which(substr(myfiles,69,74)=="Sample"))
samples1<- matrix(NA, ncol=chip_num, nrow=samp_num)
for(i in 1:chip_num)
{
    ##mod chip numb
    targets2<- read.delim(myfiles[n[i]], sep=",",header=F,skip=0,as.is=T)
    str3<-paste0(substr(targets2[1,2],1,67),"_S")
    samp_num_chip<-length(grep(str3,myfiles))
    ##mod sample numb
    str1<-substr(targets2[1,2],1,67)
    for(j in 1:samp_num_chip)
    {
        iii<-j
        str5<-paste0(str1,"_Sample",iii,".csv")
        samples1[j,i]<-str5
    }
}
samples1<-samples1[grep("2100",samples1)]
##RIN Calculator#######refined peakfinder####################
##read bioanalyzer data into a matrix called dta
##since the total RNA and mRNA assays run differently
##skip more lines for the total RNA assay
elec <- read.delim("ba_lane.txt", sep="\t",header=T,as.is=T)
qc.mat <- matrix(NA, ncol=1, nrow=nrow(elec))
dta <- matrix(NA, nrow=1060, ncol=nrow(elec))
for(i in 1:nrow(elec))
{
    x <- read.csv(as.character(elec[i,2]), header=F, skip=18, nrows=1060)
    dta[,i] <- x[,2]
}
time<- x[,1]
for(samples in 1:length(samples1))
{
    ##store max fluor in max.peaks matrix
    lad<-samples
    ti1<-17
    ti2<-65
    p.row<-length(dta[which(time==ti1):which(time==ti1+1),lad])*(ti2-ti1+1)
    peaks <- matrix(NA, nrow=p.row, ncol=1)
    max.peaks <- matrix(NA, nrow=length(seq(ti1,ti2,.05)), ncol=1)
    for (i in 1:length(seq(ti1,ti2,.05)))
    {
        ii<-seq(ti1,ti2,.05)[i]
        x0<- ii
        x0<-round(x0,digits=3)
        y0<- ii
        x<- ii+.05
        x<-round(x,digits=3)
        y<- ii
        max.peaks[i]<- max(dta[which(time==x0):which(time==x),lad])
    }
    ##human or mouse
    cent<-max(max.peaks)*.6
    ##flatworm
    ##cent<-max(max.peaks)*.05
    ##define the middle of the segments
    max.avg<-seq(ti1+.025,ti2+.025,.05)
    ##find time associated with peaks > 5 fu and fill seconds matrix with them
    ##12.5% of max fu incase chip is fu
    sizes <- which(max.peaks>cent)
    gt35<-which(max.avg[sizes]>35)
    lt55<-which(max.avg[sizes]<55)
    tfr<-which(duplicated(c(gt35,lt55))==TRUE)
    area<-cbind(max.peaks[sizes-1][c(gt35,lt55)[tfr]],
                max.peaks[sizes][c(gt35,lt55)[tfr]],
                max.peaks[sizes+1][c(gt35,lt55)[tfr]])
    time2<-cbind(max.avg[sizes-1][c(gt35,lt55)[tfr]],
                 max.avg[sizes][c(gt35,lt55)[tfr]],
                 max.avg[sizes+1][c(gt35,lt55)[tfr]])
    t1<-time2[which(area[,1]<cent),1][1]
    t2<-time2[which(area[,1]<cent),1][2]
    t3<-time2[which(area[,3]<cent),3][1]
    t4<-time2[which(area[,3]<cent),3][2]
    s18a<-ifelse(is.na(t1,-30,round(t1,digits=1)))
    s18b<-ifelse(is.na(t3,-35.05,round(t3,digits=1)))
    s28a<-ifelse(is.na(t2,-30,round(t2,digits=1)))
    s28b<-ifelse(is.na(t4,-30.05,round(t4,digits=1)))
    eta<-which(time==s18a)
    etb<-which(time==s18b)
    tea<-which(time==s28a)
    teb<-which(time==s28b)
    s18<-sum(dta[eta:etb,lad])
    s28<-sum(dta[tea:teb,lad])
    qc.mat[samples]<-  (-1*exp((s28/s18)*-1)+1)*10
}

Mostly based on that last line, the RIN seems to be associated with the ratio of the 18s peak to the 28s peak.


Major variability in speed of BLAST between proteins I'm running blast+ searches for a list of proteins using specific parameters.

This is my search

 blastp -query protein.fasta -remote -db refseq_protein -out
protein_refseq_blast.txt -outfmt '6 qseqid sseqid qstart qend length evalue' -
entrez_query 'txid2 [ORGN] or txid4751 [ORGN]' -evalue 200000 -max_target_seqs 10000

For some proteins, this runs moderately quickly, ~5-10 mins. For other proteins, this command never stops running, even after an hour or so. I'm wondering if this variability in speed is due to any part of this choice in parameters, is simply variability in speed for the BLAST servers. Is it common for BLAST searches to take this long? Is there a way to run all of them in parallel using these parameters if this is the case?

Thanks.

Try changing you evalue threshold to something like 0.001.

If you blasted a sequence against a database containing 10,000 sequences and got a hit with an evalue of 1, this would mean that, given the size of your database and length of protein, you would expect to find a match as good as (or better than) the hit returned once purely by chance. If the evalue threshold in your query was set to 0.01 (for example) then this hit would not be returned.

If you got a hit with an evalue of 200,000 this would mean you would expect to see 200,000 matches with a similar score simply by chance. If your evalue threshold were set to 200000, then all of the hits with an evalue <=200000 would be considered significant. If you set it to 0.001 then only those hits with evalues <= 0.001 would be considered significant.

This means that (depending on the size of the data base) almost every sequence in your query set will be considered as significant given your threshold, which will slow down your search hugely.

See the the section on E values here


Are there any rolling hash functions that can hash a DNA sequence and its reverse complement to the same value? A common bioinformatics task is to decompose a DNA sequence into its constituent k-mers and compute a hash value for each k-mer. Rolling hash functions are an appealing solution for this task, since they can be computed very quickly. A rolling hash does not compute each the hash value from scratch with each k-mer: rather it updates a running hash value using an update strategy and a sliding window over the data.

It's also very useful for many applications to have a k-mer hash to the same value as its reverse complement. Unless the data were generated using a strand-specific sample prep, it's impossible to distinguish a k-mer from its reverse complement, and they should be treated as the same sequence.

Are there any rolling hashes that will map reverse complements to the same value? If not, how would we develop such an algorithm?

UPDATE: Ideally the hash function would be able to support k > 32, which would be lossy unless using something larger than a 64-bit integer.

ANOTHER UPDATE: I don't think it's necessary to store both the running k-mer and its reverse complement in a single value. If storing two k-mer strings and/or two hash values makes this easier, I'm totally cool with that.

Let f and r be two integers. They always keep the k-mer on the forward and reverse strand, respectively. At a new base c in a proper 2-bit encoding, we update the two integers as follows:

f = (f<<2|c) & ((1ULL<<2*k) - 1)
r = r>>2 | (3ULL-c)<<2*(k-1)

With this updating rule, f keeps the forward strand k-mer ending at c and r is f's reverse complement. The hash of the k-mer can be min(f,r). If integer operations take constant time, computing the hashes of all k-mers of a sequence of length L takes O(L) time.

The above only works if k-mer fits a word, which means k<=32 on x86. For long k-mers, you can overload bit operators in C++. In case of 32<k<=64, a simpler solution is to split into lower and higher bits:

c1 = c&1, c2 = c>>1
f1 = (f1<<1|c1) & ((1ULL<<k) - 1)
f2 = (f2<<1|c2) & ((1ULL<<k) - 1)
r1 = r1>>1 | (1ULL-c1)<<(k-1)
r2 = r2>>1 | (1ULL-c2)<<(k-1)

This is twice as slow as the k<=32 version. It is possible to implement these lines with SSE2 intrinsics, but that is overkilling.

When you want to hash a long k-mer into fewer bits, there are a few options:

xor:      f^r
min:      min(f,r)
min+hash: hash(min(f,r))

where hash() is a generic randomization hash function, which could be Murmur3, Wang's integer hash function, etc.

The following is a microbenchmark. Here, we squeeze all k-mers into a 32-bit hash table. collision_rate = 1 - #distinct_hash/#distinct_kmer. This is not the best metric, but it somehow measures the randomness, which is better than nothing.

====================================================
 Algorithm    data  max-k  k   %collision  CPU time
----------------------------------------------------
 xor          chr11  32    31     8.5%       0.9s
 xor+Wang     chr11  32    31     9.6%       1.3s
 min+Wang     chr11  32    31     1.4%       1.3s
 min+Wang     chr11  64    31     1.4%       1.9s
 min+murmur3  chr11  inf   31     1.4%       5.7s
 min+Wang     chr11  64    51     1.5%       2.1s
 min+murmur3  chr11  inf   51     1.5%       6.8s
 min+Wang     chr11  64    63     1.5%       2.1s
 min+murmur3  chr11  inf   63     1.5%       7.5s
 radix sort   chr11  -     -      -        5.6-7.3s
----------------------------------------------------
 min+Wang     hg38   64    63    26.2%        37s
 min+murmur3  hg38   inf   63    26.2%       192s
 radix sort   hg38   -     -      -         ~138s
====================================================

Observations:

Don't use XOR. In fact, XOR hashes all palindromes to 0xffff.... This is already worrying enough.
Unless you want to work with very long k-mers, don't use generic string hash functions like FNV and Murmur. Murmur is even slower than radix sorting all hashes.
For k<=64, min+Wang is the best here. It is fast and simple to compute and has randomness nearly as good as murmur.


What is a good RNA seq normalization method that allows for across sample comparisons and between transcripts What is a good RNA seq normalization method that allows for across sample comparisons, and allows between transcripts comparisons as well? I read that TMM for example allows across sample comparisons but not between transcripts. On the other hand, I read that TPM is not good enough for between sample comparisons

Unfortunately, there doesn't seem to be any such metric. We use counts two-way normalized across genes and samples but even that is approximate as there is no way to ensure that all genes being counted are definitely expressed across all samples, and that all samples are sequenced to the exact same depth and that there are no batch effects and that biological differences are reflected with full fidelity.


How to exclude the repetition of gene-gene correlation calculation in python? I am trying to make a paired matrix of gene-gene correlation. Considering that I have a huge matrix (13000 genes and 900 samples) and for some reasons I don't want to decrease the number of my genes, my gene-correlation matrix would be 13000*13000 and my paired matrix will become 169 million *4 (Column 1: Gene 1; Column 2 : Gene 2; Column 3: Correlations; Column 4: P-values) . In this case, I have to exclude unnecessary calculations as much as I can. I have excluded the situation that Gene 1 = Gene 2. But I couldn't find a way to exclude the condition that "Column 1: Gene 1 ; Column 2: Gene 2 = Column 1: Gene 2; Column 2: Gene 1 ". To make a long story short, correlation between G1 and G2 is equal to G2 and G1. It is like calculating just lower section of diag in a symmetric matrix. I would be grateful if anybody help me in this case. I have enclosed my python codes here:

import pandas as pd
import numpy as np
import scipy
import math
import openpyxl
from openpyxl import Workbook
from scipy.stats import spearmanr

din=pd.read_csv('m_test.csv', index_col=0)

out=pd.DataFrame()
outdf=pd.DataFrame()

for g1 in din.index:
    for g2 in din.index:
        temp=din.loc[[g1, g2]]
        if g1==g2:
            next
        else:
            spR, spP=spearmanr(temp.loc[g1], temp.loc[g2])
            frame={'g1':[g1], 'g2':[g2], 'spR':[spR], 'spP':[spP]}

            out=pd.DataFrame(frame)

            outdf=pd.concat([outdf, out])

You don't need g2 to go from 0 to din.index, just from 0 to g1 - 1. This way, you'll end up calculating just for the "lower triangle":

g1-> 0   1   2   3
g2
0                    #(0,0 to 0,1-1 = nothing)
1    X               #(1,0 to 1,1-1 = 1,0)
2    X   X           #(2,0 to 2,2-1 = 2,0; 2,1)
3    X   X   X       #(3,0 to 3,3-1 = 3,0; 3,1 and 3,2)
4    X   X   X   X   #(4,0 to 4,4-1 = 4,0; 4,1; 4,2 and 4,3)