How does DeepVariant construct RGB images from DNA sequences? DeepVariant is a pipeline to call genetic variants from DNA sequencing data.
A major step, before feeding the CNN, is to translate these DNA sequences into images. It's unclear why and how Google constructs the RGB images from the DNA data. Obviously, DNA is a string over an alphabet with the characters: {A, T, C, G}.
It is even hard to understand how the mapping works based on the source code of their unit tests.
In their figure from the paper: A is Red, C is Green, G is blue, and T is Yellow (G+R), but this is still unclear how they construct the 3xNxN image.

EDIT from google's blog:

In this article we will show the six channels in a row, but in
DeepVariant they are encoded as six layers in the third dimension,
giving each tensor a shape of (100, 221, 6) corresponding to (height,
width, channels). The variant in question is always in the center of
each pileup image, here marked with a small line at the top.
Channels are shown in greyscale below in the following order:
Read base: different intensities represent A, C, G, and T.
Base quality: set by the sequencing machine. White is higher quality.
Mapping quality: set by the aligner. White is higher quality.
Strand of alignment: Black is forward; white is reverse.
Read supports variant: White means the read supports the given
alternate allele, grey means it does not.
Base differs from ref: White means the base is different from the
reference, dark grey means the base matches the reference.

Actually the paper has made it clear how they did it. You just have to read the supplementary materials closer.

In their figure from the paper: A is Red, C is Green, G is blue, and T
is Yellow (G+R), but this is still unclear how they the 3xNxN image.

In RGB, each dimension is an NxN image. Since you have three dimensions, so it's 3xNxN. The red dimension was used to encode the nucleotide bases. The green dimension was used to encode quality scores. Finally, the blue dimension was used to encode the strand information.

Obviously, DNA is a string over an alphabet with the characters: {A,
T, C, G}.

This is easy. Lot's of ways. You could do a one-hot encoding or what DeepVariant used:
def get_base_color(base):
    base_to_color = {'A': 250, 'G': 180, 'T': 100, 'C': 30}
    return base_to_color.get(base, 0)


Is it possible to use spike-ins with droplet technologies? By reading Bacher & Kendziorski (2016), I found the following statement:

  recent droplet technologies are not yet able to accommodate spike-ins

What's the reason why current droplet technologies cannot use spike-ins? By searching through the available methodologies I found this question on 10X community, but still I cannot understand what is the problem created by using spike-ins with droplet technologies.

Bacher, R. & Kendziorski, C. Design and computational analysis of single-cell RNA-sequencing experiments. Genome Biol. 17, 63 (2016).

I believe it is mostly a cost issue as spike-ins are "expensive" to use given the large number of cells that are/can be sequenced via droplet based methods (as compared to plate-based ones). The 10x link that you provided actually states that this is possible with 10x and Klein et al. used ERCC spike-ins in 953 droplets (Cell. 2015 May 21;161(5):1187-1201. doi: 10.1016/j.cell.2015.04.044.).


How to extract metadata from NCBI's experiment? I want to extract metadata from experiment SRX1596422. I have sratoolkit and Entrex direct but not sure what line of code to use. Here's the link for the experiment: experiment SRX1596422

You'll want something like:
esearch -db sra -query SRX1596422 | efetch -format runinfo

This will produce a CSV output to the screen with columns containing the meta information available in SRA:
Run,ReleaseDate,LoadDate,spots,bases,spots_with_mates,avgLength,size_MB,AssemblyName,download_path,Experiment,LibraryName,LibraryStrategy,LibrarySelection,LibrarySource,LibraryLayout,InsertSize,InsertDev,Platform,Model,SRAStudy,BioProject,Study_Pubmed_id,ProjectID,Sample,BioSample,SampleType,TaxID,ScientificName,SampleName,g1k_pop_code,source,g1k_analysis_group,Subject_ID,Sex,Disease,Tumor,Affection_Status,Analyte_Type,Histological_Type,Body_Site,CenterName,Submission,dbgap_study_accession,Consent,RunHash,ReadHash
SRR3184375,2016-02-24 00:00:00,2016-02-23 23:14:45,556018223,111203644600,556018223,200,57997,GCA_000001635.1,https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos2/sra-pub-run-7/SRR3184375/SRR3184375.1,SRX1596422,,WGS,size fractionation,GENOMIC,PAIRED,0,0,ILLUMINA,Illumina HiSeq 2000,SRP070623,PRJNA312530,,312530,SRS1306521,SAMN04492320,simple,10090,Mus musculus,MCNT_E1,,,,,female,,no,,,,,THE SCRIPPS RESEARCH INSTITUTE,SRA356784,,public,52A6F836C13F2D15755EDE3E52A18D78,240E7130CFE36405818DA5EE2A10E24C

The more relevant fields are: Subject_ID, Sex, Disease, Tumor, Affection_Status, Analyte_Type, Histological_Type, Body_Site. Note that not all of these are available for that particular sample.


Does SBOL support timing and threshold value parameters? Like the propagation delay and threshold value parameters in electronic circuits, does SBOL format support any such parameters/glyph to represent the propagation delay (time to trigger the output protein once the input is provided) and threshold value (minimum input concentration required to trigger the model's output) of genetic models?
If currently not supported, can these parameters be incorporated in the SBOL visual representation which would give reader an idea (just by looking at the schematic diagram) about the timing behavior of a model?

Note: I will use SBOL3 in my descriptions; the answer can be equivalently translated to SBOL2
At present, SBOL doesn't have any explicit recommendations for how to represent these types of timing parameters. This is because there is not yet a clear scientific consensus on whether or not such parameters are a good way to represent genetic circuits.
That said, one of the advantages of using SBOL is that its extensible RDF-based model allows additional information like this to be added into a model. Two good ways of doing so are:

If the parameters are simple and can stand alone as measurements, then they can be added to designs as a Measure object with appropriate units and type. For example, if you have biological devices with a 90 minute propagation delay, then on the Component representing the device you can add a Measure with numerical value 90, units of om:minute.  You'd likely want two types, one of SBO:0000225 (Delay), and a custom type for your own notion of propagation delay, since I don't think that's included in any ontology yet.
If the parameters are more complex and/or inter-related, e.g., a full data-sheet, then you'd likely want to instead create one or more extension classes that are used to represent the information. A library that makes this easy is sbol_factory, though one can do it with any of the SBOL libraries.

This becomes more complicated if you want to talk about behavior in different environmental contexts, but the basic principles are still the same.
Finally, as for SBOL visual: again, there's no explicit representation, but the model is extensible, and this information can be added as any sort of text or custom glyph as long as it doesn't conflict with ones already in the standard.


Hierarchical clustering for outlier detection - single cell RNASeq & WGCNA I wish to conduct WGCNA on a single cell RNASeq dataset and, when choosing the optimal beta parameter after running pickSoftThreshold am presented with a rather strange plot:

I have normalized the data using NormalizeData in Seurat using the LogNormalize method.

One of my colleagues suggested using hclust as a method for outlier detection to see potential cells/samples that could be excluded. I think this makes sense however I haven't found much information on such procedures besides the following paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5845381/.

My question would be, is this procedure just be as simple as running a hierarchical clustering on the data and seeing which samples cluster by themselves (and then removing those based on a minimum cluster size criteria) or is their more to it? Also, for those familiar with WGCNA, why does the scale-free topolgy plot generate such strange looking values, is this indicative of outliers being present in the data?

Thank you for your time.

For outlier identification I suggest using the sample network approach developed by Oldham et al. It basically amounts to constructing a inter-sample connectivity (assuming normalizedData contains the log-normalized data with genes in colunms and samples in rows)

k = colSums(cor(t(normalizedData))),

scaling the connectivities

Z.k = scale(k)

and then removing samples with high negative Z.k (e.g. Z.k < -4).

A slightly different approach is based on euclidean distances rather than correlations,

Z.k = scale(as.matrix(dist(normalizedData)))

and removing samples with high positive Z.k (since distance is a measure of dissimilarity).

Regarding the scale-free topology fit plot, I wouldn't worry about it too much, the fit index sometimes jumps around a bit. It certainly does not look like you have strong outliers or very heterogeneous data. But it is always a good idea to plot the sample clustering tree to check (1) that there are no strong clusters, and (2) that the outlier identification did something reasonable.


Getting System.IO.FileNotFoundException with MaxQuant running on Linux with mono I am getting this error with MaxQuant 1.6.14 running with mono (5.14.0.177=he1b5a44_1001).

System.IO.FileNotFoundException: Could not find file "/home/ubuntu/run-dir/combined/combinedRunInfo"

And I not sure why I am getting this error. The process has read and write permissions to that location.

...seems the mqpar file was not compatible with the maxquant version.
Make sure the mqpar.xml file shows the same version number as the MaxQuant binary.
MaxQuant is not really made for the command line. So, the errors shown are sometimes quite misleading.


Given a transcription factor, what genes does it regulate? I have a list of transcription factors and I am interested in finding out which which genes might be transcribed as a result of the formation of transcription factor complexes with that transcription factor.

Any ideas on good databases? I've seen that ENCODE hosts data from CHIPseq experiments, but I'm not sure if I can find conclusions regarding interactions from this site. I'm looking for a database with known or putative transcription/gene interactions.

Thanks in advance.

MSigDB has a collection (C3:TFT) of gene sets corresponding to transcription factor targets.

Harmonizome has functional terms for genes extracted from over a hundred publicly available resources.


Convert FPKM to RPKM I have my 24 samples analyzed by Hisat2>Stringtie, from paired end sequenced data. I have got FPKM data and am planing to convert them to RPKM. to be able to compare my samples with another already published.

1- How can I convert from FPKM to RPKM?

2- Is it an acceptable method to compare samples using RPKM between two different studies or I need more or different procedures?

I don't know if you have followed the discussions from the last years, but FPKM and RPKM are obsolete. Statisticians found that FPKM and RPKM are not the right values that can be used for proper statistics. It is a wrong way of normalizing RNA-seq data. In stead they suggest to use raw read counts, which can be used in tools like edgeR or DESeq2. Also, CPM would be better than FPKM or RPKM, for example you can use (log2) CPM values in limma with trend.

But to answer your questions.

You can't convert, you need to recalculate from raw data.
No not acceptable. To compare studies you need to correct for batch, which you can do with statistical programs such as edgeR and DESeq2.

The best advice I can give is to start with raw data, also from the already published data that you would like to compare your own experiment with. And continue with edgeR or DESeq2 or the like.


Difference in Calculation through GEO2R and difference in Article Here is microarray analysis where author uses GSE50632 series on an  Agilent microarray GPL17660 platform with 4 eligible samples which I guess are GSM1224991, GSM1224992, GSM1224993 and GSM1224994 among which two are healthy and two are patient samples.
However I found the difference in the calculations-
Through GEO2R analysis table it is-
ID  adj.P.Val   P.Value t   B   logFC   miRNA_ID_LIST   SPOT_ID
148000  1.86e-02    5.82e-04    3.47    -0.402958   3.91    hsa-miR-3195

However in the paper the calculated values are far different specially logFC value.

What could be the reason for this?

Here are two example of calculations of logFC validated in GEO2R analysis which depends if values are already log2 transformed. During GEOQuery analysis there is calculation to transform where expression set is checked for value distribution and log2 transformation is applied-
ex <- exprs(gset)
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
          (qx[6]-qx[1] > 50 && qx[2] > 0)
if (LogC) { ex[which(ex <= 0)] <- NaN
  exprs(gset) <- log2(ex) }

Example 1:
=> Sample Values are already log2 Transformed
n(TumorSamples) = 45
n(ControlSamples) = 46

a = average(TumorSamples) = 9.199053111
b = average(ControlSamples) = 4.988139348

logFC = |a - b| = 4.210914

Example 2:
=> Sample Values are not log2 transformed
n(TumorSamples) = 2
n(ControlSamples) = 2

a = average(log2(TumorSample1) + log2(TumorSample2) ... ) = (log2(84) + log2(21.8)) / 2 = 5.419287

b = average(log2(ControlSample1) + log2(ControlSample2) ... ) = (log2(624) + log2(665)) / 2 = 9.331306

logFC = |a - b| = 3.91202


Increase number of threads for GATK 4.0 HaplotypeCaller I am using GATK version 4.0, I tried to use multiple threads for calling variants using HaplotypeCaller using following command

gatk --java-options -Xmx90G  -nt 28 HaplotypeCaller -I output.bam -R wheat_ref.fa -O final.vcf

and the error is

'-nt' is not a valid command

I have also used --num_threads but again the error was the same, as written here this parameter is available in previous versions of GATK but I don't know what's its status in current version of GATK i.e GATK 4.0. Whether it is deprecated or replaced by some other parameter?

Is there any way of calling multiple threads for variant calling with GATK 4.0?

I am assuming you are using the non-SPARK version of the method. To specify the number of threads you wish to use with HaplotypeCaller, include --native-pair-hmm-threads (documentation). This will only parallelize the pair hidden Markov models (pair HMM) process. This is an argument for the method HaplotypeCaller itself, so it should be located after the method is called as opposed to one of the --java-options.

I do recommend trying out the SPARK implementations of the GATK tools. The whole method will take advantage of the threads you specify. I've been running SPARK locally (i.e. not on Firecloud) on computational servers. I've been impressed with the enhanced performance.


How to count the number of mapped read in 100-bp window from a BAM/SAM file Although I know how to get total number of mapped read using samtools flagstat (samtools flagstat file_sorted.bam) but I want to count total number of mapped read in a non-overlapping sliding window of fixed size (let's say 500 bp) only with respect to reference (hg19) (assigning each read only once by its start position).

An example of the desired output would be something like:

Chromosome    Region_Start     Region_Stop    Num_Unique_Reads
chr2          1                500            200
chr2          501              1000           202
chr2          1001             1501           230
chr2          .                 .              .
chr2          .                 .              .
chr2          .                 .              .
chr2          8001              8500           20

Can anyone suggest me how to perform that window operation?

one-liner

Here's a gritty one-liner to count the number of reads in a region if you have just one region that you want to investigate. Change the 1 in ($4 >=1) and the 500 in ($4 <=500) to set your window. Change "hg19" to your target sequence.

Note, this one-liner does not double-count reads because of uniq.

samtools view file_sorted.bam | \
    awk '{ if (($3 == "hg19") && ($4 >=1) && ($4 <=500)) { \
  print($1) }}' | sort | uniq | wc -l

python 3

If you need something a little more programmatic to get a known region you should consider trying the pysam package.

import pysam
samfile = pysam.AlignmentFile("file_sorted.bam", "rb")
region_set = set()
start = 100
stop =  500
counter = 0
for read in samfile.fetch('interesting_contig', 100, 500):
    region_set.add(read.query_name)
    counter += 1
print("{} unique reads in [{}, {}]".format(len(region_set), start, stop))
print("{} alignments in [{}, {}]".format(counter, start, stop))

For my bam file, this program outputs:

1618 unique reads in [100, 500]
2487 alignments in [100, 500]

python 3, non-overlapping windows

It is important to note that the samfile.fetch() function conforms to 1-based indexing, consistent with SAM files, unlike pysam's otherwise default 0-based indexing that matches the BAM-indexing.

import pysam
samfile = pysam.AlignmentFile("file_sorted.bam", "rb")
print(samfile.references)
print(samfile.lengths)
window = 500
print("ref\tstart\tstop\tuniq_reads")
for i in range(len(samfile.references)):
    refname = samfile.references[i]
    seqlen = samfile.lengths[i]
    for j in range(1, seqlen, window):
        stop = j+window-1 if j+window-1 < samfile.lengths[i] else samfile.lengths[i]
        region_set = set()
        for read in samfile.fetch(refname, j, stop):
            region_set.add(read.query_name)
        print("{}\t{}\t{}\t{}".format(refname, j, stop, len(region_set)))

For the small bam file I looked at the output was:

ref         start  stop  uniq_reads
OdonB       0      500   1668
OdonB       501    1000  2647
OdonB       1001   1500  129400
OdonB       1501   2000  598
OdonB       2001   2500  538
OdonB       2501   3000  872
OdonB       3001   3500  561
OdonB       3501   4000  49120
OdonB       4001   4500  48433
OdonB       4501   4806  28121


What is the pvalue from anova(glm1, glm2, test = “Chisq”) telling me in R when the glms have a negative binomial distribution? Having constructed two nested negbin glms in R, one has a dropped factor, We can do an anova with test = "Chisq" to see if there is a significant difference between the models. But what does the Pvalue as a result of this actually tell me?
anova(mnegbin1, mnegbin2, test = "Chisq")
Likelihood ratio tests of Negative Binomial Models

Response: count
                        Model    theta Resid. df    2 x log-lik.   Test    df      LR stat. Pr(Chi)
1 origin + substrate + sample 18.46502      1232       -2459.886
2          substrate + sample 18.46502      1232       -2459.886 1 vs 2     0 -1.500666e-11       1


Are crystal structure of Ligands used in volunteer computational docking like FightAIDS@home? In large Volunteer computing for molecular docking projects such as FightAIDS@home, do they use the original crystal structure of ligands in docking or unconfirmed structures of ligand?

I don't know that project's details, so cannot say officially, but undoubtedly it will be

crystal structures of protein or models of point mutants based on structures of protein.
Computer structure of the ligands based on force-field calculation

Protein structures

Firstly, HIV virus has been extensively studied by structural biologists, so most components have structures.

Secondly, docking on models is generally done only for threaded models based on closely related templates, say a human model of a mouse structure or a recently diverged pair (i.e. conserved proteins that split with the genome double duplication in bony fish). Whereas a given virus protein will likely have poor homology to a known structure and not be a good candidate.

A corner case is core point mutants... If you have a mutation in a protein it is likely going to be on the surface and not core: a model of a point mutant is likely very faithful to an empirical version —and generally you'd want to design drugs to a region that is not variable anyway as you'd want it to work universally. If you have a change in the active site to confer resistance say (i.e. in the core), there is generally a fitness loss to the organism (say lowered catalysis or thermostability) which needs to be compensated with other mutations (epistasis). As this is rarer and very interesting scientifically and medically, these resistant variants will most likely have been crystallised as the conditions for crystallization will most likely be the same.

Small molecule structures

Force field calculations for small molecules are most often valid (there computer Vs empirical RMSD is minute). The partial charges may be off, but that's not something crystals give anyway: In fact a ligand used in docking is (generally/ought to be) energy minimised beforehand. For docking getting charges wrong is the killer.

Crystallization of small molecules is used as validation of a synthesis route, so some may have structures, but they will not be used in a pipeline as a SMILES string is all you need.


Genotyping technologies do not maintain phase? Why exactly does genotyping not maintain phase? I guess I'm more confused on what's the difference between sequencing and genotyping, and as a result how phase is maintained in one, but not the other. This is an image trying to do phasing using genotyping methods.

I understand how phasing works, I just don't quite get what is the difference between genotyping and sequencing, and why one of them doesn't work.

Phasing means that you are able to distinguish the both copies of your target sequence. (Usually you have one from your mother and one from your father).

If you can do this, depends mainly on the sequencing technology you are using. For example in case of sanger sequencing, you cannot do this, because you (pre)amplify all your copies at once and detect the sequence at the same time at the same location.

In case of ngs you have read information at the end. One read represents the sequence of one original dna molecule. If you detect multiple variants in one read, you can be sure that they are from the same chromosome and thus they are in phase. Now you can go on and try to find other reads that have unambiguous overlaps with these variants on your read. If you find some, you can also say that these variants have their origin on the same chromosome.


Calculating L2 norm (Euclidean norm) to set of genes in time series I have gene expression RNA-seq data of 9 time points (0 hour, 2 hour, ..., 16 hour). I need a plot on which x axis is time point and each data point at a certain time point represents the L2 norm (Euclidean norm) of the difference between the gene expression vector at that time point and the previous time point.

This is the link of my expression data

https://drive.google.com/file/d/1crJSqVdGujy1tFyeoV6E7e-xO9SXteJZ/view?usp=sharing

If you read your matrix into R and saved it as m:

diffMatrix = m[,2:9] - m[,1:8]  # The change per-timepoint
l2n = apply(diffMatrix, 2, function(x) sqrt(sum(x^2)))  # The L2-norm

R also has a norm() function, which can directly produce the L2-norm, but it turns out that this is slower than the simple method shown here. If you click through to the other question you can find somewhat faster implementations of this, but anything will suffice for small datasets like yours.


Filtering imputed GWAS SNPs based on a MAF difference of 10% There are many posts on the web regarding QC steps pre and post-imputation.

Does applying below (new?) 10% MAF difference rule make sense, pitfalls?

Here is the process:

Get MAF for imputed set, using SNPTEST with flag -summary_stats_only
Convert imputed set to hard-calls using gtools with flag --threshold 0.9
If the MAF from step 1 and step 2 differs more than 10% than exclude the variant.

More info:

GWAS is 50K vs 50K case control samples.
This step is applied after info > 0.4 filter.

I cannot think of any principled rationale for choosing this filtering strategy.

However, I am going to take a guess that this filtering strategy is supposed to filter out SNPs for which imputation did not work well?  In that case the appropriate statistic to filter on is the INFO score as described here.

You might consider picking a higher threshold than 0.4 if you feel that you aren't filtering out enough variants.


How to convert gff to gtf? My annotation file is in .gff format. I would like to convert it to .gtf format or to know if there is a way to directly download the annotation file in .gtf format?

I am working on sequences from the P Falicparum 3D7.

You can use gffread to convert gff to gtf2, below is from the manual:

  In order to see the GTF2 version of the same transcripts, the -T
  option should be added:

gffread -E annotation.gff -T -o- | more

  The examples above also show that gffread can be used to convert a
  file between GTF2 and GFF3 file formats.


NCBI retrieve genome assembly with splicing or "modifications" I am trying to fetch genome assemblies that include splicing information.
E.g. Here, is a genome assembly of a yeast, where genes show variants like this: join(5928..7626,7685..7785,7857..7982).

I am retrieving the information from NCBI with the following commands:
efetch -db nuccore -id CP003820.1 -format tabular. I tried different formats and it changes but I don't get the join(...)

1) Am I correct and the join(...) indicates "splicing-like" modifications in yeast?
2) How can I get the join(...) (any other format is fine too) with efetch?

The page you link to is a GenBank Flat file, the default format used by GenBank. So, to get that format, which includes the join lines, use:

./efetch -db nuccore -id CP003820.1 -format gb

The join() lines give the different ways in which a sequence can be built from this gene. For example, one mRNA sequence produced by the locus you mention (CP003820) is:

  mRNA            complement(join(100..1497,1556..1767,1823..2274,
                 2322..2782,2846..3126,3206..3890,3958..5263,5322..5422,
                 5494..>5645))

The corresponding coding sequence is for the mRNA above is:

 CDS             complement(join(168..1497,1556..1767,1823..2274,
                 2322..2782,2846..3126,3206..3890,3958..5263,5322..5422,
                 5494..5645))

These aren't necessarily splicing "modifications", and aren't "variants", they are simply the way that the exons are joined to make an mRNA.

Note that since CP003820.1 is "Cryptococcus neoformans var. grubii H99 chromosome 1, complete sequence", there are multiple genes described in there and each will have its own (or multiple, if there is alternative splicing) join statements. Specifically, I see 1014 genes described in the file:

$ grep -Pc '^\s+gene\s+' CP003820.gb
1014

950 mRNAs:

$ grep -c 'mRNA.*join(' CP003820.gb
950

And 955 proteins:

$ grep -Pc 'protein_id=' CP003820.gb
955


Programmatically retrieve Accession List from NCBI I would like to programmatically retrieve the file "Accession List" from an NCBI page.
I'm running firefox, so I'm clicking the file and clicking Ctrl + Shift + E to try and figure out what the file endpoint is. But I can't seem to figure it out (I haven't done this much, I just remember doing it once for a previous project).
When I do this, I see two links popping up on my screen, as shown below:

I see that the first one has a field (I'm talking about the panel on the right) called "acc", which varies from one data dump to the next, so I'm presuming this is the link I need to use to retrieve the file (though I'm not sure), and just change the "acc" field depending on the data I want. Assuming this is right, I'm not sure how to go from here.
I've included the text on the right below for clarity.
GET
    https://www.google-analytics.com/collect?v=1&_v=j83&a=1402353855&t=event&_s=4&dl=https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA491191&o=acc_s%3Aa&ul=en-us&de=UTF-8&dt=Run Selector :: NCBI&sd=24-bit&sr=1280x720&vp=1232x601&je=0&ec=Link Click&ea=Metadata  &el=NONE&_u=SACAAEAB~&jid=&gjid=&cid=1350895533.1590937485&tid=UA-68648221-3&_gid=391407234.1594002343&cd1=run_selector&cd2=sra&cd3=&cd5=NO&cd6=run_selector:sra:&cd8=&cd9=&cd10=&cd11=&cd12=&cd13=&cd14=&cd15=&cd16=&cg1=run_selector&cg2=sra&cg3=&cg4=run_selector:sra:&z=508671370

Which accessions are you trying to get to? In the page you refer to, there is an option to download the list of accessions - one SRR run accession per line - by clicking the 'Accessions' button or the entire table by clicking the 'Metadata' button as highlighted in the screenshot below:

You can get the same data from the command line using the Entrez Direct. The command for that would be something like:
esearch -db sra -query 'PRJNA491191[bioproject]' | efetch -format runinfo

to download the table in CSV format.


When to account for the blacklisted genomic regions in ChIP-seq data analyses? We have heard in the group that it is important to keep track of and to filter artifact regions when analysing data from functional genomics experiments, especially ChIP-seq.

Here, we have seen pipelines that remove the ENCODE tracks i) before cross-correlations QCs, ii) after cross-correlations QC but before peak calling and iii) after peak calling.

We have noticed that removal of the tracks does not affect significantly cross-correlation and peak-independent QCs. However, we are not sure whether peak calling should be done on the filtered tracks or not?

Aside: Cross-correlation is largely meaningless, regardless of what some of the ENCODE folks might argue. When we process our DEEP samples we don't even look at that value.

Regardless, if you're using SPP/phantomPeakQual for cross-correlation then note that it already removes the highest peaks from your dataset before computing the cross-correlation (in fact, it can remove most of the actual peaks too, which makes one further wonder what it's actually telling you). I don't know that this is actually documented anywhere, it's something I noticed when going through the code while pondering whether to implement it in deepTools. But at least it's ignoring these regions already :)

In general, it's most convenient to just remove peaks overlapping blacklisted regions. In an ideal world you'd filter out the blacklisted reads before peak calling, but (1) this is really inconvenient (more time and disk required) and (2) I've never seen an appreciable gain in peak calling performance. In theory at least you should be losing sensitivity right around blacklisted regions if you don't remove reads overlapping blacklisted regions, but you have to ask yourself whether you want to trust such peaks anyway. For other QC steps, at least with deepTools we provide a parameter with every tool to specify a BED file of blacklisted regions to skip.

As an aside, there are many fewer blacklisted regions in more recent genome builds (GRCh38 and GRCm38 at least), so this is less of an issue in general with them.


Glocal\semi-local\Hybride Globale-Local alignment with Python I was looking for a simple way to do a glocal alignment.
The case I have is I have a small sequence which should be find in a bigger one, thus typically a glocal alignment.
Also I can not Install anything on the target machine (Python 3.6 with Biopython is already installed), all I can do is pip install.

My question is : are there any glocal alignment implementation out there ?

I was unable to found one and I don't feel proficient enough with Python yet to implement it myself.

As suggested in the comments, the built-in BioPython pairwise aligner should work for this. You can find a short sequence match in a longer sequence using a global alignment with the right score parameters (query_end_... = 0). The other scoring parameters are values that I've found to be good for aligning similar sequences.
Example code:
from Bio import Align

aligner = Align.PairwiseAligner()
aligner = Align.PairwiseAligner()
aligner.match_score = 5
aligner.mismatch_score = -9
aligner.mode = 'global'
aligner.target_internal_open_gap_score = -12
aligner.target_internal_extend_gap_score = -3
aligner.target_end_open_gap_score = -12
aligner.target_end_extend_gap_score = -3
aligner.query_internal_open_gap_score = -12
aligner.query_internal_extend_gap_score = -3

# Set the query_end penalties to zero to effectively implement a 'glocal' alignment as described in the question
# These are the default values, so you could skip the next two lines,
#    but this is relevant to the question, so I'm including them in the example.
aligner.query_end_open_gap_score = 0
aligner.query_end_extend_gap_score = 0

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATCCG")

for a in alignments:
    print(a)
print()
print()

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATACG")

for a in alignments:
    print(a)
print()
print()

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATACCG")

for a in alignments:
    print(a)


How to convert fasta file to tab delimited file I have a fasta file like

>sample 1 gene 1
atgc
>sample 1 gene 2
atgc
>sample 2 gene 1
atgc

I want to get the following output, with one break between the header and the sequence.

>sample 1 gene 1   atgc
>sample 1 gene 2   atgc
>sample 2 gene 1   atgc

If you have multi-line fasta files, as is very common, you can use these scripts1 to convert between fasta and tbl (sequence_name <TAB> sequence) format:

FastaToTbl

#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
      if (NR>1)
                     printf "\n%s\t", substr($0,2,length($0)-1)
      else
       printf "%s\t", substr($0,2,length($0)-1)
            else
                    printf "%s", $0
}END{printf "\n"}

TblToFasta

#! /usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld

     if (fld == NF-1)
      {
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls)
  {
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}

Save those in your $PATH, make them executable, and you can then do:

$ cat file.fa
>sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
>sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC
$ FastaToTbl file.fa
sequence1   ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGCTAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
sequence2   GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGCATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

And, to get the Fasta back:

$ FastaToTbl file.fa | TblToFasta
>sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
>sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

This can be a very useful trick when searching a fasta file for a string:

TblToFasta file.fa | grep 'foo' | FastaToTbl

If you really want to keep the leading > of the header (which doesn't seem very useful), you could do something like this:

$ perl -0pe 's/\n//g; s/.>/\n>/g; s/$/\n/;' file.fa
>sequence1 ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGCTAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATAC
>sequence2 GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGCATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

But that will read the entire file into memory. If that's an issue, add an empty line between each fasta record, and then use perl's paragraph mode to process each "paragraph" (sequence) at a time:

perl -pe  's/>/\n>/' file.fa | perl -00pe 's/\n//g; s/.>/\n>/g; s/$/\n/;'

1Credit to Josep Abril who wrote these scripts more than a decade ago.


Problems about value log2(IP/Input) less than zero in ChIP-seq? I have a question about the normalization for ChIP-seq. I used CPM to normalize my bam files of each IP and Input. Then I calculate the coverage of gene bodies for all genes on the genome. I have WT and different mutants. Take H3K4me1 as an example, I first calculate the matrix over gene bodies for H3K4me1 IP and Input. Then I used the matrix of H3K4me1 divided by the matrix of input. then log2(H3K4me1/Input) and plot this value over entire gene bodies. I found all the regions are lower than 0 and it seems that my IP does not work. In fact, this is not the case since I have checked the antibody's efficiency via ChIP-qPCR and it works well. Moreover, the distribution of H3K4me1 is consistent with the published data. But the problem is the profile of log2(H3K4me1/Input) is lower than 0.  I don't know how to deal with such problems? I would appreciate it if someone can give me some suggestions on that.  Thanks!

Seconding what Devon says:

Use a proper normalization technique such as TMM from edgeR. For applying this to bigwig files you could follow my tutorial at biostars.org.

My two cents beyond that:

I would not use enrichment over input as input is hard to normalize and can be sparse depending on sequencing depth. Therefore, you might get biased enrichments that do not well reflect reality, especially if counts in IP or input (of both) are low. I personally only use input for peak calling, and for plots like this then compare IPs among each other.

So the strategy would be to first normalize your bigwigs (I guess you use these for the plots?) using the strategy as linked above, and then use these robust CPMs to make the plots without any input samples, only using the IPs. You can use MA-plots to explore whether samples are properly normalized. I made a little code suggestion here in this biostars post on how to do this starting basically from a matrix of raw counts using edgeR. This is basically (from the TMM normalization standpoint) the same as in the first link, but then doing pairwise MA-plots to see whether TMM properly centers the majority of points along y-axis being somewhat zero. Normalization is key in any experiment, therefore I always recommend to use these kinds of plots for diagnostics.


Why are there missing calls in a VCF file from exome sequencing? My data is a VCF file generated from an exome sequencing variant call pipeline. I'm not very familiar with the sequencing and variant calling process. I noticed that there are some missing genotypes, which are recorded as "./." at the GT field. From googling I learned that they're not homozygous reference genotype ("0/0"), but missing calls due to some sequencing failure.

In order to make my data "cleaner", I thought it would be better to filter out the loci with missing calls if there're not too many of them. I also checked the corresponding "DP" of the loci with and without missing calls. For example:

chr1:123 GT:DP 0/0:2 1/1:4 ./.
chr1:234 GT:DP 0/0:10 1/1:11 1/1:20
chr1:345 GT:DP 0/1:40 1/1:37 0/0:78
chr1:456 GT:DP 0/1:7 0/0:23 ./.
chr1:567 GT:DP 0/1:34 1/1:39 0/0:58

In the above toy example, there are 3 people and 5 loci. I checked the mean DP of all 5 loci and I found that the mean DP of loci with missing calls (1st and 4th loci) to be significantly lower than the mean DP of loci without missing calls (2nd, 3rd and 5th loci). This happens to my real data with a lot more loci and samples than this toy example. Is that a coincidence? Or are there any specific reason that the loci with missing calls have lower coverage than the "normal" loci? Thanks!

Missing variant calls due to lack of coverage shouldn't happen in the targeted capture region and I'd think most of these would come from off-target regions where some samples had reads mapped. I'd filter out the VCF to only include on-target loci before proceeding with further analyses.

If you did your variant calling on samples separately and then merged them together, you'll have missing calls at loci called variant in one sample and homozygous reference in the other.

You can get over this by asking the variant caller to emit all sites within the target region, which would lead to a large VCF or do multi-sample calling, which is more computationally intensive.

One other way is to do variant calling on all samples separately, generating a BED file of all variant regions in the cohort and then asking the variant caller to genotype those loci across all samples before merging them together.


Annotation format design Bashing file formats is a favorite pastime in bioinformatics, and annotation file formats such as GFF and BED seem to get special attention. A lot of this frustration stems from community's shockingly inconsistent adherence to specifications and conventions, but there are also some (dare I say objectively) problematic design choices in each of these formats.

GFF (and its more common derivatives GTF and GFF3) use 1-based closed interval notation, which optimizes for human comprehension but is far inferior to 0-based half-open interval notation (such as used by BED) for computations involving interval arithmetic.
Although BED and GTF were designed for very specific use cases (visualization and gene prediction, respectively), they have been used and abused in a much wider set of contexts. For example, the BED fields related to the thick part are irrelevant if you're not plotting them in a genome browser.
BED supports a single level of feature decomposition (a feature can be broken up into blocks). GTF supports two levels (exons grouped by transcript_id, transcripts grouped by gene_id). In contrast, GFF3 supports an arbitrary number of levels, and uses parent/child relationships defined by ID and Parent attributes to declare a directed acyclic graph of features.
Data that does not fit into mandatory pre-defined fields must be relegated to optional fields or free-form attribute key/value pairs. While this flexibility is powerful, a common complaint is that "all the action" happens in these optional/free-form fields.
There is a dearth of validation tools, and those that do exist focus primarily on validating syntax and not semantics. To use an aging analogy, it's one thing to say an XML file is valid, but it's completely different to validate it against a schema. There are essentially no widely used tools that do the latter for annotation files.

If we were tasked with creating a new annotation format, and if we were guaranteed the resources needed to develop it, and interest and wide adoption from the wider community (one can dream!), what design criteria should be considered in the development of this new format? What, if anything, makes an objectively good annotation data format?

Presuming we consider "human readable", "easily parsable", and "quickly queryable" to be objectively good qualities (and if not, I worry for the future):

Text-based: It's absurdly common to want to use grep or awk on annotations. Sure, one could make variants of these that are binary-format aware, but why reinvent the wheel. Of course text files don't innately allow region-based queries of their contents, so on to point 2...

This should further explicitly be line-based. Fields would need to be tab-separated (ideally it'd use the ascii record separator character instead, but I fear that ship has long ago sailed).

A strictly defined line order: This is one of my personal pet-peeves about BED/GFF/GTF files. If you're going to make a text-based annotation format, everyone would end up ahead if said format made explicit that it should be sorted. This then allows things like tabix to be used so the "query a region" problem is then solved. But I would go further than that. The issue with things like GFF is that there are multiple inter-dependent lines and there's no strict rule about whether a parent line absolutely must come before a child line. This just makes implementing things a nightmare and more often than not a randomly ordered file will just break tools. Since GTF or GFF-style interdependent lines likely be how any annotation format works, this ordering between lines with the same start position should be made explicit in the format.


tp53 example in Momand's book I am new to bioinformatics. While reading Momand's 'Concepts in Bioinformatics and Genomics' , I understood clearly that each codon translates to a specific amino acid.

But then, he talks about the TP53 gene, and then without explanation talks about P53... Is that a naming convention? It's not explicitly explained... Does that applies for any case? ie TP90 gene encodes protein P90?

I don't understand the naming: correct me if wrong, but they named the gene after the protein's weight, which is 53 kilodaltons?
(It confuses me with the location of the gene).

Lastly I do not understand: if there is a mutation in the DNA, the mutation shall appear in only one cell I guess, how then the others cells' DNA start to mutate? or is it that all DNA in all cells are mutated instantly?

This is possibly more of a molecular genetics question appropriate to biology.stackexchange.com, also note that it's usually better on stackexchange sites to restrict your post to a single question and do some preliminary research via google or your search engine of choice. But I don't see any harm in answering it here nonetheless.

The real answer is that naming conventions for genes/proteins are incredibly inconsistent and counterintuitive. You can't really draw any conclusions about a gene/protein based on its name alone. For example, some proteins such as p53 are named for biochemical features, others are named for phenotypes observed when their gene is mutated, such as for instance the gene decapentaplegic in fruit flies. Often they involve little jokes, such as the gene called mothers against decapentaplegic which represses the decapentaplegic mutant phenotype. So it's all a bit of a mess and people kind of assume that you'll follow along, which can be confusing.

It's a little rarer that gene names differ substantially from the names of the proteins they encode, but not unheard-of. This is indeed the case with your TP53/p53 example. See a short explainer here or a longer one here. What you'll notice is that there are in fact multiple different names for the same gene/protein. What you'll often see is that the shorter name (p53) is the one that people tend to use colloquially.

For your second question regarding how mutations are propagated, the short answer is that mutations are inherited from mother to daughter cells. Cells don't generally "infect" other cells with mutations (though in biology there are always exceptions). I think that what you are interested in is the difference between "somatic" and "germ-line" mutation.

Hope that helps.


What are RIKEN genes? I am using data from the Tabula Muris Consortium.

These are gene counts from scRNA-seq of mouse cells.

There are some specific genes with name ending with suffix 'Rik' (e.g. 0610005C13Rik or 0610007C21Rik). Searching online I found that they are from RIKEN project.

However, they are different in some respects (e.g. antisense lncRNA gene, or mouse 0610007C21RIK Adenovirus).

I want to exclude them from downstream analysis, but it is unclear to me the main purpose of this collection. For which kind of analysis are they used?

My understanding is that these are cDNA sequences from the RIKEN project that don't yet have another gene name. The project carried out a comprehensive sequencing of the mouse genome, and as a result found a lot of cDNA transcripts that were previously unknown. Now they have some annotation (which can be varying degrees of functional annotation), and an accepted official identifier to allow people to talk about these transcripts when discovered in other studies.

It'd be interesting to see if any novel transcripts have been found from the Tabula Muris project. They would also need their own identifier, possibly leading to similar intrigue about their history.


Problem installing last build of Sniffles I have followed the steps described here to install Sniffles on my Mac Book Pro running macOS High Sierra.

However, I am struggling with this step:

cmake -D CMAKE_C_COMPILER=/opt/local/bin/gcc-mp-4.7 -D CMAKE_CXX_COMPILER=/opt/local/bin/g++-mp-4.7 ..

I do not have an /opt folder, so I get the following output:

-- The C compiler identification is unknown
-- The CXX compiler identification is unknown
CMake Error at CMakeLists.txt:2 (project):
  The CMAKE_C_COMPILER:

    /opt/local/bin/gcc-mp-4.7

  is not a full path to an existing compiler tool.

  Tell CMake where to find the compiler by setting either the environment
  variable "CC" or the CMake cache entry CMAKE_C_COMPILER to the full path to
  the compiler, or to the compiler name if it is in the PATH.

CMake Error at CMakeLists.txt:2 (project):
  The CMAKE_CXX_COMPILER:

    /opt/local/bin/g++-mp-4.7

  is not a full path to an existing compiler tool.

  Tell CMake where to find the compiler by setting either the environment
  variable "CXX" or the CMake cache entry CMAKE_CXX_COMPILER to the full path
  to the compiler, or to the compiler name if it is in the PATH.

-- Configuring incomplete, errors occurred!
See also "/Users/xxx/Documents/xxx/xxx/Sniffles-master/build/CMakeFiles/CMakeOutput.log".
See also "/Users/xxx/Documents/xxx/xxx/Sniffles-master/build/CMakeFiles/CMakeError.log".

Unsurprisingly, the error file says that /opt/local/bin/gcc-mp-4.7 does not exist.

So I tried:

cmake -D CMAKE_C_COMPILER=/usr/bin/gcc -D CMAKE_CXX_COMPILER=/usr/bin/gcc ..

I get another long message but the most important bit is:

CMake Error at /usr/local/Cellar/cmake/3.8.1/share/cmake/Modules/FindPackageHandleStandardArgs.cmake:137 (message):
  Could NOT find OpenMP (missing: OpenMP_C_FLAGS OpenMP_CXX_FLAGS)

I do not know what to do to solve this error.

I do have gcc here: /usr/bin/gcc and cmake version 3.8.1.

I have created an issue on the Sniffles GitHub repository.

Update 13th June 2018

I installed the latest version of gcc:

brew install gcc
Warning: gcc 8.1.0 is already installed and up-to-date
To reinstall 8.1.0, run `brew reinstall gcc`

Then, I tried to make cmake use gcc and g++ in the /usr/local/bin folder.

cmake -D CMAKE_C_COMPILER=/usr/local/bin/gcc-8 -D CMAKE_CXX_COMPILER=/usr/local/bin/g++-8 ..

This recreated the same problem:

CMake Error at /usr/local/Cellar/cmake/3.8.1/share/cmake/Modules/FindPackageHandleStandardArgs.cmake:137 (message):
  Could NOT find OpenMP (missing: OpenMP_C_FLAGS OpenMP_CXX_FLAGS)

Do I need to set the OpenMP_C*_FLAGS to something?

You need to install GCC, or a newer version of Clang.

Your compiler — though the binary is called gcc — is actually Clang, and the Clang version that’s shipping with macOS/Xcode still does not support OpenMP.

The easiest course of action is probably to install GCC via Homebrew (brew install gcc) and then configure Sniffles to use Homebrew’s GCC (i.e. setting the CMAKE_C*_COMPILER paths to /usr/local/bin/g{cc,++}-$version).


Odds ratio calculations in GWAS paper Posting here because the reddit thread I made didn't gain any traction.

I'm having a little bit of trouble understanding how the odds ratio was calculated in this table, specifically the OR marked in red:

To my understanding, an odds ratio = $$[(\text{# of cases with exposure}) / (\text{# of controls with exposure})] / [(\text{# of cases without exposure}) / (\text{# of controls without exposure})]$$

In the red-highlighted case, my odds ratio is

$$(0.661 / 0.6580) / [(1-0.661) / (1-0.6580)] = 1.0134$$

which is very off from $0.98$. Also, I know I'm using frequencies rather than absolute counts, but it shouldn't matter because if I were to multiply each frequency by the corresponding population values, all the population values would cancel out.

So, what have I done wrong?

Link to the paper

Those odds ratio aren't calculated from the raw allele frequencies. They're calculated from coefficients estimated from complex linear models with a bunch of covariates that are run on subsets of the data then meta-analysed together.

It is a bit surprising that the odds ratio is less than one when the frequency is higher in cases, but the p-value indicates that the OR is not significantly different from 1 so it may be down to statistical noise. It's possible that within populations one allele is more common in cases than in controls but that cases are primarily sampled from populations where that allele is uncommon, causing it to be more common in controls overall.  See here. Alternatively, 1/0.98 is pretty close to the value you calculated, so it's possible that the reference allele got switched around at some point.


What does final generated models with TM-score of DMPfold signify in "sequence.ema.txt" file? So, my runs have finished and got 5 final models. I predicted TM-score of one model to understand the result. The output of "predict-tmscore.sh" for "final_1.pdb" is:

1 models read from PDB file
1 models read from PDB file
nres is 377
neff is 15423
dpscore is -300.599030
dpscorenorm is -0.082424
Probability of each TM-score bin (cols 1-10), estimated TM-score (col 11) and prob TM > 0.5 (col 12) written to sequence.ema.txt
Estimated TM-score of ../final_1.pdb is 0.5

could you please help me understand what does "nres", "neff", "dpscore", and "dpscorenorm" below mean:

    nres is 377
neff is 15423
dpscore is -300.599030
dpscorenorm is    -0.082424

and it also generates a "sequence.ema.txt" file with below value:

2.0308916646172293e-05 2.0682286958617624e-06 0.04564645141363144 0.14857041835784912 0.27177685499191284 0.31814685463905334 0.19150513410568237 0.022412151098251343 0.001872337656095624 4.743059980683029e-05 0.5036169363265117 0.5339840054512024

Please let me know what these above values signify?

And,I would really appreciate your comments on what does "rawdistpred.1" file mean? There are total 9 such files. When predicitng TM-score, do I need to use "rawdistpred.1" file for all final models or I can use any "rawdistpred.*" file?

thanks!

This is described briefly in the paper. The four quantities are inputs to the estimator of model accuracy:

nres is the sequence length.
neff is an effective sequence count, or how deep the alignment was.
dpscore is how well the final model matches the initial distance prediction.
dpscorenorm is a normalised version of the above.

Don't worry too much about directly interpreting these. They are only given for completion. The important line is Estimated TM-score of ../final_1.pdb is 0.5, indicating that DMPfold just about thinks this is a decent model.

The "sequence.ema.txt" is more raw output giving the probability of the TM-score being in each bin. "rawdistpred.1" is the distance predictions from the first iteration. Always use "rawdistpred.1" as input to the TM-score predictor as this is how it was trained.


How to determine a correct promotor sequence? (Bisulfite sequencing - uni assingment) aspiring bioinformatician here, I'd like to ask you for help with an university assignment. First let me post the entire question:
Design primers for bisulfite sequencing of the promoter of the ACE2 gene in microbat  Myotis lucifugus (based on the annotation of the Broad Institute Myoluc2.0/myoLuc2). The gene is localized on the contig GL429816. As a transcription start, use the start of the ACE2 transcript variant XM_023753670.1; once you find ACE2 in the UCSC genome browser .You will find it in NCBI RefSeq track
I don't want anyone to solve it for me but I do have a problem with one aspect. I am not clear on how to determine what is a valid promotor for this gene in microbat. I do know it is the upstream sequence however I am not clear on how to tell the number of upstream bases that are the actual promotor. I did find some sources for the actual promotor sequence determined by much smarter people than myself however those are for ACE in humans so I am afraid that it is not applicable in my case.
Any help will be appreciated, even just pointing me in the right direction where I could find my answer would be enough. (And yes I tried contacting my teacher about this, but unfortunately he is not answering and the deadline is approaching rather quickly).

I think there is no fixed size. Ideally you would have something like histone ChIP-seq or ATAC-seq data to approximate the relevant upstream interval. Since this is not always possible most people go with -500bp. Others might have different suggestions though. Maybe matching known motifs of transcription factors that regulate that promoter to the upstream region can help approximating a better window. But as Sanger sequencing usually gives about 1kb of sequence I'd just take a 500bp to 1kb window and sequence that, fwd and rev primers for sequencing to have double confirmation.


How to calculate fold change in gene expression from RSEM (RNAseq) I have gene expression data from RNAseq, specifically: log2(x+1) transformed RSEM normalized counts.

How can I convert these data from multiple samples to determine fold change in gene expression?

You should use a proper statistical framework for RNA-seq dfferential analysis (which includes FC calculation). Standard tools for this are (among others) edgeR or DESeq2. You could use tximport to import RSEM outputs into R and then use its output for e.g. DESeq2. The linked manual provides example code for this. Please read the manuals of edgeR or DESeq2 towards the respective downstream analysis. Both are extensive and compreshensive.


Interpreting this PCA plot for RNA-seq I have RNA-seq from two sequencing batches; Lab technician says that he has run the RNA expression quantification two times in bathes 1 and 2 for example tumor 1 in batch 1 and tumor 1 in batch 2 , normal 2 in batch1 and normal 2 in batch 2. This is my design for DESeq2

> head(mycols)
        condition batch
N_1_305         N     1
N_1_310         N     1
N_1_337         N     1
N_1_353         N     1
T_1_305         T     1
T_1_310         T     1
> tail(mycols)
        condition batch
T_2_337         T     2
T_2_338         T     2
T_2_344         T     2
T_2_346         T     2
T_2_349         T     2
T_2_353         T     2
>

I got this PCA plot

And this is biplot of samples

In PCA plot I am seeing for instance
, T_1_337 (batch1) has been placed too close to T_2_337 (batch2)

Then I used svd for detecting hidden batch

Does this mean that there is no big batch effect between experimental runs and I can concatenate the fastqs from both batches for each sample (technical replication) or collapse technical replicates afterwards?

Please help me to interpret these
Thank you

  EDITED

Sorry @Devon, I have 4 lanes for each samples (paired end) in each of experimental runs; For concatenating fastq files can I do like this ?

cat fastq1_lane1_batch1 fastq1_lane1_batch2 fastq1_lane2_batch1   fastq1_lane2_batch2  fastq1_lane3_batch1 fastq1_lane3_batch2 fastq1_lane4_batch1 fastq1_lane4_batch2  > fastq1

cat fastq2_lane1_batch1 fastq2_lane1_batch2 fastq2_lane2_batch1   fastq2_lane2_batch2  fastq2_lane3_batch1 fastq2_lane3_batch2 fastq2_lane4_batch1 fastq2_lane4_batch2 > fastq2

Yes, you can safely concatenate the technical replicates. Odds are good that these are even the same libraries just sequenced twice, so even labeling them as replicates is a bit of a stretch. As an aside, it would be surprising if you actually had a batch effect in a situation like this. You will commonly see sequencing facilities just sequence a given sample a second time if it didn't get the number of reads requested by the client to begin with. They then just have the client concatenate the files.


Using variables with fasterq-dump? I am trying to download multiple fastq files from the SRA NCBI database which, conveniently enough, have their IDs in a range. Building a simple loop doesn't seem to do the trick, as bash appears to interpret my variable as an additional command? Help would be appreciated. This is what I have for now:
for i in {256..345};do
    sra_code = "DRR163""$i"
fasterq-dump "$sra_code" -p -v --outdir ./output/ -e 8

done
```

You have an obvious error in your shell script that should have given you an error message:
sra_code = "DRR163""$i"

Should be:
sra_code="DRR163""$i"

You cannot have spaces around the = sign in a variable assignment. The line sra_code = "DRR163""$i" means "run the command sra_code with the arguments = "DRR163""$i". Since there is no such command, you see an error.
Here's a working version of your command:
for i in {256..345}; do
    sra_code="DRR163$i"
    fasterq-dump "$sra_code" -p -v --outdir ./output/ -e 8
done


I assume by "reverse sort of variance" you mean "highest variance". Assuming you made a matrix out of that (set the row names to the first column and then remove it) and called it m:

sel = order(apply(m, 1, var), decreasing=TRUE)[1:100]

sel then contains the indices into your matrix (or the original dataframe).

BTW, I hope you're not filtering by variance before performing differential expression. That would be incredibly incorrect. Doing this sort of filtering for performing PCA or something like that can be useful, though you should really be using normalized counts (or a transformation of them).


What do the read colors in IGV mean? I was looking at a bam file in the IGV viewer and saw:

What do the different read colors mean? Why is one read a light blue color, another green, another aquamarine (?), another purple and another two orange? What do these different read colors represent? The IGV docs have sections on interpreting color by insert size, pair orientation, bisulphide mode, splice junctions and sashimi plot, but none of them seems to be describing what I'm looking at there.

More generally, where are all these default behaviors documented? Are they somewhere in http://software.broadinstitute.org/software/igv/UserGuide and I'm missing them?

Could it be that these are reads with their mates found on other chromosomes? That's the part of the IGV documentation which seems to make sense in this context:

for paired end reads that are coded by the chromosome on which their mates can be found


Way to get genomic sequences at given coordinates without downloading fasta files of whole chromosomes/genomes first? So I have a list of start and stop positions along chromosomes in different species, and I'd like to get the corresponding DNA sequence for each set of coordinates. In the past, I've just download the genome as a fasta file and then use pyfaidx to extract the sequences at the given positions. But now that I'm working with several species at once, I was wondering if there's any kind of tool in Python or R that can fetch your sequences of interest without downloading a bunch of large files. Thanks

using a http request.

if there is a DAS server, you can always use this protocol to download the xml -> fasta. see https://www.biostars.org/p/56/

$ curl -s "http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chrM:100,200"  | xmllint --xpath '/DASDNA/SEQUENCE/DNA/text()' - | tr -d '\n'
ggagccggagcaccctatgtcgcagtatctgtctttgattcctgcctcattctattatttatcgcacctacgttcaatattacaggcgaacatacctacta

or use the UCSC utility twoBitToFa which works with remote files.

If you have an NCBI accession number you can use the E-utilities with seq_start and seq_end.

$ wget -q  -O -  "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=AE014134.1&seq_start=100&seq_stop=300&rettype=fasta"
>AE014134.1:100-300 Drosophila melanogaster chromosome 2L complete sequence
TGCCAACATATTGTGCTCTTTGATTTTTTGGCAACCCAAAATGGTGGCGGATGAACGAGATGATAATATA
TTCAAGTTGCCGCTAATCAGAAATAAATTCATTGCAACGTTAAATACAGCACAATATATGATCGCGTATG
CGAGAGTAGTGCCAACATATTGTGCTAATGAGTGCCTCTCGTTCTCTGTCTTATATTACCG


Basic questions about GSEA I'm currently learning about Gene Set Enrichment Analysis (GSEA) in the hopes of using it in my analysis of differentially expressed genes, and I just had a few questions about the program, specifically about GSEAPreranked, which I need cleared up.

On the ranked list needed for GSEA input, should the list include all genes, or only those that pass a certain threshold of significance (i.e. fold change higher than 2, p value less than 0.05, etc.)? Ideally I'd like to sort the genes by fold change alone as I don't trust my p values as much, so should I only include genes with high fold changes?
I am comparing multiple conditions of disease with different treatments. Am I correct that GSEA only compares two conditions? If this is the case should I run GSEA for each control/treatment comparison? Would this be conventional?

You seem to refer to the GSEA provided by the Broad institute, (there are other GSEA algorithms).

1) You can provide whatever you wish, but if you want to know if those gene sets in which side of the ordered list are they, then provide all the list (of genes) you have.

2) GSEA analyize if the order of a given list distributes in certain way the elements of another (unordered) list. In your case you would need to compute the GSEA for each treatment vs control comparison.


How to detect mismatch before mapping in RNA-Seq data In the methods of this paper, the authors say:

  Reads were then filtered based on quality score in the UMI region. Any read with >1 low-quality base (phred <=10) were discarded. Reads with more than one mismatch in any of the three 8 nt cell barcodes were also discarded. The cDNA reads (Read 1) were then mapped to either a combined mm10- hg19 genome or the mm10 genome using STAR (44).

I understand how to filter based on mismatch after mapping but in this case it is done before mapping. Has this to be interpreted as "Reads with low quality in more than one base in any of the three 8 nt cell barcodes were also discarded"?

At each step in the SLiTing and pooling the barcodes added are not random, but are taken  from a predefined list of possible barcodes.

When the reads arrive, first you filter by reads that have >1 low quality bases in the UMI region (which is truely random).

Then you compare the sequences in the 3 cell-barcode regions against the list of allowed barcodes for each of those positions. If a read is not within an edit distance of one from one of those barcodes, it is discarded.

Finally the reads that pass these filters can be aligned.


How can I read FCS files using open source libraries? FCS is a patented data format used for storing flow cytometry data. The most recent version is FCS3.1. There is some documentation on the format, but there is no information on how to read these files. There are some R packages and a MATLAB code to read an FCS file, but I am looking for standard libraries developed either by the FCS consortium or any other group. I also wish to know if FCS is a subset of an existing standard data format that can be read by a standard library using any programming language.

Finally, I would want to convert these files to an easily readable format like HDF5.

A few years ago, I wrote a python script to convert FCS files into tab-separated format. It was far from handling all the possibilities that the format description offers, but at least it worked for some of the files produced on one of our machine: http://www.igh.cnrs.fr/equip/Seitz/en_equipe-programmes.html

The format documentation I found enabled decoding (see section 3 of the pdf you mention), but it requires reading data in binary mode.

The general idea of this format (and, I guess, many other binary formats), is that there is a header zone at the beginning of the file with a defined number of fields encoding numbers indicating how the rest of the file is structured. So a first phase is to parse this header, following the description given in the documentation of the format. The information extracted from the header tells where to find the data and how it is encoded, still according to rules described in the format documentation.

In case this may be useful, and for the record, here is the code from the above-mentioned script (after stripping comments, some of which are merely copied from the format documentation, and adding a few ones):

#!/usr/bin/env python
"""This script tries to read FCS flow cytometry data.
Format parsing inspired by information found here:
http://isac-net.org/Resources-for-Cytometrists/Data-Standards/Data-File-Standards/Flow-Cytometry-Data-File-Format-Standards.aspx
"""

import re
# To decode binary-encoded data
import struct
import sys

class Parameter(object):
    """This object represents one of the parameter types that are present in a DATA segment of a FCS file."""
    __slots__ = ("p_name", "p_bits", "p_range", "p_ampl", "parser")
    def __init__(self, p_name, p_bits, p_range, p_ampl):
        self.p_name = p_name
        self.p_bits = p_bits
        self.p_range = p_range
        self.p_ampl = p_ampl
        # Function for parsing a value of the parameter in the data segment
        self.parser = None

##############################################
# Here starts the parsing of the header part #
# which tells where the other parts are.     #
##############################################

f = open(sys.argv[1], "rb")
# The format name is encoded in 6 letters
# An ASCII letter is coded with one octet
file_format = "".join([f.read(1) for __ in range(6)])
sys.stdout.write("Format: %s\n" % file_format)
# The format descriptions reserves 4 octets that we skip
skip = f.read(4)
# 8 octet chunks encode the start and end positions
# of different parts of the data
text_start = int(f.read(8).strip(" "))
text_end = int(f.read(8).strip(" "))
data_start = int(f.read(8).strip(" "))
data_end = int(f.read(8).strip(" "))
analysis_start = int(f.read(8).strip(" "))
analysis_end = int(f.read(8).strip(" "))

if (analysis_start and analysis_end):
    sys.stderr.write("Cannot deal with ANALYSIS segment of an FCS file.\n")

####################################################
# Here starts the parsing of the "TEXT" portion    #
# which describes how the data proper is organized #
####################################################
f.seek(text_start)
# The first character in the primary TEXT segment is the ASCII delimiter character.
sep = f.read(1)
if sep not in ["_", "@"]:
    alt_sep = "_@_"
elif sep not in ["_", "|"]:
    alt_sep = "_|_"
else:
    assert sep not in ["+", "|"]
    alt_sep = "+|+"
text_segment = f.read(text_end - text_start)

fields = text_segment.split(sep)

info = {}

i = 0
while i < len(fields) - 1:
    key = fields[i]
    i += 1
    val = fields[i]
    i += 1
    # Keywords are case insensitive, they may be written in a file in lower case, upper case, or a
    # mixture of the two. However, an FCS file reader must ignore keyword case. A keyword value may
    # be in lower case, upper case or a mixture of the two. Keyword values are case sensitive.
    info[key.upper()] = val
print "%s events were detected." % info["$TOT"]
print "Each event is characterized by %s parameters" % info["$PAR"]

if info["$NEXTDATA"] != "0":
    sys.stderr.write("Some other data exist in the file but hasn't been parsed.\n")

# L - List mode. For each event, the value of each parameter is stored in the order in which the
# parameters are described. The number of bits reserved for parameter 1 is described using the
# $P1B keyword. There can be only one set of list mode data per data set. The $DATATYPE
# keyword describes the data format. This is the most versatile mode for the storage of flow
# cytometry data because mode C and mode U data can be created from mode L data.
assert info["$MODE"] == "L"

parameters = []

# indices of the parameters
p_indices = range(1, int(info["$PAR"]) + 1)
for i in p_indices:
    p_name = info["$P%dN" % i]
    p_bits =  info["$P%dB" % i]
    p_range = info["$P%dR" % i]
    p_ampl =  info["$P%dE" % i]
    parameters.append(Parameter(p_name, p_bits, p_range, p_ampl))

sys.stdout.write("The parameters are:\n%s\n" % "\t".join([par.p_name for par in parameters]))

# How are 32 bit words organized
if info["$BYTEORD"] == "4,3,2,1":
    endianness = ">"
else:
    endianness = "<"
    assert info["$BYTEORD"] == "1,2,3,4"

    # I stripped a long comment which is just a copy of the documentation
# Type of data:
if info["$DATATYPE"] == "I":
    for par in parameters:
        nb_bits = int(par.p_bits)
        assert nb_bits % 8 == 0
        nb_bytes = nb_bits / 8
        # Determine format string for unpacking (see https://docs.python.org/2/library/struct.html)
        if nb_bytes == 1:
            c_type = "B" # unsigned char
        elif nb_bytes == 2:
            c_type = "H" # unsigned short
        elif nb_bytes == 4:
            c_type = "L" # unsigned long
        elif nb_bytes == 8:
            c_type = "Q" # unsigned long long
        else:
            raise ValueError, "Number of bytes (%d) not valid for an integer (see https://docs.python.org/2/library/struct.html#byte-order-size-and-alignment)." % nb_bytes
        fmt = "%s%s" % (endianness, c_type)
        p_range = int(par.p_range)
        def parser(data):
            value = struct.unpack(fmt, data.read(nb_bytes))[0]
            try:
                assert value < p_range
            except AssertionError:
                print "Value %s higher than %d" % (str(value), p_range)
            return value
        par.parser = parser
    pass
else:
    raise NotImplementedError, "Only the parsing of integer value has been implemented so far."

out_file = open(sys.argv[2], "w")
out_file.write("#amplification_types\t" + "\t".join([par.p_ampl for par in parameters]) + "\n")
out_file.write("parameters\t" + "\t".join([par.p_name for par in parameters]) + "\n")
i = 1
##############################################
# Here starts the parsing of the data proper #
##############################################
f.seek(data_start)
while f.tell() < data_end:
    values = []
    for par in parameters:
        values.append(par.parser(f))
    out_file.write("%d\t" % i + "\t".join(map(str, values)) + "\n")
    i += 1
out_file.close()
f.close()


Filtering the expression values of a gene from GEO I'm filtering the gene expression value of a particular gene using Bioconductor in R. Following the instruction given here, I could filter the gene expression value for a gene from the dataset GDS3782.

The output is as follows for the gene of my interest,

GSM524151  "23.82"
GSM524152  "27.49"
GSM524153  "24.26"
GSM524154  "17.98"
GSM524155  "25.87"
GSM524156  "26.16"
GSM524157  "26.02"
GSM524158  "25"
GSM524159  "29.85"
GSM524160  "18.25"
GSM524161  "28.86"
GSM524162  "27.33"
GSM524163  "27.42"
GSM524164  "28.23"
GSM524165  "22.09"
GSM524166  "24.85"
GSM524167  "27.11"
GSM524168  "32.01"
GSM524169  "28.97"
GSM524170  "30.98"

The above data corresponds to disease and control which I get using

Meta(gds3782)$description

I wish to see the description of each value. For example, from the webpage, I have seen that the first ten values correspond to "disease-state" and the last ten correspond to "control". How do I view this classification while processing the data using Bioconductor in R?

You're looking for the Columns() function:

> Columns(gds3782)
     sample        disease.state           description
1  GSM524151 non-diabetic control  Value for GSM524151: beta-cells_non-diabetic condition_donor1; src: human beta-cells_non-diabetic
2  GSM524152 non-diabetic control  Value for GSM524152: beta-cells_non-diabetic condition_donor2; src: human beta-cells_non-diabetic

Note that it's probably easier to work with a more common structure like and expressionSet:

eset = GDS2eSet(gds3782)

You'd then use the pData() function to return the same thing.


How can I specify an email address in the SLURM snakemake json configuration file? I am working on a SLURM-based cluster. I am currently using the following slurm json file:

{
    "__default__" :
    {
        "A" : "overall",
        "time" : "96:00:00",
        "nodes": 1,
        "ntasks": 1,
        "cpus" : 24,
        "p" : "LONG",
        "mem": 16G,
        "output": "snakemake%A.out",
        "error": "snakemake%A.err"
    }
}

along with the following command:

snakemake --jobs 200 --use-conda --printshellcmds  --cluster-config slurm.json --cluster "sbatch -A {cluster.A} -p {cluster.p} -t {cluster.time} --output {cluster.output} --error {cluster.error} --ntasks {cluster.ntasks} --cpus-per-task {cluster.cpus} --mem {cluster.mem}"

I would like to be able to specify my email and receive an email when a job finishes. I have tried adding a line in the json file with "email": my@email.com and an option to the sbatch command --email {cluster.email} but that did not work.

Anyone knows how to do it?

As noted by the comments and the SLURM manual (not so much bioinformatics related?)

You need the following terms in your batch/queuing command:

--mail-type=ALL
--mail-user=you@email.com

I used "ALL" which will send you all reports. If you just want the finish message then you would use "END" and "FAIL" like such:

--mail-type=END,FAIL

Here are the specifics from the manual if you wish to make a different selection:

--mail-type=

  Notify user by email when certain event types occur. Valid type values
  are NONE, BEGIN, END, FAIL, REQUEUE, ALL (equivalent to BEGIN, END,
  FAIL, REQUEUE, and STAGE_OUT), STAGE_OUT (burst buffer stage out and
  teardown completed), TIME_LIMIT, TIME_LIMIT_90 (reached 90 percent of
  time limit), TIME_LIMIT_80 (reached 80 percent of time limit),
  TIME_LIMIT_50 (reached 50 percent of time limit) and ARRAY_TASKS (send
  emails for each array task). Multiple type values may be specified in
  a comma separated list. The user to be notified is indicated with
  --mail-user. Unless the ARRAY_TASKS option is specified, mail notifications on job BEGIN, END and FAIL apply to a job array as a
  whole rather than generating individual email messages for each task
  in the job array.

--mail-user=

  User to receive email notification of state changes as defined by
  --mail-type. The default value is the submitting user.


How to quickly determine mutations in a read of a sam file? After DNA sequencing, I generated a sam file through alignment of a fastq file.
Before using well known variant calling programs (eg. Annovar etc.), I want to pick some reads and know what kinds of mutations are there.

Is there any quick way of doing it?

For qualitative analysis, you're probably better off using something less granular like IGV or IGB. However, if you really want to look at a couple of reads:

If you're willing to ignore sequencing errors, you can inspect the CIGAR string or the MD tag, both of which give information on the alignment of a single read.

The CIGAR string gives details on insertions, deletions, clippings, matches and mismatches. From Genome Analysis Wiki,

  The sequence being aligned to a reference may have additional bases that are not in the reference or may be missing bases that are in the reference. The CIGAR string is a sequence of of base lengths and the associated operation. They are used to indicate things like which bases align (either a match/mismatch) with the reference, are deleted from the reference, and are insertions that are not in the reference. For example:

RefPos:     1  2  3  4  5  6  7     8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T     G  A  A  C  T  G  A  C  T  A  A  C
Read:                   A  C  T  A  G  A  A     T  G  G  C  T
With the alignment above, you get:
POS: 5
CIGAR: 3M1I3M1D5M

Most common usage of the CIGAR string uses M (match/mismatch), I (insertion), D (deletion), S (soft clipping), and H (hard clipping). Note that = (match) and X (mismatch) are available as alternatives to the less informative M, but they are less widely used.

The MD tag gives specific details on mismatches and deletions. From the SAMtools tags specification,

  The MD field aims to achieve SNP/indel calling without looking at the reference. For example, a string ‘10A5^AC6’ means from the leftmost reference base in the alignment, there are 10 matches followed by an A on the reference which is different from the aligned read base; the next 5 reference bases are matches followed by a 2bp deletion from the reference; the deleted sequence is AC; the last 6 bases are matches. The MD field ought to match the CIGAR string.

Note that neither of these will give you any idea of structural variants in short reads, and neither will be particularly readable (or helpful, due to the higher error rate) in long reads.


Understanding the significance of BLAT score calculations for read fragment alignment I'm attempting to reconstruct the read fragments that form a hit to a target, excluding/skipping gaps. BLAST's outfmt 6, and alternatives such as Diamond's M8 format only give the number of gap opens on the alignment. One can't just use the start and end positions of the alignment as indices to slice the read DNA sequence, as this will include the nucleotides that were gapped by the alignment.

PSL appeared to be the only format I could find that gave the actual start positions (and lengths) of each of the pieces that formed a hit (that is, each gap forces a new block).

The regular BLAST output is human readable and is a nightmare to parse. The XML format mentioned has a structure for each HSP, but an HSP can still include gaps:

  [...]
  <Hit_hsps>
    <Hsp>
      <Hsp_num>1</Hsp_num>
      <Hsp_bit-score>51.6</Hsp_bit-score>
      <Hsp_score>122</Hsp_score>
      <Hsp_evalue>9.8e-06</Hsp_evalue>
      <Hsp_query-from>2</Hsp_query-from>
      <Hsp_query-to>98</Hsp_query-to>
      <Hsp_hit-from>35</Hsp_hit-from>
      <Hsp_hit-to>65</Hsp_hit-to>
      <Hsp_query-frame>2</Hsp_query-frame>
      <Hsp_hit-frame>0</Hsp_hit-frame>
      <Hsp_identity>24</Hsp_identity>
      <Hsp_positive>27</Hsp_positive>
      <Hsp_gaps>1</Hsp_gaps>
      <Hsp_align-len>32</Hsp_align-len>
         <Hsp_qseq>ITAIGAGLQGPAGCEVIDAGGLLVMPGGIDTH</Hsp_qseq>
         <Hsp_hseq>IAAVGTGLE-PAGAEIIDAGGLLVMPGGIDVH</Hsp_hseq>
      <Hsp_midline>I A+G GL+ PAG E+IDAGGLLVMPGGID H</Hsp_midline>
    </Hsp>
  </Hit_hsps>
</Hit>

I'm trying to be lazy, I don't want to have to parse the Hsp_qseq and Hsp_hseq to determine exactly where my gaps are. I just want co-ordinates of where my read hits a protein. The PSL format (below) tells me this (tStarts):

match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block   blockSizes      qStarts  tStarts
        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
[...]
21      4       0       0       1       8       1       24      ++      <PROTEIN>   662     321     354     <READ>      101     2       101     2       8,17,   321,337,        2,50,

As an alternative, I've executed the BLAT aligner as a means to find alignments of translated DNA reads against a small protein database. I selected BLAT as it provided the easiest means to access information on the actual alignment blocks (with the smallest amount of parsing). That is, the PSL has a comma delimited field that enumerates each of the starting positions of a query to a target.

However, BLAT outputs a significantly greater number of alignments than alternative tools. This is likely due to its identity percentage threshold defaulting to 25% when emulating a blastx style query.
I can cope with this (apologies to my disks for the I/O), but for filtering purposes, the output PSL format appears to leave me nothing but raw match counts to work with.

The PSL gives me what I want, except now I don't have information on its significance. I know that the BLAT FAQ describes a Perl script for adding "web-based identity and score calculations". I can actually get identity percentages and scores from this script, but I cannot tell what the score actually represents. Is it comparable to bit-score? How can I recover something like an e-value from a PSL output file?

It would appear that BLAT and exonerate (as suggested by @terdon) are the only pieces of software that provide enough information for "HSPFragments" to be parsed from their outputs, at least according to BioPython's documentation, anyway:

  Most search programs only have HSPs with one HSPFragment in them,
  making these two objects inseparable. However, there are programs
  (e.g. BLAT and Exonerate) which may have more than one HSPFragment
  objects in any given HSP. If you are not using these programs, you can
  safely consider HSP and HSPFragment as a single union.

It would seem the options are sticking with BLAT and generating scores afterwards with the script described in their FAQ, or using exonerate which provides a raw score as part of its default output. What these scores actually mean I will have to investigate, but for now I think the case is closed!


GC normalization ATAC seq data Is it really required to do GC normalization ATAC seq data ?
One of the paper where they have ATAC seq data did the normaisation of GC bias after peak calling. When i looked for the library and the reason why GC bias i found this
Could it be true that genes with higher GC content are higher expressed? It has been suggested that genes that are either extremely high or extremely low expressed are under some form of selection leading to “extreme” GC content. What CQN does, is making the effect of GC content comparable across samples, and we show in [1] that this leads to improved inference. It also flattens the effect of GC content on gene expression, but we believe this is better than having the effect of GC content depend on the sample.
https://www.bioconductor.org/packages/release/bioc/vignettes/cqn/inst/doc/cqn.pdf

Now that scenario is for RNA seq would it be logical to do the same for ATAC seq?

The only reason to normalize for GC content in RNA-seq is if it differs notably between samples/groups. If that's not the case and you aren't trying to compare genes withing samples then you have no reason to try to account for GC content. The same goes for ATAC-seq, though there the danger with trying to normalize for GC content is that you then mask changes in global accessibility.
In short, normalize for things that should cause problems with your analysis and nothing else.


Transcription factor targets for Human and Mouse genes I want to get transcription factor targets for my list of genes (For Human and Mouse)

I have found ENCODE Transcription Factor Targets database (The link is given below). But I am looking for some other databases also if available. I will really appreciate if anyone can suggest any other database.

http://amp.pharm.mssm.edu/Harmonizome/dataset/ENCODE+Transcription+Factor+Targets

Thank you

JASPAR is another useful transcription factor profiling database.

  The JASPAR CORE contains a curated, non-redundant set of profiles, derived from published and experimentally defined transcription factor binding sites for eukaryotes. It should be used, when seeking models for specific factors or structural classes, or if experimental evidence is paramount.

http://jaspar.genereg.net/


Is there a tool that can perform a read-group-aware mpileup from a single file? I would like to perform a samtools mpileup from a single file that contains thousands of read groups with different SM tags. I could split the bam by read group using samtools split, and then perform the mpileup, but splitting into thousands of files and then performing an mpileup from thousands of files is crazy slow on my machine.

I could imagine an implementation using pysam, but I was hoping someone has a solution already.

Bcftools mpileup already has this feature, but the feature was never implemented for the original mpileup format, which is what I need.

An example

Contents of input.sam:

@HD VN:1.6  SO:unknown
@SQ SN:chr1 LN:1000000
@RG ID:0    SM:sample_0
@RG ID:1    SM:sample_1
r0  0   chr1    24  0   1M  *   0   0   G   I   RG:Z:0
r1  0   chr1    24  0   1M  *   0   0   G   I   RG:Z:1

When I run samtools mpileup I get this behavior:

> samtools mpileup input.sam
[mpileup] 2 samples in 1 input files
chr1    24  N   2   ^!G$^!G$    II

Whereas, what I actually want is this output:

[mpileup] 2 samples in 1 input files
chr1    24  N   1   ^!G$    I   1   ^!G$    I

For a long time I thought this was a bug in samtools mpileup, but it appears that it is a feature that was described/hinted at in the samtools documentation, but was never implemented: https://github.com/samtools/samtools/issues/599#issuecomment-604941884

I broke down and finally scratched this itch. I have implemented a tool that performs read-group-aware text-pileup from a single SAM/BAM file. The tool is called streaming_pileupy and it's available from pypi for installation with pip. After installation the command would be:

spileup input.bam sampel_names.txt

The tool is pre-alpha and missing some features. Contributions are welcome.


How to extract reads with INDELs > a given size? I'm trying to modify this https://www.biostars.org/p/253774/
To get reads with deletions > 20bp
I think this gives reads with exactly 20bp dels:
samtools view -h in.bam chr1:12340-200000 \ | awk '$1 ~ "^@" || $6 ! "20D"' \ | samtools view -b - > out.bam

is it possibe to do >?
Thanks!

There's no simple way to do that with awk, since you would need to parse the CIGAR string and iterate over the operations. In python using pysam you would iterate over the reads and filter as follows:
keep = False
for OP, LEN in b.cigartuples:
    if OP == 2 and LEN > 20:
        keep = True
        break

if keep:
    ...do stuff...


Blueprint replicates On Encode, it's pretty easy to find out replicate information about the experiments, however I am having some difficulty finding this information on BLUEPRINT DDC, at http://dcc.blueprint-epigenome.eu

Where can I find information which describes technical replicates and biological (isogenic) replicates?

Thanks in advance.

Blueprint does not have technical replicates, you will find that very few sequencing based experiments ever have technical replicates. As biological replicates: each individual is regarded as a replicate.

This makes sense for the sort of experimental design which blueprint was designed to answer questions such as "Do individuals from group A have higher or lower H3K4me3 than individuals from groups B" or "Does cell type A have higher accessibility than cell type B". In each of these cases the within group variability to be considered is encapsulated in the difference between individuals.

Overall technical variation is assumed to be similar in all samples or at least not correlated with the contrast of interest (this may well not be true in the second example design I gave) it is safe to directly measure it, as long we know its distribution.

As an aside, one of the reasons that ENCODE got away with only two reps is that their reps are effectively technical (repeated preparations from the same cell lines) and thus there is next to no biological variation (which is what our algorithms are set up to find).


Is there a way to use BioLemmatizer from Python? Lemmatization is one of the most important tasks for text normalization in natural language processing (NLP), i.e. using deep learning to understand text data.
One of the good lemmatizers I found for biomedical domain is BioLemmatizer. However, it is developed in Java and is available from maven repo. I would like to use it from Python, but there seems to be no counterpart for BioLemmatizer in Python. Is there a way I can access BioLemmatizer from Python using some module like JPype. JPype, however, supports access only for the default classes in Java.

This area is huge in Python because of machine learning text processing, in fact O'Reilly have an entire book on it . However, we don't really deal with words in Bioinfo SE,

You might want to look at the nltk (natural language toolkit) for Python, but for direct parallels with the Java package bio package (your question) I don't think it's something we would know.

You could solve this via subprocessing into Java, and passing the result back to Python. I get the idea that the core ML/ DL used in natural language is set up for Python. Google for example have fantastic tensorflow stack for this.


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

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

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

};

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

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