text:
# What is the physical significance of dot & cross product of vectors? Why is division not defined for vectors?

I get the physical significance of vector addition & subtraction. But I don't understand what do dot & cross products mean?

More specifically,

• Why is it that dot product of vectors $$\vec{A}$$ and $$\vec{B}$$ is defined as $$AB\cos\theta$$?
• Why is it that cross product of vectors $$\vec{A}$$ and $$\vec{B}$$ is defined as $$AB\sin\theta$$, times a unit vector determined from the right-hand rule?

To me, both these formulae seem to be arbitrarily defined (although, I know that it definitely wouldn't be the case).

If the cross product could be defined arbitrarily, why can't we define division of vectors? What's wrong with that? Why can't vectors be divided?

• Division is the inverse of multiplication. A vector space in which you can also multiply two vectors is called an algebra (over a field). The cross product is not a type of multiplication as it is not associative. The dot product also doesn't count as multiplication as it maps two vectors into a scalar. The Quaternions are an example of a vector space which is also an algebra. – Olaf Aug 29 '11 at 11:15

I get the physical significance of vector addition & subtraction. But I don't understand what do dot & cross products mean?

Perhaps you would find the geometric interpretations of the dot and cross products more intuitive:

The dot product of A and B is the length of the projection of A onto B multiplied by the length of B (or the other way around--it's commutative).

The magnitude of the cross product is the area of the parallelogram with two sides A and B. The orientation of the cross product is orthogonal to the plane containing this parallelogram.

Why can't vectors be divided?

How would you define the inverse of a vector such that $\mathbf{v} \times \mathbf{v}^{-1} = \mathbf{1}$? What would be the "identity vector" $\mathbf{1}$?

In fact, the answer is sometimes you can. In particular, in two dimensions, you can make a correspondence between vectors and complex numbers, where the real and imaginary parts of the complex number give the (x,y) coordinates of the vector. Division is well-defined for the complex numbers.

The cross-product only exists in 3D.

Division is defined in some higher-dimensional spaces too (such as the quaternions), but only if you give up commutativity and/or associativity.

Here's an illustration of the geometric meanings of dot and cross product, from the wikipedia article for dot product and wikipedia article for cross product:

• Cross product exists in 3D and 7D only. – renormalizedQuanta Dec 26 '15 at 4:53
• If you say that direction of cross product is the orthogonal direction of the area of parallelogram then how the direction of the vector changes when we commute the vectors. – Adesh Tamrakar Sep 5 '16 at 4:50
• Why would someone give an answer that is almost completely irrelevant to the question, I wonder. This was exactly the question in my mind as I was searching about it on Google, but this definitely is not an answer to it. I wish I could down-vote it. – ThoAppelsin Mar 26 '17 at 20:41
• @renormalizedQuanta Oh yes, I keep forgetting about 7D. I'll have to watch myself and be more careful as I'm rather fond of saying that the cross product only yields a vector because of a 3D co-incidence. Which is true, but one should also mention other special dimensions. – WetSavannaAnimal May 18 '17 at 0:38
• And as far as I know, it is kinda funny (and, last time I checked, not very simple) why this only works in 3 and 7 dimensions... I'm sure there is a thread in here about it. – Vendetta Sep 19 '18 at 16:44

The best way is to ignore the garbage authors put in elementary physics books, and define it with tensors. A tensor is an object which transforms as a product of vectors under rotations. Equivalently, it can be defined by linear functions of (sets of vectors) and (linear functions of sets of vectors), all this is described on Wikipedia.

There are exactly two tensors which are invariant under rotations:

$\delta_{ij}$ and $\epsilon_{ijk}$

All other tensors which are invariant under rotations are products and tensor traces of these. These tensors define the "dot product" and "cross product", neither of which is a good notion of product:

$V \cdot U = V^i U^j \delta_{ij}$

and cross product

$(V \times U)_k = V^i U^j \epsilon_{ijk}$

It is pointless to try to think of the cross product as a "product", because it is not associative, $(A\times B)\times C$ does not equal $A\times(B\times C)$. It is also less than useful to think of the dot-product as a product in the usual sense, because it takes pairs of vectors to numbers, and $(A\cdot B)C$ does not equal $A(B\cdot C)$, because the first points in the C direction, and the second points in the A direction.

The best way is to get used to the invariant tensors. These generalize to arbitrary dimensions, and they are much clearer, and do not require a right-hand rule (this is taken care of by the index order convention). You will not find a single physics paper which uses the cross product, with the single exception of Feynman's 1981 paper "the qualitative behavior of Yang-Mills theory in 2+1 dimensions", and even if you do, it is trivial to translate.

You can divide vectors with clifford ("geometric") algebra.

The geometric product of vectors is associative:

$$abc = (ab)c = a(bc)$$

And the geometric product of a vector with itself is a scalar.

$$aa = |a|^2$$

These are all the properties required to define a unique product of vectors. All other properties can be derived. I'll sum them up, however: for two vectors, the geometric product marries the dot and cross products.

$$ab = a \cdot b + a \wedge b$$

We use wedges instead of crosses because this second term is not a vector. We call it a bivector, and it represents an oriented plane. It can be instructive to introduce a basis to see this. $e_1 e_1 = e_2 e_2 = 1$ and $e_1 e_2 = -e_2 e_1$ capture the geometric product's properties for these orthonormal basis vectors. The geometric product is then,

$$ab = (a^1 e_1 + a^2 e_2) (b^1 e_1 + b^2 e_2) = (a^1 b^1 + a^2 b^2) + (a^1 b^2 - a^2 b^1) e_1 e_2$$

As I said, the geometric product of two vectors is invertible in Euclidean space. This is obvious from the associativity property: $a b b^{-1} = a(b b^{-1}) = a$. That $b b^{-1} = 1$ implies that

$$b^{-1} = b/|b|^2$$

It's informative to look at the quantity $a = (a b) b^{-1}$, using the grouping to decompose it a different way.

$$a = (ab)b^{-1} = (a \cdot b) b^{-1} + (a \wedge b) \cdot b^{-1}$$

The first term is in the direction of $b$, the second is orthogonal to $b$. This decomposes $a$ into $a_\parallel$ and $a_\perp$.

What others have said is right, you can't define just the vector cross product to be invertible. This decomposition should convince you--you cannot fully reconstruct a vector without information from both the dot and cross products. And as has been said, this product is not commutative.

• this is a great answer, except the use of superscripts for things that are not exponentiation. – Rob Dec 30 '18 at 3:12

If you're going to define division of vectors, you're going to have to define over what multiplication field you're going to define the division: For ordinary numbers, I think of $\frac{x}{y}$ as being the number that, when multiplied by $y$, gives $x$. So, $\frac{\vec x}{\vec y}$ would have to be the vector that, when "multiplied" by $\vec y$, gives $\vec x$. If our multiplication field is the dot product, we're already in trouble, because the dot product of two vectors is a scalar, and the definition above would therefore require $\frac{\vec x}{\vec y}$ to simultaneously be a vector and a scalar.

Similarly, if our operation is the cross product, then we know that, for any vectors $\vec x$ and $\vec y$ and any scalar $c$, we have ${\vec x} \times {\vec y} = {\vec x}\times \left({\vec y} + c{\vec x}\right)$, so this means that there are an infinite number of vectors that satisfy the property "when crossproducted by $\vec y$, gives $\vec x$". Therefore, division over the cross product is not unique.

In addition to nibot's answer: division something is finding a part of something. In case of a vector, its part has the same direction but a smaller length. So it is natural to divide vectors by numbers, not by vectors.

Those dot and cross products are not simple products because they depend not only on lengths but also on orientations. They are called correspondences between a couple of vectors and numbers or vectors.

Another, more intuitive, take on your question (mathematicians: please look away) is to think of the dot (A.B = ABcosθ) and cross (AxB = ABsinθ) products as simply ways to measure degrees of parallelism of vectors versus perpendicularity (orthogonality) of vectors, in the sense that:

• The dot product of parallel unit vectors, U.U yields a number U.U cos 0 = 1 (or a number of value A*B for vectors of arbitrary vector lengths) while of orthogonal (perpendicular) vectors it’s always zero, as cos 90°=0

• Conversely, the cross product of parallel unit vectors is of magnitude 0 (as sin 0°=0) while of orthogonal unit vectors it’s 1 (or magnitude AxB for arbitrary lengths); but in this case the result is mapped into a vector, which needs to be perpendicular to the plane defined by the input vectors A and B (there’s no obvious way to assign it a direction within the plane defined by A and B). And recognise also that the cross product thereby gains a sense of handedness as AxB = -BxA, which turns out to be useful (example below).

Of course, it’s the use of sin and cos that determines the way these measures range from 0 to 1 for unit vectors; just imagine the products' values changing as you think of the A and B vectors rotating towards or away from each other for each type of product.

As far as physical significance is concerned (and ignoring the deeper insights available via Clifford Algebra etc.), these “products” turn out to be useful in so many situations that their physical significance is often taken for granted, rather than being emphasised (and maybe this underlies your question).

For the dot product: e.g. in mechanics, the scalar value of Power is the dot product of the Force and Velocity vectors (as above, if the vectors are parallel, the force is contributing fully to the power; if perpendicular to the direction of motion, the force is not contributing to the power, and it's the cos function that varies as the length of the projection of the force vector on the velocity vector varies; so it's not at all arbitrarily defined).

For the cross product: e.g. angular momentum, L = r x p (all vectors), so it seems perfectly intuitive for the vector resulting from the cross product to align with the axis of rotation involved, perpendicular to the plane defined by the radius and momentum vectors (which in this example will themselves usually be perpendicular to each other so the magnitude of rp*sin90°= rp). And if the direction of rotation changes, the sign of the momentum vector is reversed and so the cross product vector L also changes sign (hence the usefulness of mapping the cross product into a vector).

Note however, you can also calculate the number resulting from AB*sinθ (rather than mapping it into a perpendicular vector). It’s just the area of the parallelogram that is defined by the vectors A and B in the cross product.

Incidentally, there's nothing stopping you mapping the dot product into a perpendicular vector, if so desired - but it's probably not often useful to do so in physics.

As regards division, that’s a little more technical and dealt with well by the earlier replies. There’s also some approachable discussion on https://www.quora.com/Can-we-divide-a-vector-by-a-vector-and-why

I hope this is of some help for the less-expert members.

For products you have the answers. For division I recommend you to read more about quaternions. Interpretation of vectors in terms of quaternions allows for more reach algebra than vector space itself.

A little math right here. For natural definition of division you need at least division ring (one may comment that division algebra is enough, then add octonions to my answer). There is a theorem that the only finite dimensional division rings are reals, complex and quaternions. Vectors are vector space in three dimensions. So, any division for three-dimensional vectors will be "unnatural".

The question "what does $\vec{a}/\vec{b}$ equal?" is equivalent to asking "What do you multiply with $\vec b$ to get $\vec a$?" -- the answer is a matrix, assuming the multiplication involves some contraction afterwards. Equivalently, you multiply and contract a (1, 0) tensor $b^\mu$ by a (1,1) tensor $A_\mu^\nu$ to get a (0,1) tensor $b^\nu$.

But there are multiple matrices you can multiply $\vec b$ by to get $\vec a$. In two dimensions, you need two sets of "this maps to this" (and the knowledge that the mapping is linear) to pin down what the linear mapping is. In general in $n$ dimensions, you need $n$ such vectors -- so instead of dividing vectors, you divide sets of vectors -- these are called matrices.

## protected by Qmechanic♦Jul 29 '13 at 16:36

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).


text:
ChIP-seq Analysis

library(memes)
library(magrittr)
library(ggplot2)
suppressPackageStartupMessages(library(GenomicRanges))

See package website for full vignette

The Bioconductor build system does not have the MEME Suite installed, therefore these vignettes will not contain any R output. To view the full vignette, visit this article page on the memes website at this link

Introduction

In this vignette, we’ll explore using memes to deeply analyze a set of ChIP-seq peaks to identify motifs to explain differences in transcription factor binding, and consequences to chromatin accessibility at these ChIP peaks.

We will use a dataset from (Nystrom, 2020) which investigated the binding of the transcription factor E93 during Drosophila wing development. These data are useful because they contain multiple groupings with interesting motif properties between groups.

The ChIP experiments were performed under two conditions: ectopic expression of E93 during early development, and under wild-type conditions in a late stage of development. Peaks are annotated based on whether they are found in both conditions (“Endogenous” peaks), only in the ectopic expression condition (“Ectopic” peaks), or only found under wild-type conditions (“Orphan” peaks).

A key question from this observation is whether there are motifs which can distinguish these three categories of binding from eachother. To investigate this question, we will first look for enrichment of known transcription factor motifs within each binding category using the AME tool.

To start, ensure memes can detect your local install of the MEME Suite (see vignette("install_guide") for help).

check_meme_install()

Prepare peaks for analysis

A subset of annotated E93 ChIP summits are shipped with the memes package and can be loaded as follows:

data("example_chip_summits", package = "memes")
head(example_chip_summits, 5)

The ChIP summits have 3 additional annotation columns: id, peak_binding_description, and e93_sensitive_behavior. id is a unique identifier for each peak, peak_binding_description describes the condition in which E93 is observed to bind to this region (described in detail below), and e93_sensitive_behavior describes the change in chromatin accessibility observed at this site following E93 expression. These annotations will be used throughout this vignette to identify motifs that are associated with differences in E93 binding, and differences in chromatin accessibility at E93 bound sites.

The summit of a ChIP peak represents the single base-pair position with the highest amplitude ChIP signal within a ChIP peak. Typically, motif analysis gives the best results when performed in a small window (100 - 200 bp) around the summit. We can resize our summits to use a 100bp window using plyranges:

# These data use the dm3 reference genome
dm.genome <- BSgenome.Dmelanogaster.UCSC.dm3::BSgenome.Dmelanogaster.UCSC.dm3

# Get sequences in a 100bp window around the peak summit
summit_flank <- example_chip_summits %>%
# this ensures we take 50bp up & downstream of the summit for a total width of 100bp
plyranges::anchor_center() %>%
plyranges::mutate(width = 100)

Determinants of ectopic and orphan binding

One approach to identify enriched motifs is to search for enrichment of known transcription factor (TF) motifs within sequences of interest. The Fly Factor Survey is a database of known Drosophila TF motifs that we will use throught this vignette. These data are also packaged with memes:

meme_db_path <- system.file("extdata/flyFactorSurvey_cleaned.meme", package = "memes", mustWork = TRUE)

Pre-filtering database for expressed transcription factors

A key consideration when using a database of known transcription factors is whether to include motifs for TFs which are not expressed. On one hand, including all motifs allows maximum discovery of potential regulators, but on the other hand may detect factors which are not relevant under your experimental conditions. Another consideration for many tools is that increasing the number of tested motifs increases the multiple testing penaly, so limiting the motif database to a smaller number of relevant candidates can be a way to improve statistical power. There is no one right answer to this question. For this vignette, I will give an example for how one could go about filtering a database to consider only expressed factors. For other analyses, this question will need to be carefully considered.

First, we load the motif database into R using read_meme. To more easily manipulate the motif metadata, we can convert this to a data.frame using to_df.

library(universalmotif)
meme_db <- read_meme(meme_db_path) %>%
to_df()

To simplify the vignette, we’ve preproccessed RNAseq data from our tissue of interest and shipped them with memes. These are FPKM normalized counts for each transcription factor gene at early and late developmental times.

data("example_rnaseq", package = "memes")
head(example_rnaseq, 4)

Choosing a method to decide whether a gene is expressed or not is a non-trivial topic. DO NOT BLINDLY DO WHAT I DO HERE. To keep this vignette on topic, I will use an overly simplistic cutoff of FPKM > 5 to define expressed genes. In this example, I will consider any gene that exceeds the cutoff at either timepoint as expressed.

expressed_genes <- example_rnaseq %>%
# For each gene, keep only those with max(FPKM) greater or equal to 5.
dplyr::group_by(symbol) %>%
dplyr::filter(max(fpkm) >= 5)

Finally, I will filter the full motif database to select only those motifs corresponding to expressed gene.

meme_db_expressed <- meme_db %>%
# the altname slot of meme_db contains the gene symbol
# (this is database-specific)
dplyr::filter(altname %in% expressed_genes$symbol) Below is a comparison of the number of motifs in the database before and after filtering. #Number of motifs pre-filtering: nrow(meme_db) #Number of motifs post-filter: nrow(meme_db_expressed) Using this filtering strategy, we’ve decreased our multiple-testing penalty by ~54%. As a convenience function, we can set a default motif database to use for each memes function by setting the R option meme_db to a universalmotif object. This can be done as follows: # to_list() converts the database back from data.frame format to a standard universalmotif object. options(meme_db = to_list(meme_db_expressed, extrainfo = FALSE)) Alternatively, you can pass a universalmotif object to the database parameter for any memes function and this will override the global option. Examination of binding categories with AME Next, to test for enrichment of known motifs between E93 binding categories, we first collect the sequences of each category like so: by_binding <- summit_flank %>% # Get a list of chip peaks belonging to each set split(mcols(.)$peak_binding_description) %>%
# look up the DNA sequence of each peak within each group
get_sequence(dm.genome)

The data are returned as a Biostrings List, where each list entry represents the sequences of each E93 binding category.

head(by_binding)

Finally, we test each set of sequencing using AME with the runAme function. This will run AME on each set of input sequences in the Biostrings List.

ame_by_binding <- by_binding %>%
runAme

This returns a list object where each entry is the AME results for each category. For example, here is what the results for the ectopic peaks look like:

head(ame_by_binding$ectopic, 5) Visualizing AME results The plot_ame_heatmap() function provides a quick way to visualize AME results. It is built on top of ggplot2, so all ggplot2 functions can be used to further modify the plot. By default, it uses the -log10(adjusted p-value) as the heat values. See the documentation (?plot_ame_heatmap) for additional notes on customization. We can visualize the top 10 hits from the ectopic binding results as follows: ame_by_binding$ectopic %>%
dplyr::filter(rank %in% 1:10) %>%
plot_ame_heatmap(group_name = "Ectopic Sites")  +
ggtitle("Top 10 AME Hits in Ectopic Sites")

To plot results from multiple runs together, they must first be joined into a single data frame. The ame_by_binding object is a list whose names correspond to the E93 binding category. The list can be combined into a data.frame using dplyr::bind_rows. Setting .id = "binding_type" creates a new column binding_type that contains the names from the ame_by_binding list. In this way, the ame_res data.frame contains all AME results for each run, which can be distinguished by the binding_type column.

ame_res <- ame_by_binding %>%
dplyr::bind_rows(.id = "binding_type")

It is possible to aggregate results from multiple runs into a heatmap by setting the group parameter in plot_ame_heatmap(). This will stratify the y-axis by the entries in the column passed to group. This results in clustering motifs with shared hits across categories along the x-axis for quickly visualizing motifs that may be shared or distinct between multiple groups.

ame_res %>%
plot_ame_heatmap(group = binding_type)

Reducing redundant motif hits

Another key consideration for the above visualization is that in the FlyFactorSurvey database we used, different TFs can have multiple motif entries in the database which are all detected separately by AME. Here, when returning the top 5 hits from each group, you can see, for example, that a motif matching “Aef1” is reported 2 times within the top 5 hits of the ectopic sites. In this situation, it makes sense to summarize the data at the TF level, instead of the motif level. Note: There may be exceptions to this if, for example, a TF has multiple DNA binding sequences it can recognize, in which case having multiple hits may reflect a biological property of your sequences. You will have to handle this on a case-by-case basis for interesting hits and different motif databases. Here, we can see that at least for Aef1, the consensus sequences are very similar.

ame_res %>%
dplyr::filter(binding_type == "ectopic", rank %in% 1:5) %>%
head(5) %>%
dplyr::select(binding_type, rank, motif_id, motif_alt_id, consensus)

How to solve this problem will vary with different motif databases (For details on how to pre-process a motif database, see vignette("tidy_motifs")). In the tidy version of the FlyFactorSurvey database, the altnames of the motifs are set to the transcription factor gene symbol. This information is included in the AME results as the motif_alt_id column. In order to reduce the occurrance of redundant hits in the heatmap, we select the hit for each TF with the lowest p-value (i.e. the most significant hit) as follows:

ame_res %>%
# perform the next dplyr operation on each TF within each binding type
dplyr::group_by(binding_type, motif_alt_id) %>%
# within each binding type, select the TF hit with the lowest adjusted p-value
dplyr::filter(adj.pvalue == min(adj.pvalue)) %>%
plot_ame_heatmap(group = binding_type, id = motif_alt_id) +
labs(y = "Binding Category",
x = "Transcription Factor Motif")

AME Heatmap Visualization

If you have read vignette("core_ame"), the “normalized rank” heatmap was introduced as an alternative visualization method useful when AME produces very large numbers of hits, or when p-values are on very different scales between groups. The ame_compare_heatmap_methods() function can be used to visualize why.

Below is a comparison of the distribution of values when using -log10(adj.pvalue) (A) vs normalized ranks (B). Because there are relatively few hits in the results (~30), and the number of hits between groups varies more than the -log10(p-value) distributions, the “normalize” method will produce a misleading heatmap relative to the -log10(p-value) map.

ame_res %>%
dplyr::group_by(binding_type, motif_alt_id) %>%
dplyr::filter(adj.pvalue == min(adj.pvalue)) %>%
ame_compare_heatmap_methods(group = binding_type)

A third option for adjusting this visualization is to cap the heatmap scale values at a certain value. A good rule of thumb for selecting this cutoff is to view the cumulative distribution plot above in (A), and select a value on the x-axis which captures a majority of the data. Here, we see that -log10(adj.pvalue) < 10 captures >90% of the entopic & ectopic hits, but only ~50% of the orphan hits. By selecting 10 as the heatmap scale cap, we can improve the dynamic range of the signal for values <10, while preserving the information that many orphan sites generally have higher scores than the other categories. The heatmap scale can be capped at 10 by setting scale_max = 10. Below is a comparison of the 3 different visualizations. Notice how the intensity of the orphan vs other categories hits changes.

best_ame_hits <- ame_res %>%
dplyr::group_by(binding_type, motif_alt_id) %>%
dplyr::filter(adj.pvalue == min(adj.pvalue))

pval_heatmap <- best_ame_hits %>%
plot_ame_heatmap(group = binding_type, id = motif_alt_id) +
labs(x = NULL,
title = "-log10(adj.pvalue)")

norm_heatmap <- best_ame_hits %>%
plot_ame_heatmap(group = binding_type, id = motif_alt_id, value = "normalize") +
labs(x = NULL,
title = "normalize")

pval_scaled_heatmap <- best_ame_hits %>%
plot_ame_heatmap(group = binding_type, id = motif_alt_id, scale_max = 10) +
labs(x = NULL,
title = "-log10(adj.pvalue) (scale capped at 10)")
cowplot::plot_grid(pval_heatmap,
norm_heatmap,
pval_scaled_heatmap,
ncol = 1,
labels = "AUTO")

De-novo motif similarity by binding

Discovery of de-novo motifs within binding categories can be another way to identify meaningful motifs that distinguish between ectopic, entopic, and orphan sites which does not rely on known motif information. The MEME tool, Dreme, discovers short de-novo motifs in input sequences. We begin the analysis by searching for 5 de-novo motifs per binding category and shuffling the input sequences to use as the control set.

We pass the same by_binding list as the input to runDreme which will run Dreme on each set of sequences independently. The results are returned in list format, which we combine into a data.frame using dplyr::bind_rows as above.

dreme_by_binding <- by_binding %>%
runDreme("shuffle", nmotifs = 5) %>%
dplyr::bind_rows(.id = "binding_type")
# The above code chunk takes a long time to run.
# memes is packaged with the results of this run in the "example_dreme_by_binding" dataset
# which can be loaded as follows:
data("example_dreme_by_binding", package = "memes")
dreme_by_binding <- example_dreme_by_binding %>%
dplyr::bind_rows(.id = "binding_type")

Next, we want to compare the de-novo motifs between each binding category to identify motifs which could distinguish the three groups. One way to interrogate this is to examine the correlation score between each motif.

Rename the motifs to indicate the binding category they were discovered in (for visualization purposes).

dreme_by_binding_renamed <- dreme_by_binding %>%
dplyr::mutate(name = paste(binding_type, seq, sep = "_")) %>%
# update_motifs updates the information in the special motif column
update_motifs()

Create the correlation heatmap:

# Set the color values for the heatmap scale
cols <- colorRampPalette(c("white", "dodgerblue4"))(255)

# This is for adding the colored annotation blocks indicating group membership
# to the heatmap
anno.df <- dreme_by_binding_renamed %>%
dplyr::select(name, binding_type) %>%
tibble::remove_rownames() %>%
tibble::column_to_rownames("name")

dreme_by_binding_renamed %>%
# Convert to universalmotif format
to_list() %>%
# Compute the pearson correlation for each motif with all other motifs
universalmotif::compare_motifs(method = "PCC") %>%
# Plot the correlation matrix along with the annotations
pheatmap::pheatmap(color = cols,
# This sets the heatmap range to be from 0-1
breaks = seq(0, 1, by = 1/255),
annotation_col = anno.df,
# the cutree options are just cosmetic to add some spacing
# between some of the clusters
cutree_rows = 6,
cutree_cols = 6,
show_colnames = FALSE)

This heatmap reveals that most motifs discovered in each group are highly similar to the motifs found in other groups.

Test de-novo motif enrichment using AME

The above analysis suggests the motif content of the different binding categories are highly similar in sequence composition. To extend these analyses, we can use AME to test for motif enrichment of the de-novo discovered motifs within each binding category, and determine whether the motifs detected in one category are indeed enriched in another. To do this, we can provide the de-novo motifs as the AME database to test for their enrichment in each sequence category.

runAme() allows using a runDreme() results object as the database input by passing it within a list(). Naming the list() entry produces an informative motif_db name in the results data.frame.

ame_denovo_by_binding <- by_binding %>%
runAme(database = list("denovo_binding_motifs" = dreme_by_binding_renamed)) %>%
dplyr::bind_rows(.id = "binding_type")

Plotting the heatmap of results reveals that indeed a majority of the de-novo motifs discovered within a single category are detected in all 3 categories, supporting the conclusion that orphan, ectopic, and entopic sites are highly similar in sequence content.

ame_denovo_by_binding %>%
plot_ame_heatmap(group = binding_type, scale_max = 10)

However, there are 2 interesting motifs which distinguish orphan and ectopic sites from entopic sites. To help identify which TFs these motifs might belong to, we can use TomTom to match them to known TF motifs.

First, we select the motif id’s that are not found in entopic sites.

entopic_motifs <- ame_denovo_by_binding %>%
dplyr::filter(binding_type == "entopic") %>%
dplyr::pull(motif_id)

ame_denovo_binding_unique <- ame_denovo_by_binding %>%
dplyr::filter(!(motif_id %in% entopic_motifs))

Next, we use the motif id’s from the unique AME results to select those entries in the Dreme results object, and run TomTom on that subset.

dreme_by_binding_unique <- dreme_by_binding_renamed %>%
dplyr::filter(name %in% ame_denovo_binding_unique$motif_id) %>% runTomTom(dist = "ed") Finally, we visualize the TomTom results to identify candidate TFs driving the presence of ectopic and orphan sites. dreme_by_binding_unique %>% view_tomtom_hits(3) %>% cowplot::plot_grid(plotlist = ., nrow = 1, labels = "AUTO", byrow = TRUE) Motifs in opening vs closing sites Another category of ChIP peaks in this dataset are ones that overlap regions with dynamic changes to chromatin accessibility in response to E93 expression. These peaks are annotated in the e93_sensitive_behavior column which indicates whether sites have higher accessibility after ectopic E93 expression (“Increasing”), lower accessibility after ectopic E93 expression (“Decreasing”), or if the accessibility does not change (“Static”). To identify motifs associated with Increasing or Decreasing accessibility in response to ectopic E93 expression, we first grab the DNA sequences of the Increasing, Decreasing, and Static sites by splitting the ChIP peaks on the e93_sensitive_behavior column. # split by response to E93 binding by_sens <- summit_flank %>% split(mcols(.)$e93_sensitive_behavior) %>%
get_sequence(dm.genome)

Proper selection of the control set of sequences is critical to identifying biologically relevant motifs. Because we are interested in motifs associated with Increasing or Decreasing behavior, we want to exclude detection of motifs that are simply associated with E93 binding. To do this, we can use the “Static” category of E93 ChIP peaks, because these theoretically represent E93 binding sites without motifs that influence chromatin accessibility. To start, we will search for de-novo motifs within dynamic E93 binding sites using Dreme.

The by_sens list has three entries: Increasing, Decreasing, and Static. When calling runDreme on a sequence list, you can pass the name of a list entry as the control argument to use that set as the background sequences for the remaining sets of sequences. Below we search for de-novo motifs in Increasing and Decreasing peaks using the Static peaks as the background set.

# By setting "Static" as the control, runDreme will run:
# - Increasing vs Static
# - Decreasing vs Static

dreme_by_sens_vs_static <- runDreme(by_sens, "Static")

The results are returned as a list with two entries: Increasing, and Decreasing. These can be joined into a single data.frame using dplyr::bind_rows.

dreme_results <- dreme_by_sens_vs_static %>%
dplyr::bind_rows(.id = "e93_response")

Next, we match the de-novo motifs with known TF motifs using TomTom. The TomTom results will be added as new columns to the input data. TomTom data columns begin with best_match_ summarizing the top hit, and the tomtom column is a nested data.frame storing the full results for each input motif.

denovo_dynamic_motifs <- dreme_results %>%
runTomTom()

We can visualize some of the top TomTom hits for each de-novo motif using view_tomtom_hits(). Passing a number to view_tomtom_hits() shows that many matches ranked in descending order. Passing no value will show all hits. It’s a good idea to inspect all hits visually. This can be an extremely revealing step, as it will give you insight into which motifs are generally similar to your input motif. You may also discover proteins with identical motifs, in which case you may need to decide which TF is the better candidate. For brevity, I’ll only show the top 3 hits for each motif.

denovo_dynamic_motifs %>%
# Plot the top 3 matches
view_tomtom_hits(3) %>%
# Arrange into figure panels
# You don't need to do this if you're just exploring your data, I do it here to make a pretty figure
cowplot::plot_grid(plotlist = ., labels = "AUTO")

After manually inspecting the above hits, we may decide that a lower-ranked motif is a better match by eye, or because of domain-specific knowledge that makes a different TF a better candidate. We can swap the best matched motif using force_best_match(). This will update the best_match_ columns and re-rank the tomtom nested data column to set a lower ranked match to the top hit. The updates best_match_ columns will reflect the values corresponding to the new match.

As an example, I’ll set the top hit for m03_AKGG to the lower ranked pho_SANGER_10 since it has weaker consensus sequence in the regions that don’t overlap the de-novo motif.

# use force_best_match to use update the best match info
# %<>% is a shortcut from magrittr for updating the data in-place
denovo_dynamic_motifs %<>%
force_best_match(c("m03_AKGG" = "pho_SANGER_10"))

Replotting the top hits, you can see that “pho_SANGER_10” is now listed as the top hit for “m03_AKGG”.

denovo_dynamic_motifs %>%
view_tomtom_hits(1) %>%
cowplot::plot_grid(plotlist = ., labels = "AUTO")

Finally, we can make publication quality figures using cowplot. This code is a little complicated, but hopefully gives an idea for how you can use R graphics to generate nice motif figures!

A striking result from this is that the de-novo motifs associated with Decreasing sites match the E93 motif itself.

denovo_dynamic_motifs %>%
# Create a new label that indicates whether the motif is found in Increasing or Decreasing sites
dplyr::mutate(label = paste0(e93_response, " in response to E93")) %>%
# Sort the top hits first
dplyr::arrange(rank) %>%
# Run the plot command separately for Increasing and Decreasing sets
split(.$label) %>% # For each set: purrr::imap(~{ # Plot the top hit for each motif in the set as a panel figure top_hits <- view_tomtom_hits(.x, 1) %>% cowplot::plot_grid(plotlist = ., nrow = 1, labels = "AUTO") # Create a figure title from the label we created above title <- cowplot::ggdraw() + cowplot::draw_text(.y) # Combine the title & the motif hits figure cowplot::plot_grid(plotlist = list(title, top_hits), ncol = 1, rel_heights = c(0.1, 1) ) }) Scanning for motif matches using FIMO Because we discover a de-novo motif from DREME that matches to E93 in Decreasing sites, this suggests that the E93 motifs in Decreasing sites have different properties relative to E93 motifs in Increasing or Static sites. For example, the E93 motif may have higher similarity to the canonical sequence, may be present in greater number, or may have different positioning in Decreasing sites relative to Increasing or Static sites. To test these questions directly, we can identify individual E93 motif matches within each sequence and investigate their properties. The FIMO tool is used to identify the positions of motif matches in a set of sequences. It will return the coordinates, sequence, and a score (higher values are more similar to the consensus) for each motif found in the input sequences. In order to scan for the E93 motif within target sequences, we first need a copy of the E93 motif. Although we could import the motifs from our local meme database, it is also possible to use motifs pulled from a MotifDb query as follows. I’ll use that approach here as an example, but if this were a real analysis, I’d take the result from our local database. e93_flyfactor <- MotifDb::MotifDb %>% # Query the database for the E93 motif using it's gene name MotifDb::query("Eip93F") %>% # Convert from motifdb format to universalmotif format universalmotif::convert_motifs() %>% # The result is a list, to simplify the object, return it as a single universalmotif .[[1]] # Rename the motif from it's flybase gene number to a more user-friendly name e93_flyfactor["name"] <- "E93_FlyFactor" # Display the cleaned motif e93_flyfactor Just like runDreme() and runAme(), runFimo() takes a Biostrings::XStringSet as input. Using get_sequence() to generate these for DNA sequences is preferred because it creates sequence inputs that can be parsed into genomic coordinates by Fimo. fimo_res <- summit_flank %>% get_sequence(dm.genome) %>% runFimo(motifs = e93_flyfactor, thresh = 1e-3) The results of runFimo() are returned as a GRanges object containing the positions of each motif discovered in the input sequences. The best way to integrate these data with our input peaks is to use the plyranges suite of tools for performing overlaps and joins between GRanges objects. fimo_res Counting the number of motifs per peak plyranges extends dplyr-like syntax to range objects. Here we add a count for each motif per peak in the n_motifs column. We also add a column has_motif which will be a binary indicator of whether a peak contains any motifs. summit_flank %<>% # Note, if running FIMO with multiple motifs, this solution will not work # as it will count all motifs within the fimo-results without splitting by motif_id plyranges::mutate(n_motifs = plyranges::count_overlaps(., fimo_res), has_motif = n_motifs > 0) First, we want to determine whether E93 sensitive sites are more likely to have E93 motifs in certain response types. Here we can see that sensitive decreasing sites are more likely to have E93 motifs than sensitive increasing or insensitive static sites. summit_flank %>% data.frame %>% dplyr::mutate(has_match = dplyr::if_else(has_motif, "Match", "No Match")) %>% ggplot(aes(e93_sensitive_behavior)) + geom_bar(aes(fill = forcats::fct_rev(has_match)), position = "fill") + scale_fill_manual(values = c("Match" = "firebrick", "No Match" = "Black")) + labs(fill = "E93 Motif Match", y = "Fraction of Sites", x = "Response to E93 binding") To investigate whether E93-sensitive sites that increase, decrease, or remain static have different numbers of E93 motifs, we plot the fraction of sites with each number of motifs per group. Here we can see that in addition to being more likely to contain an E93 motif, sensitive decreasing sites are more likely to to contain 2 or more matches, where 10% contain at least 2 motifs. summit_flank %>% # currently, group operations are faster as data.frames, so we convert to data.frame data.frame %>% dplyr::group_by(e93_sensitive_behavior, n_motifs) %>% dplyr::count() %>% dplyr::group_by(e93_sensitive_behavior) %>% dplyr::mutate(frac = n/sum(n)) %>% ggplot(aes(n_motifs, frac)) + geom_line(aes(color = e93_sensitive_behavior), size = 1) + labs(y = "Fraction of Sites", x = "Number of E93 Motifs", color = "Response to E93 Binding") + theme_bw() Finally, we want to assess whether the quality of E93 motifs is different between sensitivity categories. To examine this, we need to determine which motifs are found in which peaks. We use plyranges::join_overlap_intersect to return motif entries appended with peak-level metadata, like the peak id each motif is found within. # return position of each motif match w/ peak metadata intersect <- fimo_res %>% plyranges::join_overlap_intersect(summit_flank) We use the FIMO score as a proxy for quality, where higher scores are better matches to the motif. Here we examine only the best match (highest score) motif within each peak. best_motifs <- intersect %>% # group by ChIP peak id plyranges::group_by(id) %>% # Keep only motifs w/ the highest score within a peak plyranges::filter(score == max(score)) %>% plyranges::ungroup() Finally, visualize the FIMO scores for the best E93 motif within each ChIP peak across categories. This reveals that Decreasing sites appear to have higher FIMO scores, suggesting the E93 motifs within Decreasing sites are higher quality E93 motifs than those in Increasing or Static sites. best_motifs %>% data.frame %>% ggplot(aes(e93_sensitive_behavior, score)) + geom_boxplot(aes(fill = e93_sensitive_behavior), notch = TRUE, width = 0.5) + guides(fill = "none") + labs(x = "Response to E93", y = "FIMO Score") + theme_bw() We can test whether the distribution of FIMO scores between multiple groups is statistically significantly different using an Anderson-Darling test. For a pairwise comparison (in case you only have 2 groups to test), a Kolmogorov-Smirnov test is commonly used instead. The PMCMRplus::adAllPairsTest() performs the Anderson-Darling test and corrects for multiple testing. This reveals that Decreasing FIMO scores are significantly different from Increasing and Static scores (p < 0.05), but Increasing and Static scores are not significantly different from eachother (p > 0.05), supporting the conclusion that E93 motifs in Decreasing sites are higher quality than those in Increasing or Static sites. best_motifs %>% data.frame %>% dplyr::mutate(behavior = factor(e93_sensitive_behavior)) %>% PMCMRplus::adAllPairsTest(score ~ behavior, ., p.adjust.method = "fdr") We can visualize the differences in E93 motif quality between groups by reconstructing PCM’s from the matched motif sequences. Here, it is clear that although each category contains an E93 motif match, the sequence content of those motifs differs between categories, providing a nice visual representation of the statistical analysis above. best_motifs %>% # Get sequences of each matched motif add_sequence(dm.genome) %>% # Split by E93 response category split(mcols(.)$e93_sensitive_behavior) %>%
# Convert from GRangesList to list() to use purrr::imap
as.list() %>%
# imap passes the list entry as .x and the name of that object to .y
# here, .y holds the e93_sensitive_behavior name, and .x holds the motif positions
purrr::imap(~{
# This builds a new motif using the matched sequences from each group
create_motif(.x$sequence, # Append the response category to the motif name name = paste0("E93_", .y)) }) %>% # Prepend E93 Fly Factor motif to beginning of motif list so it is plotted at the top c(list("E93_FlyFactor" = e93_flyfactor), .) %>% # Plot each motif view_motifs() We can also build a more complex representation of these data by making a heatmap of the sequence variability for each matched sequence using plot_sequence_heatmap(). This function will plot each sequence as a row of a heatmap where each column represents the position relative to the motif start coordinate. The cells are colored by the DNA sequence found at that position for each sequence. Using this type of plot it becomes even more clear how variability at position 6 differs between Decreasing and Increasing sites, where Decreasing sites have a very consistent “C” at position 6, while Increasing sites are more likely to have non-C bases at this position. best_motifs %>% # Get sequences of each matched motif add_sequence(dm.genome) %>% data.frame %>% # Split by E93 response category split(.$e93_sensitive_behavior) %>%
# Extract the "sequence" column from each category
purrr::map("sequence") %>%
# Plot a heatmap for each category
plot_sequence_heatmap(title_hjust = 0.5) %>%
# Arrange into a grid
cowplot::plot_grid(plotlist = ., nrow = 1)

Centrality of E93 motif

Next we want to visualize whether the E93 motif has altered positioning across the response categories. To do this we add the metadata for each nearest motif to our peak summits. Setting distance = TRUE in plyranges::join_nearest adds a distance column indicating the distance to the nearest joined region.

summit_nearest_e93_motif <- summit_flank %>%
plyranges::anchor_center() %>%
plyranges::mutate(width = 1) %>%
plyranges::join_nearest(fimo_res, distance = TRUE)

We can visualize these data as a cumuluative distribution plot. In this example, it does appear that decreasing sites tend to be closer to the peak summit than increasing or static sites.

summit_nearest_e93_motif %>%
data.frame %>%
# Limit our analysis to peaks that contain motifs
dplyr::filter(has_motif == TRUE) %>%
ggplot(aes(distance)) +
stat_ecdf(aes(color = e93_sensitive_behavior), size = 1, pad = FALSE) +
labs(x = "Distance to nearest E93 Motif",
y = "Fraction of Sites",
color = "E93 Response") +
theme_linedraw()

Conclusion

This is just one example of how memes can integrate multiple data sources to perform deep analysis of genomics data. Here we went from a set of annotated ChIP peaks do multifaceted interrogation of motif content. Each approach has advantages and disadvantages. The most important thing is to critically evaluate the data at every step: visualize de-novo motif hits with view_motifs, compare matches using view_tomtom_hits, and think carefully about which sequences to use as control regions.

References

Nystrom SL, Niederhuber MJ, McKay DJ. Expression of E93 provides an instructive cue to control dynamic enhancer activity and chromatin accessibility during development. Development 2020 Mar 16;147(6).

Session Info

sessionInfo()

format: plain


text:
## Expectations (16/365)

Sometimes expectations and random variables appear in so many ways that I find it a little confusing at times. The book I am following has a neat diagram that helps quite a bit.

## Basic Expectation

Start with a distribution $P(\cdot)$ for sample space $\Omega$ and a random variable $\theta$. The basic expectation tells us the value $\theta$ is likely to take on average.

$\displaystyle E\theta = \sum_\omega \theta(\omega) P(\omega) \\ \text{Notice that } E\theta = E[\theta I_\Omega]$

Note that $\theta$ induces a decomposition $D_1, \dots, D_n$ of $\Omega$ as follows

$\displaystyle E\theta = \sum_{i=1}^n x_i P(D_i) \text{ where } D_i = \{ \omega | \theta(\omega) = x_i \}$

## Conditional Expectation

Instead of $P(\cdot)$ we could be given a distribution $P(\cdot | D)$ with respect to an event $D$. The expectation over this is the value $\theta$ is likely to take on average conditioned on the event $D$ (i.e. restricted to).

$\displaystyle E(\theta | D) = \sum_\omega \theta(\omega) P(\omega | D) \\ \text{Notice that } E(\theta | D) = \frac{E[\theta I_D]}{P(D)}$

Let’s suppose that we have an event $A$ and its conditional probabilities $P(A | D_1), \dots, P(A | D_n)$ where $\{ D_i \}$ is a decomposition of $\Omega$. We can write this as a random variable that takes on the value $P(A | D_i)$ whenever $\omega \in D_i$.

$\displaystyle P(A | \mathcal{D})(\omega) = \sum_{i=1}^n P(A | D_i) I_{D_i}(\omega)$

Now that we have this random variable, we can once again take its expectation which in this case gives the probability of $A$ on average conditioned on an event from $\mathcal{D}$. This has a special name and is called the total probability.

$\displaystyle E P(A | \mathcal{D}) = \sum_{i=1}^n P(A | D_i) P(D_i) = P(A)$

## One More Generalization

Suppose now that we have a random variable $\theta$ inducing the decomposition $A_1, \dots, A_m$ where $A_j = \{ \omega | \theta(\omega) = x_j \}$. We also have conditional probabilities $P(A_j | D_i)$. We can certainly take the following expectation from before

$\displaystyle E(\theta | D_i) = \sum_j x_j P(A_j | D_i)$

which is the expectation of $\theta$ conditioned on $D_i$. We now do this for the entire decomposition $\mathcal{D}$ to arrive at this random variable

$\displaystyle E(\theta | \mathcal{D})(\omega) = \sum_i E(\theta | D_i) I_{D_i}(\omega)$

which takes on a conditional expectation of $\theta$ at each $\omega$. We can now generalize the total probability formula such that it gives the the expectation of $\theta$ on average conditioned on an event from $\mathcal{D}$.

$\displaystyle E(A | \mathcal{D}) = P(A) \\ \text{generalizes to} \\ EE(\theta | \mathcal{D}) = E\theta$

The lesson here is to always translate the things we do with probabilities to expectations.

This entry was posted in Uncategorized and tagged , . Bookmark the permalink.


text:
Browse Questions

# Which one of the following is non-ideal solution?

$CCl_4 + CHCl_3$
Hence (D) is the correct answer.


text:
# double periodic functions with real variables

$$u(x+X,y)=e^{iky}u(x,y)\ u(x,y+Y)=u(x,y)$$

This is a quasi double periodic boundary condition. x and y are real numbers. I’m wondering whether there exists a general formula of real variable functions under such boundary condition, just like elliptic functions. I’d appreciate it very much.


text:
# Proving the convergence of the maximum of Uniform Distribution

I have a random sample of size $$X_1, X_2, .., X_n$$ following $$U(0,2)$$. I need to prove that $$X_{(n)}$$ which is the maximum ordered statistics will converge to $$2$$ in probability and almost surely.

I know that I can just prove this to converge to 2 almost surely and from there I can make a statement that since it converges almost surely, it will converge in probability and in distribution too.

But the definitions of convergence in probability and almost sure convergence looks identical to me. I could prove that this this maximum ordered statistic will converge to $$2$$ in probability using definition as follows:

$$\lim_{n \to \infty} P(|X_n-2|<\epsilon) = \lim_{n \to\infty} P(-\epsilon+2

The CDF of the largest observation is $$(\frac{x}{2})^n$$. Hence, the above probability can be written as:

$$\lim_{n\to\infty} P(-\epsilon+2

Hence, I could see that the largest observation will converge in probability.

How do I go about proving the almost sure convergence. I know the definition. I don't want to write it down again but I cannot separate between this definition and the one I have read for almost sure. It just puts the limit inside the probability statement unlike this definition which puts it before probability.

$$\sum_{n=1}^{\infty}P[|X_{(n)}-2|>\epsilon]=\sum_{n=1}^{\infty}P[X_{(n)}<2-\epsilon] =\sum_{n=1}^{\infty}\left(\frac{2-\epsilon}2\right)^n=\frac{2-\epsilon}{\epsilon}$$
• Hello John. Thanks for the answer. I understand the defition there and also understand your steps. Now, we need to show that the result is less than $\infty$. How can we do that? Since, $\epsilon$ is a small number and hence $\frac{2}{\epsilon}$ will be large. Something large minus 1 is still large and how do we know that it is less than $\infty$ – userNoOne Feb 20 at 14:36


text:
# Fusing Operations¶

Another approach we can use to help improve tensor computations is to specialize commonly used combinations of operators. Combining operations can eliminate unnecessary intermediate tensors. This is particularly helpful for saving memory.

There is a lot of ongoing work on how to do this operator fusion automatically. However, for very common operators. it is worth just writing these operators directly and speeding them up.

In minitorch we will do this fusion by customizing special operations. If we find it useful, we will add these operations in Numba.

## Example: Matrix Multiplication¶

Let's consider a matrix multiplication example. In past modules we have done matrix multiplication by applying a broadcasted zip and then a reduce. For example, consider a tensor of size (2, 3, 1) and a tensor of size (3, 2). We first zip these together with broadcasting to produce a tensor of size (2, 3, 2):

And then reduce it to a tensor of size (2, 1, 2), which we can view as (2, 2):

This computation is correct, but it has two issues. First, it needs to create a tensor of size (2, 3, 2) in the intermediate step. These tensor can be an order of magnitude larger than either of the tensors that were passed in. Another more subtle problem is that the reduce operator may need to call save_for_backwards on this tensor, which means it stays in memory throughout all of our forward computations.

An alternative option is to fuse together these two operations. We will create a single matmul operations using python @ operator. We can then skip computing the intermediate value and directly compute the output. We can do this by writing a specialized tensor Function with a forward function that directly produces the output and a backward that produces the required gradient.

This allows us to specialize our implementation of matrix multiplication. We do this by 1) walking over the output indices, 2) seeing which indices are reduced, and 3) then seeing which were part of the zip. Here's what that looks like in our diagram.

Given how important this operator is, it is worth spending the time to think about how to make each of these steps fast. Once we know the output index we can very efficiently walk through the input indices and accumulate their sums.

A fast matrix multiplication operator can be used to compute the forward. To compute backward we need to reverse the operations that we computed. In particular this requires zipping grad out with the alternative input matrix. We can see the backward step by tracing the forward arrows:

It would be a bit annoying to optimize this code as well, but again luckily, to compute the backward step , we can reuse our forward optimization. In fact we can simply use the following identity from matrix calculus that tells use backward can be computed as a transpose and another matrix multiply.

$\begin{split}\begin{eqnarray*} f(M, N) &=& M N \\ g'_M(f(M, N)) &=& grad_{out} N^T \\ g'_N(f(M, N)) &=& M^T grad_{out} \end{eqnarray*}\end{split}$


text:
# Can I modify Ctrl— to insert Indexed expressions?

For me, Indexed[x, n] is a much more useful and versatile construct than Subscript[x, n]. I would therefore like to modify Ctrl-- to insert Indexed expressions instead of Subscripts. Is this possible?

Note that I do, occasionally, need to work with code containing Subscripts, so simply redefining Subscript[x__] := Indexed[x] is not an option for me.

• May I ask, why do you prefer one to another? Where do you used such indexed expressions: in Input cells? In DisplayFormula cells? Somewhere else? What is the advantage? – Alexei Boulbitch Feb 18 '15 at 9:14
• @AlexeiBoulbitch I use them in input cells, primarily as a substitute for Part's ugly double-bracket notation. They're also useful for integration over multidimensional regions, e.g., Integrate[Indexed[x, 1]^2 Indexed[x, 2]^2, Element[x, Rectangle[]]]. – David Zhang Feb 19 '15 at 5:59
• @Mr.Wizard This my question is to the whole community, by I just directed it to Mr.Wizard to address one. Previously it was known that Mma "does not like" to operate with expressions with subscripts. It was a general knowledge that it may occasionally evaluate an expression containing subscribed variables correctly, but may also return a mistake. So we preferred to use, say, x1 , rather than Subscript[x,1] . The question is, if this is not a problem for calculations with indexed variables, as OP states? – Alexei Boulbitch Feb 19 '15 at 8:50
• @AlexeiBoulbitch Could I ask you for an example of Mathematica giving incorrect results in computations with Subscripts? I've never run into any such problems with either Subscript or Indexed, but my use cases of both are admittedly limited. – David Zhang Feb 19 '15 at 9:06
• I never tried anything with Indexed before, in fact I did not know about its existence until now. So this I cannot. With Indexed, however, I and other people tried many times. When I started with Mma3 it was simply a "no go". It seemed that in the last versions (9 and 10) it is not the case, but I recently tried to make calculations in the 10.0.2 with the subscripts and it, indeed, made several operations correctly, but then returned an erroneous result which immediately turns into the right one, as soon as I removed the subscripts. I did not save this example though. – Alexei Boulbitch Feb 19 '15 at 9:16

CreatePalette @ PasteButton @ Indexed[\[SelectionPlaceholder], \[Placeholder]]

AppendTo[
CurrentValue[$FrontEndSession, "InputAliases"], "idx" -> TemplateBox[{"\[SelectionPlaceholder]", "\[Placeholder]"}, "IndexedDefault"] ]  For a system-wide keyboard shortcut you can copy KeyEventTranslations.tr from SystemFiles\FrontEnd\TextResources\(your OS) within $InstallationDirectory to the like path within $UserBaseDirectory and replace the line containing "KeypadSubtract" with: Item[KeyEvent["KeypadSubtract", Modifiers -> {Control}], FrontEndExecute[{ FrontEndNotebookWrite[FrontEndInputNotebook[], TemplateBox[{"\[SelectionPlaceholder]", "\[Placeholder]"}, "IndexedDefault"], After] }]],  This will allow insertion of the same template as above. I am still seeking a way to emulate the default Subscript behavior where a template is inserted for a blank line, but the first element is auto-filled if an expression exists. With the hint from Kuba to use NotebookApply we can get nearly perfect emulation with: Item[KeyEvent["KeypadSubtract", Modifiers -> {Control}], FrontEndExecute[{ FrontEndSelectionMove[FrontEndInputNotebook[], All, Word]; FrontEndNotebookApply[FrontEndInputNotebook[], TemplateBox[{"\[SelectionPlaceholder]", "\[Placeholder]"}, "IndexedDefault"]] }] ],  Related: • @David The first snippet creates a Palette Notebook -- look for a small window with a button; on my system it appears in the upper right hand corner. Yes, it should be possible to map this to a key combination, however (1) I need to modify my configuration to test it, and (2) I was trying to think of a way to emulate the (default) behavior with Subscript, rather than inserting the double placeholder object. I shall update this answer later with whatever I can get working. – Mr.Wizard Feb 18 '15 at 0:09 • @David Please see update – Mr.Wizard Feb 18 '15 at 0:39 • @Mr.Wizard You can use this: SelectionMove[nb, Previous, Character]; SelectionMove[nb, All, Character]; NotebookApply[nb, TemplateBox[{"\[SelectionPlaceholder]", "\[Placeholder]"}, "IndexedDefault"]] – Kuba Feb 18 '15 at 9:02 • @Kuba Thanks! NotebookApply is what I needed. However I get better behavior from All, Word as shown in the present answer; can you confirm? – Mr.Wizard Feb 19 '15 at 2:14 • @Kuba This is fantastic! I'm surprised you found a way to emulate the selection behavior of Subscript. I have, however, noticed something strange: while I can override KeypadSubtract and _ just fine, I can't seem to modify the behavior of -. Is this intended? – David Zhang Feb 19 '15 at 5:32 One way to do what it seems me you would like to is to use the InputAliases function by evaluating in your notebook the following expression:  SetOptions[InputNotebook[], InputAliases -> {"ind" -> "Evaluate[Indexed[\[SelectionPlaceholder],\[Placeholder]]]"}]  After that as soon as you type Esc+ind+Esc placeholders to type in an indexed expression will appear in the place of the cursor. This is, however, boring, since you will need to evaluate such an operation in each notebook before its use. One may also do it once for ever as follows: 1. Having any of your notebooks open go to Menu/Format/OptionInspector (alternatively press Shift+Ctrl+O). This opens the OptionInspector. 2. In its slot entitled "Show option values" choose Global Options. 3. In the slot "Lookup" type InputAlaises, or alternatively in the left window choose Edit Options 4. In the right window choose InputAlaises and click on the gray box at the right margin of the line. 5. This opens a dialog box with a list of the already existing aliases. Press the button "Add". 6. This opens the "Edit setting rules" dialog box enabling one to introduce a new rule. By default x->y stays there. 7. Remove x and insert the key combination that you need, say, in. (do not put simply i, since such an aliase already exists, I checked). 8. Instead of y insert TemplateBox[{"\[SelectionPlaceholder]", "\[Placeholder]"}, "IndexedDefault"]. Click the "OK" button on the both dialog windows, and "Apply" on top of the OptionInspector window. Done. Now as soon as you type anywhere in any of your notebooks Esc+in+Esc you will get the placeholders enabling you to type in the expression which will appear to be indexed. Have fun! • One can add a persistent input alias far more easily; simply evaluate AppendTo[ CurrentValue[$FrontEnd, "InputAliases"], "idx" -> TemplateBox[{"\[SelectionPlaceholder]", "\[Placeholder]"}, "IndexedDefault"] ]. For more information see: (72989) – Mr.Wizard Feb 19 '15 at 16:19


text:
# Concerning aerobic respiration by animals in the light vs dark – be able to explain if an expected difference in

###### Question:

Concerning aerobic respiration by animals in the light vs dark – be able to explain if an expected difference in rate of respiration (CO2 concentration) exists, and why or why not?

#### Similar Solved Questions

##### Use the regression line to predict for birth weight; in grams; of a baby born to an 18- year-old mother: Enter your answer to the nearest whole number:3,250
Use the regression line to predict for birth weight; in grams; of a baby born to an 18- year-old mother: Enter your answer to the nearest whole number: 3,250...
##### What is the wavelength of the radio transmission?
FM station KCOM transmits simultaneously, at a frequency of 99.9 MHz, from two broadcast towers placed precisely 31.5 m apart along a north-south line.a. What is the wavelength of KCOM's transmissions?b. Suppose you stand 90.0 m due east of the point halfway between the two towers with your port...
##### Explain how the Fed injecting money into the economy can bring on the liquidity effect How...
Explain how the Fed injecting money into the economy can bring on the liquidity effect How does it increase and or decrease interest rates?...
##### 9+.2 Pincl_+he Shepe c? Exercise Lne 40_422 greph The] e Hangen t Fua ckon 4+ Ine 2 f Le Gien pamn +36- â‚¬90)
9+.2 Pincl_+he Shepe c? Exercise Lne 40_422 greph The] e Hangen t Fua ckon 4+ Ine 2 f Le Gien pamn + 36- â‚¬ 90)...
##### Question: Let U span{X, = [1,0,0j7 , Xz = (-2,1,0j7 , Xs = [1,-1,4F7}Use the Gram-Schmidt algonthm t0 find an orthogonal basis Of U{EJ EL [} {E] [] [} {EJGH [H} {LH] [UJ EH {EL[H] EJ}
Question: Let U span{X, = [1,0,0j7 , Xz = (-2,1,0j7 , Xs = [1,-1,4F7} Use the Gram-Schmidt algonthm t0 find an orthogonal basis Of U {EJ EL [} {E] [] [} {EJGH [H} {LH] [UJ EH {EL[H] EJ}...
##### Fill in the blanks.$(x-2)^{2}+(y+5)^{2}=9:$ center ( __________ ), ( __________ ); radius __________
Fill in the blanks. $(x-2)^{2}+(y+5)^{2}=9:$ center ( __________ ), ( __________ ); radius __________...
##### Chapter 12 Assignment Assignment: Chapter 12 Assignment Assignment Score: 0.00% Save Submit Assignment for Grading Questions...
Chapter 12 Assignment Assignment: Chapter 12 Assignment Assignment Score: 0.00% Save Submit Assignment for Grading Questions Exercise 12.07 (Self-Test) Algorithmic <Question 2 of 4 Check My Work (7 remaining) eBook Video The following table contains observed frequencies for a sample of 200. Colum...
##### 8) While the hydroxides and carbonates of sodium and potassium are readily soluble in water; the corresponding salts of magnesium and calcium are only sparingly soluble. Explain why this is So (2) 9) The alkali and few alkaline earth metals impart colour to a flame a property used in the flame test experiment; however; Be and Mg do not impart colour to a flame. Give reason (2) 10) Write a short note on biological relevance of sodium and calcium H) In aqueous solution, Beryllium ion becomes attac
8) While the hydroxides and carbonates of sodium and potassium are readily soluble in water; the corresponding salts of magnesium and calcium are only sparingly soluble. Explain why this is So (2) 9) The alkali and few alkaline earth metals impart colour to a flame a property used in the flame test ...
##### The following transformations are applied to the function yWrite the expression for thisfunction_No T-axis reflection_ Shifted units down; Shifted units to the right:f(c)This function exhibit: [Solecl an answorsince the base Soloctan answor
The following transformations are applied to the function y Write the expression for this function_ No T-axis reflection_ Shifted units down; Shifted units to the right: f(c) This function exhibit: [Solecl an answor since the base Soloctan answor...
##### Hi-Way Ltd. Apportionment of Shared Cost You are working as a Management Accountant in Hi-Way Ltd....
Hi-Way Ltd. Apportionment of Shared Cost You are working as a Management Accountant in Hi-Way Ltd. a manufacturing company, it has two production departments. Machining and Assembly, and two service departments, Repairs and Quality Control. The following information is available about production ove...
##### Is this correct QUESTION 40 The following factor can prevent the attachment of C5b67 and eventual...
is this correct QUESTION 40 The following factor can prevent the attachment of C5b67 and eventual formation of the MAC attack: a. C3b b. S protein O c. cia O d. properdin...
##### Are there any equation for 2d motion and if present what are they?
Are there any equation for 2d motion and if present what are they?...
##### 1. Dave currently has 82.5% on the homework, 9/18 on Test 1. 28/40 on Text 2....
1. Dave currently has 82.5% on the homework, 9/18 on Test 1. 28/40 on Text 2. 36/40 on Test 3, and did 4/5 Quizzes. Assume the total course grade is calculated just as in your course. The tests are weighted the same and the quizzes are assigned only completion grades. The marks for the course are as...
Check my work A portfolio manager invested $2,300,000 in bonds. In one year, the market value of the bonds dropped to$2,280,000. The interest payments during the year totaled $92,000. a. What was the manager's total rate of return for the year? (Round your intermediate calculations to 4 decimal... 5 answers ##### Graphlng the results from klnetlcs experlments wlth enzyme Inhlbltors The followlng klnetlc data were obtalned for an enzyme In the absence of Inhlbltor (1) and In the presence of two dlfferent Inhlbltors (2) and (3) at 5 mM concentratlon: Assume [ET] Is the same In each experlment:VS]k(mM{v(umolmL-sedv (umolmL-sedvumol/mL-sed1124.35.5202913516012 4018Plot these data a5 double-reclprocal Llneweaver-Burk plots and use your graph to answers a. and b Determlne Vmax and Km for the erzyme b.Deduce th Graphlng the results from klnetlcs experlments wlth enzyme Inhlbltors The followlng klnetlc data were obtalned for an enzyme In the absence of Inhlbltor (1) and In the presence of two dlfferent Inhlbltors (2) and (3) at 5 mM concentratlon: Assume [ET] Is the same In each experlment: VS] k(mM{v(umolm... 1 answer ##### 3. R, S Configurations Assign R, S configurations to each indicated chirality center in the molecules... 3. R, S Configurations Assign R, S configurations to each indicated chirality center in the molecules below. 14. OH 15. COOH ---H CH NH2 HỌN H -CH3 HO alanine norepinephrine - 18. 16. HOC / OHH HV CH; НС . H7 HO COH CH2 dihydrocarvone tartaric acid... 1 answer ##### Describe the barriers, cells, and proteins involved in the different phases of the immune response Describe the barriers, cells, and proteins involved in the different phases of the immune response... 5 answers ##### In Exercises$19-36,$determine whether the equation represents$y$as a function of$x .$x^{2}+y=4$
In Exercises $19-36,$ determine whether the equation represents $y$ as a function of $x .$ $x^{2}+y=4$...
##### 3.111 3.6 m IN 3.111 Two cables exert forces of 90 kN each on a truss...
3.111 3.6 m IN 3.111 Two cables exert forces of 90 kN each on a truss of weight W = 200 KN. Find the resultant of these three forces acting on the truss and the point of intersection of its line of action with line AB. 30 3.6 m 90 KN 3.6 m Fig. P3.111...
##### Lysergic acid diethylamide (LSD), a potent hallucinogen and potential treatment for depression and Alzheimer' $disease contains three - nitrogen types: nitrogen (B), the nitrogen embedded within the 6-membered ring (C), and the indole nitrogen (4) Analogous to pyrrole, the indole nitrogen (A) lone pair of electrons are part of an aromatic system. From left to right; rank the nitrogen'$ B, and C) in order of increasing basicity and provide rationale invoking hybridization; resonance_ a
Lysergic acid diethylamide (LSD), a potent hallucinogen and potential treatment for depression and Alzheimer' \$ disease contains three - nitrogen types: nitrogen (B), the nitrogen embedded within the 6-membered ring (C), and the indole nitrogen (4) Analogous to pyrrole, the indole nitrogen (A) ...
##### F1 Points]DETAILSSCALCET? 5. 5 039Ealuate the indefinite integral: (Use â‚¬ for the constant ofintegraton EnQ 0 12 CosRELEEKAHMHRHLM
F1 Points] DETAILS SCALCET? 5. 5 039 Ealuate the indefinite integral: (Use â‚¬ for the constant ofintegraton EnQ 0 12 Cos RELEEKA HMHRHLM...
##### Find the exact value of the trigonometric cxpression without the usc of calculatorSuppose that the termlnal side of angle lies Quadrant" and the termninaol angleEQuadrant III M cCs 0 =and sin 0find the exact value of cos (a
Find the exact value of the trigonometric cxpression without the usc of calculator Suppose that the termlnal side of angle lies Quadrant" and the termnina ol angleE Quadrant III M cCs 0 = and sin 0 find the exact value of cos (a...
##### A copper wire has a cross-sectional area of 8.20Ã—10^(âˆ’7)m^2. Itcarries a current of 1.67 A. Find the electric-field magnitude inthe wire.
A copper wire has a cross-sectional area of 8.20Ã—10^(âˆ’7)m^2. It carries a current of 1.67 A. Find the electric-field magnitude in the wire....
##### Exercise 7.17 An economist reports a set of parameter estimates, including the coefficient esti- mates βι...
Exercise 7.17 An economist reports a set of parameter estimates, including the coefficient esti- mates βι 1.0, β2-0.8, and standard errors s(A)-0.07 and s(A) 0.07. The author writes The estimates show that Bi is larger than B2. 2 (a) Write down the formula for an asymptotic 95% confid...


text:
Question: CalcNormFactors in Edge r
0
6 months ago by
phickner20 wrote:

Can anyone explain the "CalcNormFactors" normalization process in edge r?  I have RNAseq data from insect antennae. I do not have a genome assembly for my organism, but I did assemble the transcriptome using trinity. I then retrieved my transcripts of interest (chemosensory genes) from the de novo assembly. I then mapped my reads to ~300 chemosensory genes. I used htseq count to count the reads. I provided the library sizes to edge r. I am concerned that since my counts represent a small proportion of all of the reads, normalization using CalcNormFactors may not be appropriate. I can't seem to find out exactly how the normalization is performed. Any help would be appreciated.

Thanks,

Paul

edger • 187 views
modified 5 months ago by Gordon Smyth37k • written 6 months ago by phickner20
2
6 months ago by
United States
James W. MacDonald50k wrote:

The short answer is 'don't do that'.

One issue with RNA-Seq is differences in library size. If you took two aliquots of the same exact sample and sequenced twice, one time getting 20M reads, and the other getting 10M reads, the expectation for the counts/gene would be that you would get twice the number of counts for each gene in the first aliquot than the second (because you have twice the reads). You want to account for those differences, because in this silly example you already know that the genes were expressed at the same exact level in both samples, and the difference in counts is due only to the library size.

It's not that simple when you are using biological replicates. All things equal, if you have twice the reads in one sample than the other, you expect about twice the counts/gene. But all things aren't equal, and you may have some genes that are really highly expressed in one sample, which may contribute more than expected to the total library size (this is called compositional bias), in which case you would want to find a set of genes that aren't obviously changing expression bigly between the two samples. The TMM normalization that is the default in calcNormFactors is intended to do that - estimate normalization offsets for each library, based on a reasonable set of genes.

If you subset to some arbitrary set of genes, you may well lose all the genes that aren't really changing, and that would be useful for calculating differences in library size. In addition, all the other genes that you may not be interested in are quite useful for estimating things like the relationship between gene counts and variance (which you need for linear modeling). So in general you should keep all the genes that appear to be expressed at a reasonable level in your data set all the way through the model fit. After that, you can filter out the genes you care about.

I don't know how to go from a de novo transcriptome assembly that has over 400,000 contigs to a transcriptome that contains ~18,000 transcripts in its entirety. In my opinion, mapping back to the Trinity output will not produce reliable data. Furthermore, how do I generate a gff for an unknown transcriptome? My gene families of interest are highly divergent, often with many paralogs, and require manual annotation, so the typical pipelines are usually of little use. Any suggestions would be appreciated.

0
5 months ago by
Gordon Smyth37k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth37k wrote:

I wouldn't run calcNormFactors for your data. You can run an edgeR analysis on your 300 genes, omitting the calcNormFactors step. However you have to be sure to set the library sizes to the total number of reads you have for each samples, including reads not mapped to any of your 300 chemosensory genes.

For example, you might use

y <- DGEList(counts=counts, lib.size=n)

where counts is your count matrix with 300 rows and n is a vector giving the sequencing depth for each sample.

Compared to a standard RNA-seq analysis with a reference genome, you will lose some of the protection that calcNormFactors would normally give you, but the DE analysis will still be ok provided your RNA samples are of a consistent quality.


text:
# Thread: A trig. question

1. ## A trig. question

For range $0<=A<=360$, solve:
$3cos^2 2A=1+sin A$

Much help is appreciated. Thanks

2. Originally Posted by smmxwell
For range $0<=A<=360$, solve:
$3cos^2 2A=1+sin A$

Much help is appreciated. Thanks
Substitute $\cos^2 A = 1 - \sin^2 A$ and re-arrange into an equation that looks something like the one in this thread: http://www.mathhelpforum.com/math-he...sin-x-3-a.html

That thread tells you what to do next.

3. Yeah. I know about that. Just that it's sin^2 2A. Even if i use the double angle formula, I still can't get anywhere. I know about factoring a q.e. but I just can't rearrange it.

4. Originally Posted by smmxwell
Yeah. I know about that. Just that it's sin^2 2A. Even if i use the double angle formula, I still can't get anywhere. I know about factoring a q.e. but I just can't rearrange it.
My mistake. I misread the question. Are you expected to get exact solutions?

5. Yeah. The answer is A= 30, 150, 221 deg. 49', 318 deg. 11'.
I did my own steps and arrived at
$3cos^2 2A=(cos^2 A)/(1-sin A)$

Even though $cos^2 2A$ is definitely not equal to $cos^2 A$, still, if I eliminate them, I manage to find $sin A= -2/3$, which would yield the last two readings. So, any help?

btw, I just came back from a long break, so I can't really think straight now. I found another question that needs help too.

Solve $tan x + 2 sin x=0$ for range=1 circle.
How do you go further from $sin x +sin 2x= 0$

6. Originally Posted by smmxwell
Yeah. The answer is A= 30, 150, 221 deg. 49', 318 deg. 11'.
I did my own steps and arrived at
$3cos^2 2A=(cos^2 A)/(1-sin A)$

Even though $cos^2 2A$ is definitely not equal to $cos^2 A$, still, if I eliminate them, I manage to find $sin A= -2/3$, which would yield the last two readings. So, any help?

btw, I just came back from a long break, so I can't really think straight now. I found another question that needs help too.

Solve $tan x + 2 sin x=0$ for range=1 circle.
How do you go further from $sin x +sin 2x= 0$
I will solve a question is near to your see how I solve it

$\sqrt{3}tan x =\frac{sinx }{cosx}$

$\sqrt{3} \frac{sinx}{cosx} + 2sinx =0$

$\sqrt{3}sinx + 2sinx cos x =0$

$sinx(\sqrt{3}+2cosx)=0$

$sinx=0......x=n\pi$

$\sqrt{3}+2cosx=0......cosx=-\frac{\sqrt{3}}{2}......x=....$

your's is like this

7. I'm sorry but I don''t really follow. You mean you're giving me a similar question which has the same steps? Either way, I don't see how you get $
\sqrt{3} \frac{sinx}{cosx} + 2sinx =0
$
after $
\sqrt{3}tan x =\frac{sinx }{cosx}
$

the 2 question I posed weren't related. they are different questions which I both don't know how to do.

8. Ok I do a mistake in my solution

${\color{red}\sqrt{3}tan x =\frac{sinx }{cosx}}$ I forgot to multiply the answer with $\sqrt{3}$

${\color{blue}\sqrt{3}tan x =\sqrt{3}\left(\frac{sinx }{cosx}\right)}$ the correct one

I will solve two of them Step by Step so look

$\sqrt{3}tan x + 2 sin x =0$

$\sqrt {3} \left(\frac{sinx}{cosx}\right) + 2 sinx =0$

$\sqrt{3} sinx + 2 sinx cosx =0$

$sinx ( \sqrt{3} + 2 cos x ) =0$

sinx =0 so x=n(pi)
cosx = -square(3)/2

your's

$tanx + 2 sinx =0$

$\frac{sinx}{cosx} + 2sinx =0$

$sinx + 2 sinx cosx =0$

$sinx (1+ 2cosx ) =0$

$sinx =0 ......$

$1+2cosx =0.....cosx=\frac{-1}{2}$

it is ok or not

9. Yes. I get the 2nd one. But I still don't understand the first one;

$
3cos^2 2A=1+sin A
$

10. Originally Posted by smmxwell
Yes. I get the 2nd one. But I still don't understand the first one;

$
3cos^2 2A=1+sin A
$
I will not solve it completely I will give you hints just

$3cos^2(2A)=1+sinA$

you should know that

$cos(2A)=1-2sin^2(A)$ find the

$(1-2sin^2(A))^{2}$

then plug it in your equation you will have an equation

you will have

$3(1-4sin^2(A)+4sin^4(A))=1+sin(A)$

$3-12sin^2(A)+12sin^4(A)=1+sin(A)$

rearrange it and let t=sin(A) you will have

$12t^4-12t^2-t+2=0$ this for you

11. ## Trig equations

Hello smmxwell
Originally Posted by smmxwell
For range $0<=A<=360$, solve:
$3cos^2 2A=1+sin A$

Much help is appreciated. Thanks
If the answers you gave are correct then there's a mistake in the question. None of them satisfies $3\cos^2 2A=1+\sin A$. (Try it and see!) However:

$3\cos2A = 1+\sin A$

$\Rightarrow 3(1-2\sin^2A)=1+\sin A$

$\Rightarrow 6\sin^2A +\sin A -2 = 0$

$\Rightarrow (2\sin A -1)(3\sin A +2)=0$

$\sin A = 0.5, -0.6667$

$A = 30^o, 150^o, 211.81^o, 318.19^o$

Grandad

12. Yeah thats the anwser. I was trying to type in degrees and minutes.
$
\sin A = 30^o, 150^o, 211.81^o, 318.19^o
$

@Amer. I don't think quartic eq. is in my syllabus. I'm sure there's go to be a more elegant and simpler way?

@ Grandad,

The steps are for
$
3\cos2A = 1+\sin A
$

Whereas, the question wants $
3\cos^2 2A=1+\sin A
$

your answer is correct though..

13. Originally Posted by smmxwell
Yeah thats the anwser. I was trying to type in degrees and minutes.
$
\sin A = 30^o, 150^o, 211.81^o, 318.19^o
$

@Amer. I don't think quartic eq. is in my syllabus. I'm sure there's go to be a more elegant and simpler way?

@ Grandad,

The steps are for
$
3\cos2A = 1+\sin A
$

Whereas, the question wants $
3\cos^2 2A=1+\sin A
$

your answer is correct though..
so you gave me a wrong equation that is not my fault I solve

$3cos^2(2A)=1+sinA$ since you gave me it

but you did not gave me

$3cos(2A)=1+sinA$

in math usually there is more than one way to solve questions that refer to the solver of the problem

14. No. The question is $
3cos^2(2A)=1+sinA
$
and only that.

Ok. Now I understand Grandad's post. The question is wrong. Alright. That means the question intended should be $
3\cos2A = 1+\sin A
$
. Thanks a lot guys. Lol. I just spent the whole day trying to solve a wrong question.


text:
# Math Help - Complex Variable question

1. ## Complex Variable question

I am not sure if I am posting this in the correct section, please do let me know if I am in the wrong section.
I need help with the question below

Prove that $\overline{z}^k = \overline{z^k}$ for every integer k, where we assume that $z\not=0$ when $k$ is negative.

2. ## Re: Complex Variable question

This follows from the fact that $\overline{1/z}=1/\overline{z}$ and $\overline{z}^k = \overline{z^k}$ for positive $k$. The latter fact follows by induction from $\overline{z_1z_2}=\overline{z}_1\cdot\overline{z}_ 2$.


text:
## Monday, March 31, 2014

### 99 Clojure Problems – 27: Group the Elements of a Set into Disjoint Subsets

1. In how many ways can a group of 9 people work in 3 disjoint subgroups of 2, 3 and 4 persons? Write a function that generates all the possibilities.
2. Generalize the above predicate in a way that we can specify a list of group sizes and the predicate will return a list of groups.

Note that we do not want permutations of the group members; i.e. ((Aldo, Beat), ...) is the same solution as ((Beat, Aldo), ...). However, we make a difference between ((Aldo, Beat), (Carla, David), ...) and ((Carla, David), (Aldo, Beat), ...).

You may find more about this combinatorial problem in a good book on discrete mathematics under the term "multinomial coefficients".

### Example

(deftest p27-multinomial-coeffients
(is (=
'((("Peter" "Paul") ("Mary"))
(("Peter" "Mary") ("Paul"))
(("Paul" "Mary") ("Peter")))
(group [2 1] ["Peter" "Paul" "Mary"]))))

### Solution

The mulinomial coefficients are a generalisation of the binomial coefficients that we saw in the last problem. $${n \choose k_1, k_2, \ldots, k_m} = \frac{n!}{k_1!\, k_2! \cdots k_m!}$$ I used this to validate some aspects of my solution as you can see in the tests. Showing the correctness of your solution in a traditional unit test gets quite unwieldy when you have, as in this example, 1260 combinations to check.

My solution to part a) of the problem is based on the idea to reuse the enumeration of all k-combinations of n elements from problem 26. We then take all possible combinations of four people out of the group of nine and remove these people from the set of all people by using Clojures set functions. We continue to do that for the 3-combinations and the 2-combinations and add them into a result vector.

The repeated set creation/difference code is neither elegant nor performant but a pragmatic choice given that the combinations function was supposed to work with/return lists. Another area for improvement is the calculation of the last k-combinations. It is redundant as its result is just the remainder of the original input sequence.

The solution to part b) is a bit more involved but very closely related in terms of its structure. It's a recursive solution that starts with calculating the combinations for the first group and mapping the result to a recursive self-call. In each recursive call we remove the elements/persons already "used" from the input set. We also remove the k for which we just calculated the combinations from the list of group sizes. The whole operation is wrapped in a mapcat to flatten the result into just one level as you can see in the example.

## Sunday, March 30, 2014

### 99 Clojure Problems – 26: Generate the Combinations of K Distinct Objects Chosen From the N Elements of a List

"In how many ways can a committee of 3 be chosen from a group of 12 people? We all know that there are C(12,3) = 220 possibilities (C(N,K) denotes the well-known binomial coefficients). For pure mathematicians, this result may be great. But we want to really generate all the possibilities."

### Example:

(deftest p26-combinations
(is (= '((0 1) (0 2) (1 2)) (combinations 2 (range 3))))
(is (= 220 (count (combinations 3 (range 12))))))

### Solution:

Here is my thinking, which is inspired by SICPs counting change example: If you want to enumerate all combinations of 2 distinct objects chosen from the list (0 1 2) you could fix the first element to be the first element of the list. That leaves you with a list of two elements and one more object to choose. There are only two combinations for C(2, 1): the two remaining elements of the list. Join them to the first elements we have already chosen and we end up with (0 1) and (0 2). Do the same thing again: Enumerate all combinations of two distinct objects from the remaining list (1 2), there is just one: the list itself is the answer.

We can generalize these observations into a solution: combinations of k distinct objects chosen from a list of n elements are the first element of the list conjoined to the combinations of k-1 distinct objects from the rest of the list, plus all combinations of k distinct objects from the list without its first element (again: the rest or tail of the list).

## Wednesday, March 19, 2014

### Example:

I am a bit at a loss when it comes to providing a meaningful test case that illustrates the problem, as I did in the previous posts in this series.

Based on my current understanding we would have to show that the randomness of the results matches a uniform random distribution. But lets look at the solutions first. I will come back to that question at the end.

### Solution:

There are at least three solutions to this problem that I am aware of.

#### 1. Idiomatic/Java Interop

The first one is to reuse the idiomatic variant of the random-select, which is based on java.util.Collections#shuffle, which in turn implements AFAIK the well known Fisher-Yates shuffle algorithm. While a functional purist might sniff at the mutation involved I think, this is perfectly fine in this case as the mutability is encapsulated in a functional/immutable shell. Your sequence is translated into a mutable collection to use the efficient algorithm based on swapping. The result is then transformed back into a persistent data structure: a pattern often used to get the best of both worlds.

#### 2. Naive Functional

Solution number two is functional but not very efficient. If I am not mistaken the running time of this solution should be O(n^2). We are calling remove-at n times, which calls split-at. Split-at has linear complexity.

#### 3. Perfect Functional Shuffle

The third solution is nothing I came up with myself. I merely followed a hint given by Phil Gold. It pointed to Oleg Kiselyov's post about perfect functional shuffle in Haskell.

Implementing it in Clojure turned out to be quite a lot of work. I have never done any Haskell and I am not sure if I got all of if right. It is probably best you have a look at Oleg's text first.

I did not read up much on Haskell syntax, but nonetheless: it started to make sense after looking at the code for a while even without knowing the language. The biggest difference with effects on the structure of the program: Haskell comes with continuation passing style built right in— something I could not or did not want to recreate in Clojure.

The gist of the solution is to build a binary tree out of the input sequence. The tree allows for efficient extraction of branches and is by its very nature at most ceil(log2 n) deep. The overall complexity of the solution should be therefore O(n*log n). Again, read the original text for more detail.

I decided to use core.match and add with it the first dependency to the project. This allowed me to stick closer to the original solution which makes heavy use of Haskell's pattern matching.

Translating the essential two functions was pretty straightforward and core.match offered all I needed. The syntax is a bit clunkier (notably the guard clauses) than Haskell's succinctness but it's not in a different ballpark alltogether.

The first interesting function is called pair-leaves. It takes a sequence of nodes—initially all leaves—(represented by vectors of the form [:node|:leaf count value]) and joins them together to form a binary tree. We need four pattern matches to do that: matching two leaves, a node and a leaf, a leaf and a node and finally one matching two nodes.

This gives us the binary tree. We now take a sequence of random numbers (each being an independent sample of a uniform random distribution [0...n-1]) and extract the numbers at the positions indicated by those random indices from the tree.

We need six pattern matches to extract from the tree. (Oleg needs five as Haskell apparently has something like a multiple guard with sub-branches which spares him writing the last match.) The patterns are:

1. The first element is always a leaf on the left.
2. Any but the first element when our tree consists of two leaves: the right leaf.
3. Any but the first element where our tree is bigger than two leaves: extract recursively to the right.
4. Any but the first element where n + 1 is equal to (or greater) than the count of the root node and we have a leaf as the right subtree. This is actually the biggest deviation from Oleg's code and possibly a mistake on my part.
But as I understand it, it is due to the lack of CPS in my solution. Oleg's continuations take the flow of execution back up into the shuffle function where he deals with the case of a single leaf. My solution has to be able to deal with the leaf case inside the extract-tree — by weakening the guard here.
5. The last two are closely related: a tree with a non-leaf sub-tree to the left and none of the above: if the number we are looking for is in the left subtree (as indicated by its count on the node) we recursivly extract to the left, otherwise to the right.

### Verifying the Results

What is left to do is to verify the property of uniform random distribution. The assumtions about time complexity as outlined by Oleg seem obvious to me. Let's look at a histogram:

This histogram left me wondering whether I had made a mistake. It does seem to show a quite clearly discernible repetitive sawtooth pattern. I created the data with this snippet of Clojure:

(def input '[a b c d e])
(def perms (vec (combo/permutations input)))
(def samples (repeatedly 10000 (fn [] (.indexOf perms (random-permute-functional input)))))

Compare this with a histogram of 10000 results generated by the first solution based on Fisher-Yates:

The most likely explanation is a mistake in my implementation or my reasoning. I have not found it yet, but I am quite keen on any hints. Could it be that perfect functional shuffle is not quite perfect? Probably not.

### Update:

Of course not. Andy Fingerhut was so kind to point out my mistake to me: You have to choose the random numbers uniformly distributed within [0,n-1] for the first number, [0, n-2] for the second number and so on. Otherwise you will bias your shuffle towards whatever is the rightmost branch of your tree, as every number > n where n is the size of your current tree will extract the rightmost element. You can find the corrected version here. The histogram looks much better now too: