# Modeling of the ALMA HL Tau Polarization by Mixture of Grain Alignment and
Self-scattering 111Released on June, 10th, 2019

Tomohiro Mori The Institute of Astronomy, the University of Tokyo, 2-21-1
Osawa, Mitaka, Tokyo 181-8588, Japan  Akimasa Kataoka National Astronomical
Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

(Received June 1, 2019; Revised January 10, 2019)

###### Abstract

Dust polarization at (sub)millimeter wavelengths has been observed for many
protoplanetary disks. Theoretically, multiple origins potentially contribute
to the polarized emission but it is still uncertain what mechanism is dominant
in disk millimeter polarization. To quantitatively address the origin, we
perform radiative transfer calculations of the mixture of alignment and self-
scattering induced polarization to reproduce the 3.1 mm polarization of the HL
Tau disk, which shows azimuthal pattern in polarization vectors. We find that
a mixture of the grain alignment and self-scattering is essential to reproduce
the HL Tau 3.1 mm polarization properties. Our model shows that the
polarization of the HL Tau at 3.1 mm can be decomposed to be the combination
of the self-scattering parallel to the minor axis and the alignment-induced
polarization parallel to the major axis, with the orders of $\sim 0.5\%$
fraction for each component. This slightly eases the tight constraints on the
grain size of $\sim 70~{}\micron$ to be $\sim 130~{}\micron$ in the previous
studies but further modeling is needed. In addition, the grain alignment model
requires effectively prolate grains but the physics to reproduce it in
protoplanetary disks is still a mystery.

polarization — protoplanetary disks — stars: individual (HL Tau) — techniques:
interferometric

††journal: ApJ††facilities: ALMA

##  1 Introduction

The presence of dust polarization of protoplanetary disks at
submillimeter/millimeter wavelengths have been indicated for several
protostars and T Tauri stars with the James Clerk Maxwell Telescope (JCMT)
single-channel polarimeter (Akeson & Carlstrom, 1997; Tamura et al., 1995,
1999). However, for two decades, spatially resolved polarization has not been
obtained due to the lack of the sensitivity and spatial resolutions of the
instruments (Hughes et al., 2009, 2013) until the first detection of the
millimeter polarization for the HL Tau disk by (Stephens et al., 2014).
Nowadays, ALMA dust polarization observations have enabled us to obtain the
polarization for various class II protoplanetary disks, showing the various
polarization morphology and polarization fraction (Kataoka et al., 2016b,
2017; Stephens et al., 2017; Hull et al., 2018; Bacciotti et al., 2018; Ohashi
et al., 2018; Dent et al., 2019; Harrison et al., 2019; Mori et al., 2019).

One difficulty in interpreting the disk polarization is that multiple
mechanisms of polarization are at work at millimeter wavelengths. One
mechanism is the grain alignment, where elongated dust grains are aligned to
some directions and emit polarized thermal dust emission. The grain alignment
mechanism has been long and well studied, especially the radiative torques
(RATs) (e.g., Lazarian & Hoang, 2007). Here we only discuss which direction
the grains are aligned to because it reflects the polarization orientation.
The directions of the grain alignment may represent magnetic fields (Davis &
Greenstein, 1951; Cho & Lazarian, 2007), radiation (Tazaki et al., 2017), or
gas flow (Gold, 1952; Lazarian & Hoang, 2007; Kataoka et al., 2019). Since the
grain alignment has been mainly studied in the condition of the interstellar
medium, further discussion of the alignment physics in the protoplanteary
disks is necessary. Another mechanism that would produce polarization is the
self-scattering, where the dust grains scatter and polarize the incoming
thermal emission from other dust grains (Kataoka et al., 2015). Dust grains
with the size comparable to observed wavelengths efficiently scatter and
polarize the thermal emission. The interpretations of the observations depend
on each target. However, the interpretation itself is one of the major issues
of the millimeter polarization of protoplanetary disks.

To interpret the dust polarization, it is essential to disentangle the
polarization mechanisms. In this paper, as one step toward complete
understanding of the polarization, we investigate the effects of superposition
of the two major polarization mechanisms, which are the self-scattering and
alignment. Especially, we model the 3.1 mm polarization of the HL Tau disk,
which has been interpreted purely as alignment-induced polarization (Kataoka
et al., 2017), but also pointed out to be contaminated by self-scattering
(Yang et al., 2019).

HL Tau is a class I/II protostar in the Taurus star-forming complex at a
distance of 140 pc from the solar system (Rebull et al., 2004; Robitaille et
al., 2007). A notable feature of the surrounding disk is the concentric multi-
ringed structure with a spatial scale of $\sim$100 au (ALMA Partnership et
al., 2015). The polarization morphology at millimeter wavelengths on the HL
Tau disk shows strong wavelength dependence (Stephens et al., 2014; Kataoka et
al., 2017; Stephens et al., 2017). The polarization morphology at 870
$\micron$ is parallel to the disk minor axis, while that at 3.1 mm shows
almost the azimuthal pattern. The polarization pattern with the vectors
parallel to the minor axis is well reproduced by the self-scattering model
(Yang et al., 2016; Kataoka et al., 2016a). Therefore, the $870~{}\micron$
polarization has been interpreted by the self-scattering. In contrast, the
origin of the 3.1 mm polarization has been interpreted by dust alignment
because the self-scattering fails to explain the azimuthal pattern. It has
been first interpreted by grain alignment with the radiative anisotropy
(Tazaki et al., 2017; Kataoka et al., 2017), but slight inconsistency in the
morphology has also been pointed out (Yang et al., 2019).

The switching of the polarization mechanisms between 870 $\micron$ and 3.1 mm
must be explained. The reason why the self-scattering disappear at 3.1 mm is
an expected behavior because polarization fraction quickly drops if dust
grains are larger than the wavelengths divided by $2\pi$ (Kataoka et al.,
2015). However, this requires the grain size to be $\sim$ 100 $\micron$, which
is not consistent with the general SED analysis, which requests millimeter to
centimeter dust grains (e.g., Beckwith & Sargent, 1991; Ricci et al., 2010).
Carrasco-González et al. 2019 analyzed the millimeter SED of the HL Tau disk
with the effects of scattering-induced intensity reduction (Liu, 2019; Zhu et
al., 2019) and obtained grain size to be millimeter, which is larger than
$\sim$ 100 µm. In contrast, the reason why alignment-induced polarization are
observed at 3.1 mm but disappear at 870 $\micron$ is still unresolved. One
idea to explain the transition is that the grain size is comparable to the
wavelengths, and their intrinsic polarization is strong enough to be detected
at 3.1 mm but too weak at 870 $\micron$. This happens because of the
transition from Rayleigh to Mie regime. If the grain size is much smaller than
the wavelengths, the intrinsic polarization does not have strong wavelength
dependence. If the grain size is much larger than the wavelengths, in
contrast, the polarization fraction drops. Note that if the grain size is
comparable to the wavelengths, polarization direction becomes perpendicular to
the grain major axis because of the Mie regime (e.g., Guillet et al., 2020).
In this way, the wavelength dependence between millimeter and submillimeter
wavelengths are the key to understanding the polarization.

Toward the full understanding of the wavelength-dependent polarization of the
HL Tau, in this paper, we focus on the modeling the 3.1 mm polarization by the
mixture of the self-scattering and grain alignment. Yang et al. (2019) raised
the issue that neither of aerodynamic nor radiative flux induced polarization
can explain the 3.1 mm polarization and speculated that possible contamination
of the self-scattering component. We quantitatively investigate how strong the
self-scattering contamination is even at the 3.1 mm band and leave the full
modeling of the three wavebands for future studies. Note that this does not
solve the issue that no alignment model can perfectly explain the
polarization, but it may change the wavelength dependence of the self-
scattering components in the HL Tau polarization, which may have an impact on
the grain size constraints.

Figure 1: ALMA Band 3 ($\lambda$ = 3.1 mm) polarimetric observations of the HL
Tau disk, which has been previously reported by Kataoka et al. 2017 and
Stephens et al. 2017. The left panel shows the Stokes I intensity in color,
the central panel shows the polarized intensity, and the right panel shows the
polarization fraction. The solid contours in the three panels are equivalent
and represent the intensity with the levels of (10, 25, 60, 150, 366, 900)
times the noise level, which is $\sigma_{I}$ (= 34 $\mu$Jy beam-1). The
overlaid line segments represent the polarization vectors where the polarized
intensity is larger than 3 times the noise level of the polarized intensity,
which is $\sigma_{PI}$ (= 5.3 $\mu$Jy beam-1). Note that the length of the
polarization vectors is set to be equal. The beam has a size of 0.″51 $\times$
0.″41 and position angle of -6$\fdg$76, which is represented in the bottom
left as an ellipse.

##  2 Model

We briefly describe the outline of the radiative transfer models in this
study. First, we construct a dust disk model that reproduces the Stokes I
emission obtained with the 3.1 mm polarization observation. Then, on the
constructed disk model, we take into account the effect of the grain alignment
and self-scattering in calculating polarization.

###  2.1 Disk Model

Figure 1 shows the Stokes I intensity, polarized intensity, and polarization
fraction, all of which are overlaid with the polarization vectors. The notable
features of the 3.1 mm polarization are (1) azimuthal polarization morphology
and (2) azimuthally uniform polarization fraction at $\sim$ 1%.

We construct an axisymmetric disk model that reproduces the stokes I image
obtained on the polarization observation at 3.1 mm (Kataoka et al., 2017;
Stephens et al., 2017). By assuming that the Stokes I emission at 3.1 mm is
optically thin, we construct the radial profile of the optical depth at 3.1 mm
with the form of

$\rm\tau_{model}(r)=\tau_{0}\>\Bigl{(}\frac{r}{\rm
1\>au}\Bigr{)}^{-\gamma_{1}}\exp\Bigl{[}\Bigl{(}\frac{r}{r_{c}}\Bigr{)}^{-\gamma_{2}}\Bigr{]},$
(1)

with the temperature profile of $\rm T(r)=310(r/1AU)^{-0.57}$ K (Okuzumi &
Tazaki, 2019). We note that both the dust surface density and temperature are
assumed to have smooth radial distributions without any ring and gap
structures observed with the higher spatial resolution (ALMA Partnership et
al., 2015) for simplicity.

We adopt the parameter set of ($\rm\tau_{0},r_{c},\gamma_{1},\gamma_{2}$) =
(1.0, 129.2, 6.6, 9.0), which reproduces the Stokes I emission. Figure 2 shows
the images of observed and modeled Stokes I produced assuming the parameter
set and residual.

We derive the dust surface density profile $\rm\Sigma_{d}(r)$ by calculating
$\rm\Sigma_{d}(r)$ = $\rm\tau_{model}/\kappa_{abs}$, where $\rm\kappa_{abs}$
is the absorption opacity. As we will describe below, the value of
$\kappa_{abs}$ is a function of $a_{max}$ and $\lambda$ when constituent
materials are determined. For the fiducial value, we adopt
$\rm\kappa_{abs}(a_{max}=100\>\micron,\lambda=3.1\>mm)$ = 0.0788 cm2 g-1 in
calculating the dust surface density.

The vertical dust density is assumed to be Gaussian density distribution with
a dust scale height $h_{d}$ such that
$\rm\rho_{d}=\Sigma_{d}/\sqrt{2\pi}h_{d}\exp(-z^{2}/2h_{d}^{2})$. We also
assume that $h_{d}$ is same as the gas pressure scale height.

While the residuals are still significant, they show axisymmetric
distribution. For this study, we aim at reproducing the continuum emission
with an axisymmetric model for simplicity, and thus we do not further search
the parameter sets.

Figure 2: The left panel shows the Stokes $I$ image of the observations, the
central panel shows the model image, and the right panel shows the residual
map where the observed Stokes $I$ subtracted with the model image. The
contours of the left and the central panels are the same as Fig.1 and the
contour of the right panels show the residual levels of (-10, -25) $\times$
$\sigma_{I}$

.

###  2.2 Dust Model

We use a grain composition model developed by Birnstiel et al. (2018), where
the grains are the mixture of silicate, troilite, organics, and water ice. The
refractive index used in the calculation is as follows: Draine (2003) for
astronomical silicate, Henning & Stognienko (1996) for troilite and refractory
organics, and Warren & Brandt (2008) for water ice. We compute the mixture of
them with the effective medium theory using the Maxwell-Garnett rule (e.g.,
Bohren (1983); Miyake & Nakagawa (1993)). We assume the grains have a power-
law-size distribution with a power of $q$ = -3.5 (Mathis et al., 1977) with
the maximum grain size of $a_{max}$, while the minimum grain size is fixed at
$a_{min}$ = 0.05 $\mu$m. We assume the same grain population throughout the
disk.

To reproduce polarization, we consider both alignment and scattering-induced
polarization. For the alignment-induced polarization, we consider both prolate
and oblate grains. 222Grains are believed to be spinning when they are
aligned. The spinning motion makes the grains to be effectively prolate or
oblate ellipsoidal bodies. We call prolate/oblate grains for simplicity but
this does not rule out any complex shapes of the dust grains. The motivation
to consider the both models is the observed polarization morphology. There are
two ways to reproduce azimuthal polarization vectors; effectively prolate
grains with their major axis to be parallel to the disk azimuthal direction
and effectively oblate grains with their minor axis parallel to the disk
radial direction. We assume the Rayleigh regime (Lee & Draine, 1985) to
calculate the intrinsic polarization. This approximation holds when the size
parameter $x=2\pi a/\lambda\ll 1$. We assume the grain size up to 200
$\micron$ at the wavelength of 3.1 mm, which corresponds to the size parameter
of $x=0.42$, which justifies the use of Rayleigh regime. We assume the axes
ratio of the grains to be $\alpha=0.9$ for prolate grains and $\alpha=1.1$ for
oblate grains and that the alignment efficiency to be $\epsilon=0.6$. We
assume that the oblate grains are aligned with their short axis parallel to
the radial direction of the disk while the prolate grains are aligned with
their long axis parallel to the azimuthal direction of the disk. We ignore the
effects of temperature for simplicity.

As pointed out by Guillet et al. (2020), to discuss the wavelength dependence
of the alignment-induced polarization, the Mie regime is essential to be
considered. In this paper, however, we do not treat the Mie regime for
simplicity and we focus on the single band observations on the HL Tau at 3.1
mm. In the cases of further modeling at 870 $\micron$ and at 1.3 mm, the Mie
regime treatment is required because the size parameter of $200~{}\micron$
grains becomes $x=1.44$ at 870 $\micron$ and $x=0.97$ at 1.3 mm.

The scattering-induced polarization is included by turning on the scattering
in the radiative transfer simulations. We treat the full phase function and
take into account the multiple scattering. The controlling parameter of the
scattering is the maximum grain size, $a_{max}$. The self-scattering
polarization is the most effective when the maximum grain size is $a_{max}$
$\sim$ $\lambda/2\pi$ (Kataoka et al., 2015).

We assume spherical dust grains for computing self-scattering while we assume
oblate or prolate grains for computing alignment-induced polarization. These
are not consistent but we neglect the effects of scattering by non-spherical
grains for simplicity. The angle-averaged absorption opacity of oblate/prolate
grains is set to be the same as the spherical grains.

We also note a technical treatment. The radiative transfer code RADMC-3D
allows users to compute the alignment-induced polarization only when the dust
scattering is turned on. To reproduce the alignment-induced polarization
without scattering, therefore, we use $a_{max}$ = 50 $\mu$m with the
scattering turned on, where the scattering-induced polarization is negligible.

##  3 Results

###  3.1 Alignment-only Model

First, we see how the beam dilution affects the appearance of the polarization
vectors. We use the same disk model of HL Tau but assume that the disk is
observed with a face-on view, $i=0^{\circ}$. The left panels of Figure 3
demonstrates the effects of beam dilution with the face-on view. To reproduce
the azimuthal polarization pattern, the prolate grains are assumed to be
aligned with their long axis to be parallel to the azimuthal direction, while
the oblate grains are with their short axis to be parallel to the radial
direction. Both models generally shows the azimuthal pattern in polarization
vectors, but the vectors at around the central region have certain inclination
angles, which is not expected in the case without beam dilution.

Next, we see the effects of changing the disk inclination angle on the
appearance of the polarization vectors and polarized intensity. The center and
right columns of Figure 3 show the case of the inclination angles of 45∘, and
60∘. The polarization vectors are almost in the azimuthal directions in the
inclination angles of 45∘ and 60∘ for both the prolate and oblate grains,
except the central regions where the beam dilution affect the morphology.

Yang et al. (2019) proposed that, if the disk inclination is non-zero, we can
distinguish the alignment models between effectively prolate and oblate grains
by investigating if the polarization vectors are in the azimuthal directions
(i.e., azimuthal pattern) or are with their normal direction to be toward the
disk center (i.e., circular pattern). However, as we see in Figure 3, the
difference between elliptical and circular patterns is hard to be
distinguished because of the beam dilution.

In contrast, the polarized intensity shows prominent differences between the
prolate and oblate grain models; the prolate grain models show higher
polarized intensity (i.e., brighter color in Figure 3) along the minor axis
than the major axis, while the oblate models show higher polarized intensity
along the major axis. The reason can be understood by a simple geometry.
Prolate grains along the disk minor axis always shows the longest axis in the
observed projected plane even if the disk inclination is changed, while those
along the disk major axis shows less polarization fraction if the disk
inclination changes to be closer to edge-on view because of the projection
effects. The situation is completely opposite in the case of oblate grains.
Oblate grains along the disk major axis always shows the longest axis to the
observers while those along the minor axis becomes fainter in higher
inclination angles because of the projection effects. Therefore, investigating
the polarized intensity pattern is a promising way to distinguish the
mechanisms between prolate and oblate models if the polarization is purely due
to the grain alignment.

Now we come back to the case of the HL Tau disk. As shown in the central panel
of Figure 1, the polarized intensity of the HL Tau disk at 3.1 mm wavelength
shows little variation in the azimuthal direction. This is different from the
models that we discussed, where the prolate grain model shows that polarized
intensity is brighter along the disk minor axis while the oblate model shows
it is brighter along the disk major axis. This illustrates that alignment-only
models may not be able to reproduce the observed polarized intensity pattern.
This motivates us to add another mechanism of polarization, self-scattering,
which we will discuss in the next chapter.

We note that the general properties described above hold for the polarization
fraction as shown in Figure 5. In addition, the general behavior and
differences between prolate and oblate grain models is also discussed by Yang
et al. (2019) with semi-analytical models.

Figure 3: The six panels show the polarized intensity of models with prolate
and oblate grains with different disk inclinations. The top three panels show
the results of the prolate grains with their long axis parallel to the
azimuthal direction and the disk inclinations are $i$ = 0∘ (left), 45∘
(center), and 60∘ (right). The bottom three panels show the same as the top
but for oblate grains. The levels of contours and the line segments are the
same as Fig.1 but for each model. The synthesized beam is represented as the
ellipse at the bottom left of each panel.

###  3.2 Mixture of Alignment and Scattering Models

Next, we discuss a possibility that a mixture of the alignment and self-
scattering models can reproduce the HL Tau polarization at 3.1 mm. The self-
scattering polarization generally produces the polarization vectors parallel
to the minor axis, and the fraction is a strong function of the grain size.
Therefore, we set the maximum grain size $a_{max}$ to be the the major
parameter and changes it to reflect the strength of the self-scattering
polarization. We use the same disk models of the alignment-induced
polarization of the prolate and oblate cases, and add the self-scattering-
induced polarization.

Figure 4 shows the polarized intensity of the prolate and oblate grain models
with $i$ = 45∘. The contributions of the self-scattering model increases by
increasing $a_{max}$ from left to right panels while that of the alignment
models are fixed.

As expected, the oblate + scattering models do not reproduce the observation.
Both polarization of the alignment of oblate grains and self-scattering
produces stronger polarization along the major axis. As a result of
superposition, the azimuthal variation becomes stronger, which is not
consistent with the observation as shown in Figure 1, which shows little
azimuthal variation in polarized intensity. Instead, the prolate + scattering
model shows less azimuthal variation than the oblate + scattering model. This
is because the prolate grain alignment produces strong polarization along the
disk major axis while the self-scattering does along the disk minor axis. As a
result of their superposition, the polarized intensity shows less azimuthal
variation, which is consistent with the observations.

Now the question is how strong is the contributions of the alignment and
scattering. We focus on the prolate + scattering model and discuss it. The top
three panels of Figure 4 show the case of different grain sizes, which
corresponds to difference of contribution of self-scattering. As demonstrated
in the figure, the $a_{max}$ = 100 $\mu$m case is dominated by the alignment-
induced polarization because the polarized intensity is stronger along the
minor axis, while the $a_{max}$ = 200 $\mu$m is by the self-scattering because
it is stronger along the major axis. The $a_{max}$ = 130 $\mu$m case shows the
half and half contributions between alignment and scattering, which reproduce
little azimuthal variation as observed. Therefore, we conclude that prolate +
scattering model with $a_{max}$ = 130 $\mu$m best reproduces the observed HL
Tau polarization at 3.1 mm. Note that we also successfully reproduce the
observed polarization fraction with the mixture model $a_{max}$ = 130
$\micron$ as shown in Figure 5.

Figure 4:  The six panels show the polarized intensity of models with mixture
of the grain alignment and self-scattering. The upper three panels show the
case of the combination of the aligned prolate grain and self-scattering. The
prolate grain model is fixed through the three panels while the self-
scattering-induced polarization is changed by setting the maximum grain size
to be $a_{max}$ = 100 $\micron$ (left), 130 $\micron$ (center), and 200
$\micron$ (right). Note that the contribution of self-scattering is weaker for
the left panel and stronger for the right panel. The bottom panels are the
same as the top panels but for oblate grains.

##  4 Discussion

The main results of the radiative transfer calculations are (1) the HL Tau 3.1
mm polarization can be reproduced with the mixture of the alignment and self-
scattering models and (2) prolate grain model can explain the observations but
oblate model does not. In this section, we discuss the implications on the
grain properties from both alignment and scattering theories as well as the
caveats of our models.

###  4.1 What can reproduce effectively prolate shape?

As discussed also by Yang et al. (2019), the effectively prolate grain model
is more preferable for reproducing the HL Tau polarization but it is
theoretically counter-intuitive. We usually assume that grains are spinning if
they are aligned, and they are observed as effectively oblate grains.

One possible speculation for explaining the prolate grain model is that grains
are not spinning in the HL Tau disk. The only mechanism that does not assume
grain spinning so far would be the mechanical alignment with a supersonic gas
flow, which is so-called Gold alignment, where the long axis of needle-like
grains is parallel to the gas flow onto dust grains (Gold, 1952). However, the
Gold alignment assumes that the gas speed on the dust grains is supersonic
whereas that is generally subsonic in protoplanetary disks (e.g., Cho &
Lazarian, 2007). All the discussion above is phenomenological and we need
further development of general grain alignment theory in protoplanetary disks.

Furthermore, we note that, even if the Gold mechanism is at work, the
predicted polarization pattern is not consistent with the observations. The
direction of the polarization would be parallel to the gas flow on the frame
of dust grains (Kataoka et al., 2019), which is a strong function of the
Stokes number, which denotes how well dust grains couple to ambient gas. The
grain radius $a$ is found to be $a_{max}$ = 130 $\micron$ in our modeling,
which gives the Stokes number to be St = 5 $\times$ 10-3 with the assumption
of the gas surface density of $\rm\Sigma_{g}$ $\sim$ 20 g cm-2. With such a
small number of St much smaller than unity, the directions of gas velocity on
grains would be in the radial directions (see Figure 2 in Kataoka et al.,
2019). Note that this discussion holds even for ring-like disks (Mori et al.,
2019). However, our modeling has found that the directions of the relative
velocity are dominantly in the azimuthal directions, which is incompatible
with the theory.

###  4.2 How Efficient is the Scattering at 3.1 mm?

The size of dust grains has been constrained by modeling the wavelength
dependence of the contribution of the scattering-induced polarization. In the
case of HL Tau disk, previous studies have assumed that the Band 7 (0.87 mm)
polarization is dominated by self-scattering while the Band 3 (3.1 mm) has no
contribution of self-scattering, which results in the tight constraint on the
grain size to be 70 $\micron$ in size (Kataoka et al., 2016a; Yang et al.,
2016; Kataoka et al., 2017). However, we have revealed that there is a
contribution of self-scattering even for the 3.1 mm polarization. Therefore,
we discuss if our modeling changes the interpretation of the grain size in the
HL Tau disk.

To estimate the scattering component at each wavelength, as conducted in
Kataoka et al. (2017), we compute the total polarization across the HL Tau
disk, which is given by

$\displaystyle PF_{total}=\left\\{\begin{array}[]{ll}1&\rm{(minor)}\\\
-1&\rm{(major)}\\\ \end{array}\right\\}\times{\sqrt{(\Sigma Q)^{2}+(\Sigma
U)^{2}}}{\Sigma I},$ (4)

where $\Sigma$I, $\Sigma$Q, and $\Sigma$U are the Stokes I, Q, and U emissions
integrated for the whole disk. Here, we explicitly define the sign of the
total polarization fraction; $PF_{total}$ is positive if the position angle of
the total polarization, which is calculated from $\Sigma$Q, and $\Sigma$U, is
closer to the disk minor axis, while $PF_{total}$ is negative if the position
angle is closer to the disk major axis. This is motivated to be able to
estimate the canceling out effects. In the case of axisymmetric protoplanetary
disks, the position angle of polarization is usually parallel to either the
disk minor or major axes because of the symmetry. If this is the case,
defining sign with respect to the disk minor and major axes helps to see which
components dominate the total polarization.

Table 1 summarizes the computed total polarization fraction for the
observations and models. We update the total polarization fraction presented
in Kataoka et al. (2017), who used the SMA and CARMA data, with the ALMA data
in Band 6 and 7 while its upper limit in ALMA 3 is presented the same as that
study. $PF_{total}$ is calculated from the final images as shown in Figure 4
with the equation (4). We also show $PF_{scat,model}$ and $PF_{align,model}$,
which represent the total polarization of the cases where either scattering or
alignment model is calculated. As shown in the table 1, the total polarization
fraction is nearly equal to the sum of those of the scattering and alignment
components. This illustrates the powerful estimation of the total polarization
by considering the positive and negative signs of each polarization
components.

Table 1:  Summary of the total polarization fraction of the models and observations. The first column represents the model names. $\lambda$ is the wavelengths. $PF_{total}$ represents the total polarization fraction, which can be decomposed to the scattering component, $PF_{scat}$ and the alignment component, $PF_{align}$. We define positive values of the polarization fraction for the component parallel to the disk minor axis while negative values for that parallel to the disk major axis.  |  $\lambda$ [mm] |  $PF_{total}(\simeq PF_{scat}+PF_{align})$ [%] |  $PF_{scat}$ [%] |  $PF_{align}$ [%]
---|---|---|---|---
oblate | 3.1 | 1.07 | 0 | 1.07
oblate + scattering ($a_{max}$ = 100 $\micron$) | 3.1 | 1.25 | 0.20 | 1.07
oblate + scattering ($a_{max}$ = 130 $\micron$) | 3.1 | 1.47 | 0.42 | 1.07
oblate + scattering ($a_{max}$ = 200 $\micron$) | 3.1 | 2.23 | 1.21 | 1.07
prolate | 3.1 | -0.55 | 0 | -0.55
prolate + scattering ($a_{max}$ = 100 $\micron$) | 3.1 | -0.35 | 0.20 | -0.55
prolate + scattering ($a_{max}$ = 130 $\micron$) | 3.1 | -0.12 | 0.42 | -0.55
prolate + scattering ($a_{max}$ = 200 $\micron$) | 3.1 | 0.68 | 1.21 | -0.55
observation | 3.1 | $-0.1<PF_{total}<0.1$ | $\\--$ | $\\--$
observation | 1.3 | $0.551\pm 0.003$ | $\\--$ | $\\--$
observation | 0.87 | $0.619\pm 0.004$ | $\\--$ | $\\--$

As shown in the table 1, oblate grain model always shows the positive value
while the prolate grain model shows the negative values, where positive value
corresponds to the polarization parallel to the disk minor axis while negative
value corresponds to that parallel to the disk major axis. In contrast, the
self-scattering always shows positive values because of the disk geometry. To
reproduce the tight upper limit on the absolute value of the total
polarization fraction to be less than 0.1 %, the oblate grain model is not
favored and the combination of the prolate model and self-scattering is the
key.

Here, we emphasize our big assumption that we do not know the wavelength
dependency of the alignment-induced polarization, and assume that there is no
contribution from alignment polarization at the wavelengths of the Band 6 and
7. To reproduce this, we need to calculate the intrinsic polarization in the
Mie regime (Guillet et al., 2020; Kirchschlager et al., 2019), which
computationally costs.

By assuming that there is no contribution of alignment induced polarization at
Bands 6 and 7, we confirm the decreasing trend of the observed total
polarization fraction from 0.87 mm to 3.1 mm, which has been discussed by
Kataoka et al. (2017) partially with the SMA and CARMA results from (Stephens
et al., 2014). Furthermore, our study revealed that the non-detection of total
polarization at 3.1 mm is explained by the combination of the alignment (-0.55
%) and the self-scattering ($+0.42$ %), which requires modification on the
wavelength-dependence of the self-scattering polarization. However, the
contribution of self-scattering to the total polarization at 3.1 mm, which is
estimated to be 0.42 % in this study, is still smaller than the 0.87 mm and
1.3 mm total polarization fractions. This means that the decreasing trend of
the self-scattering polarization still holds, which supports the idea of $\sim
100$ micron-size grains because the wavelength, at which the self-scattering
polarization peaks is, 0.87 mm or shorter (see Kataoka et al., 2016a, 2017;
Yang et al., 2016).

Still, there are significant differences in the size measurements between
millimeter-wave polarization and spectral indices. Latest modeling for the
continuum emission of the HL Tau disk obtained at multiple wavelengths from
0.9 mm to 13 mm found the radial gradient of the grain sizes from $a_{max}$ =
1.5 mm in the inner region to $a_{max}$ = 500 $\micron$ in the outer region,
which is at least four times larger than that in this study (Carrasco-González
et al., 2019). Plenty of differences between the two studies, for example, the
observation wavelengths, angular resolutions, and modeling methods can affect
the size measurements and partly explain the discrepancy. Another primary
difference is assumed constituent materials of the dust grains, which
determine refractive indices. Because the absolute values of the
absorption/scattering opacities and their wavelength dependence are strongly
sensitive to the assumption on the composition, this can affect the
measurement conducted in both of the studies (Testi et al., 2014). To
reconcile the discrepancy, further modeling that assumes the same grain
composition and disk parameters should be conducted.

###  4.3 Caveats

In the radiative transfer calculations, the dust models are not fully
consistent between the alignment and scattering models. When we calculate the
polarization of aligned elongated grains, we set an angle-dependent absorption
opacities. However, we assume spherical dust grains that have the same
absorption opacities for elongated grains for calculating their scattering
because of the difficulties in calculating scattering opacities of non-
spherical dust grains.

##  5 Conclusions

The origin of the polarized emission of the HL Tau disk at 3.1 mm has been
debated in the previous works. We have focused on the 3.1 mm polarization of
HL Tau, which shows the azimuthal pattern of polarization vectors, and have
conducted radiative transfer simulations for the 3.1 mm HL Tau polarization.
Our main findings are as follows.

  1. 1.

By radiative transfer simulations, we confirmed that alignment-induced
polarization cannot reproduce the 3.1 mm polarization feature of the HL Tau
disk as discussed by Yang et al. (2019) with a semi-analytical model.

  2. 2.

Next, we run models of the mixture of the grain alignment and self-scattering
and found that the combination of the effectively prolate grain emission and
self-scattering can reproduce the observations. However, the combination of
the effectively oblate grain emission and self-scattering cannot reproduce the
polarization. Our model shows that the total polarization fraction at 3.1 mm
can be understood by the combination of 0.42 % polarization due to the self-
scattering, which is parallel to the disk minor axis, and 0.55 % of the
alignment-induced polarization, which is parallel to the disk major axis.

  3. 3.

We found that the maximum size of the scattering dust grains is $a_{max}$
$\sim$ 130 $\micron$ to reproduce the azimuthally uniform polarization
fraction. This size is consistent with them obtained on the polarization
observation in Band 6 and 7. However, the size is a few times smaller than
that was measured by the multi-wavelength continuum analysis (Carrasco-
González et al., 2019). We confirmed that effectively prolate grain models are
more preferable than the oblate grain models to reproduce the 3.1 mm
polarization of HL Tau. However, the physical model to reproduce the situation
is still uncertain.

We note that while the radiative transfer calculations have been performed for
explaining the 3.1 mm polarization of the HL Tau disk, further modeling of the
multi-wavelength polarization with the results at 870 $\micron$ and 1.3 mm is
required in future studies. Modeling of the multiband polarization
observations will provide the spectrum of alignment and scattering efficiency,
which enable us to further constrain the grain sizes and shapes.

We gratefully thank Ian Stephens for providing us the analyzed ALMA Band 3
polarization data, which has enabled us to conduct the modeling work. We are
also deeply grateful to Takashi Miyata, Takafumi Kamizuka, and Ryou Ohsawa for
fruitful discussions for our study. This work is supported by JSPS KAKENHI
Numbers 18K13590 and 19H05088. ALMA is a partnership of ESO (representing its
member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST
and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the
Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO,
and NAOJ.

##  Appendix A The Polarization Fraction in the Mixture Models (Alignment +
Self-scattering)

Figure 5 shows polarization fraction mixture of the grain alignment and self-
scattering. As well as Figure 4 in Section 3, the model with prolate grains
and $a_{max}$ = 130 $\micron$ (top-central panel) successfully reproduces the
azimuthally uniform polarization fraction with $\sim$1$\\--$2 %, whereas other
models predict significant azimuthal variations.

Figure 5: Same as Figure 4 but for polarization fraction.

## References

  * Akeson & Carlstrom (1997) Akeson, R. L., & Carlstrom, J. E. 1997, ApJ, 491, 254, doi: 10.1086/304949
  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3, doi: 10.1088/2041-8205/808/1/L3
  * Bacciotti et al. (2018) Bacciotti, F., Girart, J. M., Padovani, M., et al. 2018, ApJ, 865, L12, doi: 10.3847/2041-8213/aadf87
  * Beckwith & Sargent (1991) Beckwith, S. V. W., & Sargent, A. I. 1991, ApJ, 381, 250, doi: 10.1086/170646
  * Birnstiel et al. (2018) Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45, doi: 10.3847/2041-8213/aaf743
  * Bohren (1983) Bohren, C. F. 1983, Tellus Series B Chemical and Physical Meteorology B, 35B, 65, doi: 10.1111/j.1600-0889.1983.tb00008.x
  * Carrasco-González et al. (2019) Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71, doi: 10.3847/1538-4357/ab3d33
  * Cho & Lazarian (2007) Cho, J., & Lazarian, A. 2007, ApJ, 669, 1085, doi: 10.1086/521805
  * Davis & Greenstein (1951) Davis, Leverett, J., & Greenstein, J. L. 1951, ApJ, 114, 206, doi: 10.1086/145464
  * Dent et al. (2019) Dent, W. R. F., Pinte, C., Cortes, P. C., et al. 2019, MNRAS, 482, L29, doi: 10.1093/mnrasl/sly181
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241, doi: 10.1146/annurev.astro.41.011802.094840
  * Gold (1952) Gold, T. 1952, MNRAS, 112, 215, doi: 10.1093/mnras/112.2.215
  * Guillet et al. (2020) Guillet, V., Girart, J. M., Maury, A. J., & Alves, F. O. 2020, arXiv e-prints, arXiv:2001.08400.  https://arxiv.org/abs/2001.08400
  * Harrison et al. (2019) Harrison, R. E., Looney, L. W., Stephens, I. W., et al. 2019, ApJ, 877, L2, doi: 10.3847/2041-8213/ab1e46
  * Henning & Stognienko (1996) Henning, T., & Stognienko, R. 1996, A&A, 311, 291
  * Hughes et al. (2013) Hughes, A. M., Hull, C. L. H., Wilner, D. J., & Plambeck, R. L. 2013, AJ, 145, 115, doi: 10.1088/0004-6256/145/4/115
  * Hughes et al. (2009) Hughes, A. M., Wilner, D. J., Cho, J., et al. 2009, ApJ, 704, 1204, doi: 10.1088/0004-637X/704/2/1204
  * Hull et al. (2018) Hull, C. L. H., Yang, H., Li, Z.-Y., et al. 2018, ApJ, 860, 82, doi: 10.3847/1538-4357/aabfeb
  * Kataoka et al. (2016a) Kataoka, A., Muto, T., Momose, M., Tsukagoshi, T., & Dullemond, C. P. 2016a, ApJ, 820, 54, doi: 10.3847/0004-637X/820/1/54
  * Kataoka et al. (2019) Kataoka, A., Okuzumi, S., & Tazaki, R. 2019, ApJ, 874, L6, doi: 10.3847/2041-8213/ab0c9a
  * Kataoka et al. (2017) Kataoka, A., Tsukagoshi, T., Pohl, A., et al. 2017, ApJ, 844, L5, doi: 10.3847/2041-8213/aa7e33
  * Kataoka et al. (2015) Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78, doi: 10.1088/0004-637X/809/1/78
  * Kataoka et al. (2016b) Kataoka, A., Tsukagoshi, T., Momose, M., et al. 2016b, ApJ, 831, L12, doi: 10.3847/2041-8205/831/2/L12
  * Kirchschlager et al. (2019) Kirchschlager, F., Bertrang, G. H. M., & Flock, M. 2019, MNRAS, 488, 1211, doi: 10.1093/mnras/stz1763
  * Lazarian & Hoang (2007) Lazarian, A., & Hoang, T. 2007, ApJ, 669, L77, doi: 10.1086/523849
  * Lee & Draine (1985) Lee, H. M., & Draine, B. T. 1985, ApJ, 290, 211, doi: 10.1086/162974
  * Liu (2019) Liu, H. B. 2019, ApJ, 877, L22, doi: 10.3847/2041-8213/ab1f8e
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591
  * Miyake & Nakagawa (1993) Miyake, K., & Nakagawa, Y. 1993, Icarus, 106, 20, doi: 10.1006/icar.1993.1156
  * Mori et al. (2019) Mori, T., Kataoka, A., Ohashi, S., et al. 2019, ApJ, 883, 16, doi: 10.3847/1538-4357/ab3575
  * Ohashi et al. (2018) Ohashi, S., Kataoka, A., Nagai, H., et al. 2018, ApJ, 864, 81, doi: 10.3847/1538-4357/aad632
  * Okuzumi & Tazaki (2019) Okuzumi, S., & Tazaki, R. 2019, ApJ, 878, 132, doi: 10.3847/1538-4357/ab204d
  * Rebull et al. (2004) Rebull, L. M., Wolff, S. C., & Strom, S. E. 2004, AJ, 127, 1029, doi: 10.1086/380931
  * Ricci et al. (2010) Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15, doi: 10.1051/0004-6361/200913403
  * Robitaille et al. (2007) Robitaille, T. P., Whitney, B. A., Indebetouw, R., & Wood, K. 2007, ApJS, 169, 328, doi: 10.1086/512039
  * Stephens et al. (2014) Stephens, I. W., Looney, L. W., Kwon, W., et al. 2014, Nature, 514, 597, doi: 10.1038/nature13850
  * Stephens et al. (2017) Stephens, I. W., Yang, H., Li, Z.-Y., et al. 2017, ApJ, 851, 55, doi: 10.3847/1538-4357/aa998b
  * Tamura et al. (1999) Tamura, M., Hough, J. H., Greaves, J. S., et al. 1999, ApJ, 525, 832, doi: 10.1086/307953
  * Tamura et al. (1995) Tamura, M., Hough, J. H., & Hayashi, S. S. 1995, ApJ, 448, 346, doi: 10.1086/175965
  * Tazaki et al. (2017) Tazaki, R., Lazarian, A., & Nomura, H. 2017, ApJ, 839, 56, doi: 10.3847/1538-4357/839/1/56
  * Testi et al. (2014) Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 339, doi: 10.2458/azu_uapress_9780816531240-ch015
  * Warren & Brandt (2008) Warren, S. G., & Brandt, R. E. 2008, Journal of Geophysical Research (Atmospheres), 113, D14220, doi: 10.1029/2007JD009744
  * Yang et al. (2016) Yang, H., Li, Z.-Y., Looney, L., & Stephens, I. 2016, MNRAS, 456, 2794, doi: 10.1093/mnras/stv2633
  * Yang et al. (2019) Yang, H., Li, Z.-Y., Stephens, I. W., Kataoka, A., & Looney, L. 2019, MNRAS, 483, 2371, doi: 10.1093/mnras/sty3263
  * Zhu et al. (2019) Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJ, 877, L18, doi: 10.3847/2041-8213/ab1f8c


# Comparisons of Graph Neural Networks on Cancer Classification Leveraging a
Joint of Phenotypic and Genetic Features

David Oniani
Mayo Clinic
Rochester, MN 55902
<EMAIL_ADDRESS>
&Chen Wang
Mayo Clinic
Rochester, MN 55902
<EMAIL_ADDRESS>
&Yiqing Zhao
Mayo Clinic
Rochester, MN 55902
<EMAIL_ADDRESS>
&Andrew Wen
Mayo Clinic
Rochester, MN 55902
<EMAIL_ADDRESS>
&Hongfang Liu
Mayo Clinic
Rochester, MN 55902
<EMAIL_ADDRESS>
&Feichen Shen
Mayo Clinic
Rochester, MN 55902
<EMAIL_ADDRESS>
###### Abstract

Cancer is responsible for millions of deaths worldwide every year. Although
significant progress has been achieved in cancer medicine, many issues remain
to be addressed for improving cancer therapy. Appropriate cancer patient
stratification is the prerequisite for selecting appropriate treatment plan,
as cancer patients are of known heterogeneous genetic make-ups and phenotypic
differences. In this study, built upon deep phenotypic characterizations
extractable from Mayo Clinic electronic health records (EHRs) and genetic test
reports for a collection of cancer patients, we evaluated various graph neural
networks (GNNs) leveraging a joint of phenotypic and genetic features for
cancer type classification. Models were applied and fine-tuned on the Mayo
Clinic cancer disease dataset. The assessment was done through the reported
accuracy, precision, recall, and F1 values as well as through F1 scores based
on the disease class. Per our evaluation results, GNNs on average outperformed
the baseline models with mean statistics always being higher that those of the
baseline models (0.849 vs 0.772 for accuracy, 0.858 vs 0.794 for precision,
0.843 vs 0.759 for recall, and 0.843 vs 0.855 for F1 score). Among GNNs,
ChebNet, GraphSAGE, and TAGCN showed the best performance, while GAT showed
the worst. We applied and compared eight GNN models including AGNN, ChebNet,
GAT, GCN, GIN, GraphSAGE, SGC, and TAGCN on the Mayo Clinic cancer disease
dataset and assessed their performance as well as compared them with each
other and with more conventional machine learning models such as decision
tree, gradient boosting, multi-layer perceptron, naive bayes, and random
forest which we used as the baselines.

_K_ eywords Cancer Classification $\cdot$ Phenotypic and Genetic
Characterization $\cdot$ Graph Neural Networks

##  1 Introduction

Cancer is a type of disease in which abnormal cells divide without control and
can invade nearby tissues and other parts of the body via the blood and lymph
systems [1, 2]. Nowadays, cancer is responsible for millions of deaths
worldwide every year. Although significant progress has been achieved in
cancer medicine, many issues remain to be addressed for improving cancer
therapy. The goal of Precision Oncology is to enable oncologist practitioners
to make better clinical decisions by incorporating individual cancer patients’
genetic information and clinical characteristics [3]. Precision Oncology is
expected to improve selection of targeted therapies, tailored to individual
patients and ultimately improve cancer patients’ outcomes. However, patient
heterogeneity needs to be considered while providing cancer treatments. It is
important to take patient subgroups into consideration in which the mechanism
of diseases with the same group are more likely to be homogeneous [4].
Therefore, appropriate cancer patient subgrouping is the prerequisite for
selecting appropriate treatment plan, as cancer patients are of known
heterogeneous genetic make-ups and phenotypic differences.

Advances in cancer care have demonstrated the benefits of Individualized
Medicine research, and the premise of Precision Oncology is to refer similar
cancer outcomes to similar patients with similar features [5]. Phenotypic
characterizations are one of the key features to group patients by analyzing
their clinical observations [6]. Taking breast cancer as an example, specific
sub-phenotypes have been leveraged to make decision making for therapies
selection in practice [7]. Moreover, phenotypic heterogeneity amongs 7
subtypes of triple-negative breast cancer could be observed and analyzed [8].
With the advance of machine learning techniques, computational deep
phenotyping is widely used and described as the process to analyze phenotypic
abnormalities precisely and comprehensively, in order to receive a better
understanding of the natural history of a disease and its genotype-phenotype
associations. For example, the DeepPhe information model is able to extract
cancer-related phenotypes as categorical features by mining clinical
narratives or structured datasets [9]. Beeghly-Fadiel et al. have extended
DeepPhe to ovarian cancer [10]. Besides categorical phenotypic features,
embeddings representations of phenotypes are also significant to deep
phenotyping. HPO2Vec+ provides a way to learn node embeddings for the Human
Phenotype Ontology (HPO) and be able to represent patients with vectorized
representations for individual-level phenotypic similarity measurement [11,
12]. Maldonado et al. presented an approach that uses Generative Adversarial
Networks (GANs) [13] for learning UMLS embeddings and demonstrated the use
case on clinical prediction tasks [14].

In addition to phenotypic characterizations, features contained in genetic
test reports provide crucial evidence for the decision making in practice.
Genetic tests for targeted cancer therapy detect mutations in the DNA of
cancer cells. Being aware of whether the cancer has a particular mutation can
help guide the type of treatment. [15]. Specifically, pathogenic variants
demonstrate the genetic alteration that increases an individual’s
susceptibility or predisposition to a certain cancer type, and variant of
uncertain significance (VUS) indicates a genetic variant that has been
identified but without knowing too much details regarding significance to the
function of an organism. There exist some studies that use machine learning to
support variant interpretation in a clinical setting. For example, the LEAP
model classifies variants based on features collected from functional
predictions, splice predictions, clinical observation data and so on [16].
Some guidelines to define variant pathogenicity for somatic and VUS
classification were established [17, 18, 19]. Mut2Vec provides a way to
represent cancerous mutations with distributed embeddings, in order to support
cancer classification and drug sensitivity prediction [20].

Graph Neural Networks (GNNs) [21] represent a series of deep learning methods
that applied on graph domain. GNNs has been largely adopted for graph data
characterization and analysis due to the promising performance. In this study,
we constructed a graph leveraging phenotypic characterization and genetic
features collected from patients and then compared performance on cancer type
classification using various exisiting GNNs. Specifically, we compared
performance amongst 8 GNNs, including Attentive Graph Neural Network (AGNN),
Chebyshev Graph Neural Network (ChebNet), Graph Convolutional Network (GCN),
Graph Attention Network (GAT), Graph Isomorphism Network (GIN), Graph SAmple
and aggreGatE (GraphSAGE), Simplifying Graph Convolutional Network (SGC), and
Topology Adaptive Graph Convolutional Network (TAGCN).

##  2 Materials

We collected a cancer patient cohort (n=794) with accessible Mayo Clinic EHRs
and foundation genetic test reports. Seven types of cancer are included in
this cohort: lung cancer (n=223), prostate cancer (n=66), breast cancer
(n=53), ovarian cancer (n=91), pancreas cancer (n=104), colon/rectum cancer
(n=140), and liver cancer (n=107).

###  2.1 Extract phenotypic features from the Mayo Clinic EHRs

Given patient IDs retrieved from the cohort, we first applied an existing Mayo
Clinic NLP-as-a-service implementation [22] to extract all the diagnostic-
related terms from clinical notes and radiology report, normalized by the
Unified Medical Language System (UMLS) [23]. We then used a previously
developed Human Phenotype Ontology (HPO)-based annotation pipeline [24] to
further map the UMLS terms to more specific HPO normalized phenotypic terms
[25]. Specifically, for each cancer patient, we collected their phentoypic
information within the last twelve months of the first mention of cancer
diagnoses.

###  2.2 Collect genetic feature from the FoundationOne CDx genetic test
reports

FoundationOne CDx is a next-generation sequencing (NGS) based assay that
identifies genomic findings within hundreds of cancer-related genes. The
foundationOne genetic test report is designed to provide physicians with
clinically actionable information. The report for each patient is represented
in the PDF format. We first used a python library named pdfminer to parse all
the PDF report into text format, for the convenience of data processing. We
then parsed out genes alteration and VUS sections for each patient. In
particular, we treated mutations mined from the genes alteration section as
pathogenic mutation, and we named mutations mined from the VUS section as VUS
mutation. We used multi-hot encoding to represent genetic features for each
patient, where 1 indicates the existence of pathogenic mutation or VUS, and 0
otherwise.

##  3 Methods

We applied eight different machine learning models and assessed their
performance on Mayo Clinic cancer disease dataset. The main goal was to
conduct cancer disease classification tasks and evaluate the performance of
the models as well as to compare the models with each other.

The models included:

  * •

Attentive Graph Neural Network (AGNN) [26]

  * •

Chebyshev Graph Neural Network (ChebNet) [27]

  * •

Graph Convolutional Network (GCN) [28]

  * •

Graph Attention Network (GAT) [29]

  * •

Graph Isomorphism Network (GIN) [30]

  * •

Graph SAmple and aggreGatE (GraphSAGE) [31]

  * •

Simplifying Graph Convolutional Network (SGC) [32]

  * •

Topology Adaptive Graph Convolutional Network (TAGCN) [33]

We also compared them to the baseline models such as decision tree, gradient
boosting, multi-layer perceptron, naive bayes, and random forest. All of the
models were fine-tuned. We present the results in the results section.

###  3.1 Data Preprocessing and Graph Generation

For phenotypes, we first removed some general phenotypic terms mined from
EHRs, including cyst, pain, carcinoma, neoplasm, symptoms, disease, minor
(disease), and sarcoma - cateogry (morphologic abnormality). We then filtered
out phenotypic terms with low frequency across all the patients, in order to
largely avoid noises. By heuristics, we set the frequency threshold for
selecting phenotypes as 20. Similarly, for genetic features, we only kept the
ones that appear more than 10 times across all the patients.

Associations between the patients and the features were captured to build the
feature graph. Specifically, we made a link between phenotypic feature and
genetic feature if both of them are associated with the same patient. As a
result, the feature graph can be defined as $G=(N,L)$, where $N$ represents
the set of phenotypic/genetic nodes and $L$ represents the set of links/edges
between these nodes. We viewed graph $G$ as an undirected graph that bridges
mutations and phenotypes based on information collected from heterogeneous
cancer patients in clinical practice.

###  3.2 GNN Models

####  3.2.1 AGNN

AGNN is a GNN model where the relations between the nodes are denoted by the
differentiable attention mechanism. The whole model is end-to-end trainable,
as all the functions in AGNN are parameterized by neural networks. Given any
two phenotypic/genetic feature representations $f_{i}$, $f_{j}$, and their
edge representation ${e}_{i,j}$, the message aggregation is done as follows:

$\textbf{m}_{i}^{k}=\sum_{f_{j}\in\mathcal{N}_{i}}\textbf{m}_{j,i}^{k}=\sum_{f_{j}\in\mathcal{N}_{i}}M(\textbf{h}_{j}^{k-1},\textbf{e}_{i,j}^{j-1}),$
(1)

where $M()$ is a message function, $\textbf{m}_{i}^{k}$ is a received message,
$\mathcal{N}_{i}$ represents all the neighbor nodes for $f_{i}$, and $k$
indicates the number of steps. Given a state update function $T()$, the node
representation is updated as follows:

$\textbf{h}_{i}^{k}=T(\textbf{h}_{i}^{k-1},\textbf{m}_{i}^{k}),$ (2)

where $h_{i}^{0}=f_{i}$. After $k$ iterations of aggregation,
$\textbf{h}_{i}^{k}$ captures the relations within the $k$-hop neighborhood of
node $f_{i}$.

####  3.2.2 ChebNet

ChebNet is a spectral-based graph learning algorithm. It is also a highly
efficient version of GCNN [34]. Its convolutions are defined in a similar
fashion as they are usually defined on images, where a pixel value is
multiplied by the value of the kernel. In ChebNet, spectral convolutions are
defined as the multiplication of a signal by a kernel. Kernel is comprised of
Chebyshev polynomials. The polynomials are derived from the diagonal matrix of
Laplacian eigenvalues. The kernel can be represented as the truncated
expansion of order $K-1$ as follows:

$g_{\theta}(\Lambda)=\sum_{k=0}^{K-1}\theta_{k}T_{k}(\tilde{\Lambda}),$ (3)

where $g_{\theta}$ is a kernel applied to $\Lambda$, the diagonal matrix of
Laplacian eigenvalues, and $\theta\in\mathbb{R}^{K}$ is the vector of
Chebyshev coefficients. $\tilde{\Lambda}$ represents the diagonal matrix of
the scaled Laplacian eigenvalues. $K$ and $k$ represent the largest and the
smallest order neighborhoods respectively.
$T_{k}(\tilde{\Lambda})\in\mathbb{R}^{n}\times\mathbb{R}^{n}$ is the Chebyshev
polynomial of order k. Finally, the kernel is the sum of Chebyshev polynomials
from $k$ to $K-1$.

Hence, the ChebNet convolution is:

$g_{\theta}(\Lambda)\circ
x=\sum_{k=0}^{K-1}\theta_{k}T_{k}(\tilde{\Lambda})x,$ (4)

where x is the signal. In this equation, $\theta$’s are the parameters to be
trained. It should also be noted that the representation of $\tilde{\Lambda}$
maps the eigenvalues in the interval $[-1,1]$. Such representation of the
convolution greatly decreases the number of computations to be performed.

####  3.2.3 GCN

GCN can be viewed as a simplified version of ChebNet. Its convolution is
defined as follows:

$g_{\theta^{\prime}}(\Lambda)\star
x=\sum_{k=0}^{K-1}\theta_{k}^{\prime}T_{k}(\tilde{L})x,$ (5)

where $g_{0}$ is a kernel which is multiplied by the graph signal $x$.
$\theta_{k}^{\prime}$ is a vector of Chebyshev coefficients. $K$ is the $K$th
order neighborhood and $k$ is the closest order neighborhood. $T$ is the
Chebyshev polynomial applied to $\tilde{L}$.

ChebNet and GCN are similar. The biggest difference is in the choice for the
value of K. In the GCN, the layer wise convolution is limited to K = 1. This
helps avoid overfitting on the graph’s local neighborhood.

####  3.2.4 GAT

The Graph Attention Network or GAT is a non-spectral learning method which
utilizes the spatial information of the node directly for learning. This is in
contrast to the spectral approach of the GCN which mirrors the same basics as
the CNN. On the other hands, GAT, similar to GCN, uses node self-attention
features and neighbor features for training the model. It also uses multi-head
attention (similar to models such as BERT [35]).

The first step performed by the Graph Attention Layer is to apply a linear
transformation - Weighted matrix W to the feature vectors of the nodes.

GAT uses attention coefficients for getting the relative importance of
neighboring features with respect to each other. For neighboring nodes $f_{i}$
and $f_{j}$, the attention coefficient is calculated as follows:

$e_{i,j}=a(\textbf{W}\overrightarrow{h_{i}},\textbf{W}\overrightarrow{h_{j}}),$
(6)

where
$a=\mathbb{R}^{F^{\prime}}\times\mathbb{R}^{F^{\prime}}\rightarrow\mathbb{R}$
is shared attention mechanism, $h_{i},h_{j}$ are input features, and W is a
weighted matrix.

For making the attention coefficients easily comparable across different
feature nodes, GAT normalizes them across all choices of $f_{j}$ by applying
the softmax function. Finally, the aggregation is done as follows:

$\overrightarrow{h_{i}^{\prime}}=\sigma\bigg{(}\frac{1}{K}\sum_{k=1}^{K}\sum_{f_{j}\in\mathcal{N}_{i}}\alpha_{ij}^{k}\textbf{W}^{k}\overrightarrow{h}_{j})\bigg{)}$
(7)

####  3.2.5 GIN

GIN is another graph neural network that updates the phenotypic/genetic node
representations using the multi-layer perceptron according to the following
equation:

$h_{i}^{(k)}=\textnormal{MLP}^{k}\bigg{(}\big{(}1+\epsilon^{(k)}\big{)}\cdot
h_{i}^{k+1}+\sum_{u\in\mathcal{N}(i)}h_{j}^{(k-1)}\bigg{)},$ (8)

where $h_{i}^{(k)}$ is the feature vector of node $f_{i}$ at the $k$-th
iteration/layer, $\mathcal{N}(i)$ is the set of nodes adjacent to $f_{i}$, and
MLP is a multi-layer perceptron. Aggregation is done as follows:

$h_{G}=\textnormal{CONCAT}\bigg{(}\textnormal{READOUT}\bigg{(}\bigg{\\{}h_{i}^{(k)}i\in
G\bigg{\\}}\bigg{)}\mid k=0,1,\dots,K),$ (9)

where CONCAT aggregates the embeddings of feature nodes. In this study, we
adopted mean aggregation operation for GIN to learn phenotypic and genetic
feature representations .

####  3.2.6 GraphSAGE

GraphSAGE is a framework for inductive representation learning. It is suited
for generating low-dimensional embeddings for nodes and for graphs that grow
in size over time. Since GraphSAGE is inductive by nature, generating vector
representation for new nodes is a lot faster compared to transductive
techniques such as DeepWalk [36]. In this study, each phenotypic/genetic
feature node is represented by the aggregation of its neighborhood.
Specifically, we used GraphSAGE with the mean aggregator. Given any feature
node $f_{i}$, the equation for the mean aggregation is defined as:

$h_{k}^{v}\leftarrow\sigma(W\cdot\textnormal{MEAN}(\\{h_{i}^{k-1}\\}\cup\\{h_{i}^{k-1},\forall
i\in\mathcal{N}(i)\\})$ (10)

####  3.2.7 SGC

SGC is a simplified version of GCN where nonlinearities are successively
removed and weight matrices between consecute layers are collapsed. It is
comprised of a graph-based preprocessing step followed by a standard multi-
class logistic regression, which is suitable for our task on multi-class
cancer classification.

####  3.2.8 TAGCN

TAGCN is a graph neural network where each neuron in the graph convolutional
layer is connected only to a local region (local phenotypic/genetic feature
nodes and links between them) in the vertex domain of the input data volume,
which is adaptive to the graph topology. It uses a set of fixed-size learnable
filters to perform convolutions on graphs. In this study, the topologies of
these filters are adaptive to the topology of the entire feature graph when
they scan the graph to perform convolution.

##  4 Evaluation

We evaluated the results using metrics including accuracy, precision, recall,
and F1 score. For further performance assessment, ROC curve was also generated
for each model. Equations below describe these metrics:

$\displaystyle Accuracy=\frac{TP+TN}{TP+TN+FP+FN}$ (11) $\displaystyle
Precision=\frac{TP}{TP+FN}$ (12) $\displaystyle Recall=\frac{TP}{TP+FN}$ (13)
$\displaystyle F1=\frac{2\times Precision\times Recall}{Precision+Recall}$
(14)

##  5 Results

Table 1 shows the metrics for the GNN models. Additionally, it shows the
metrics for the baseline models including decision tree, gradient boosting,
multi-layer perceptron, naive bayes, and random forest. Our goal was to also
compare the results with those of more traditional machine learning
techniques.

Table 1: Evaluation results for the eight GNN models and five baseline models. Model | Accuracy | Precision | Recall | F1 Score
---|---|---|---|---
AGNN | 0.877 | 0.890 | 0.889 | 0.885
ChebNet | 0.901 | 0.902 | 0.920 | 0.907
GAT | 0.753 | 0.726 | 0.684 | 0.697
GCN | 0.802 | 0.855 | 0.787 | 0.804
GIN | 0.827 | 0.821 | 0.809 | 0.805
GraphSAGE | 0.901 | 0.926 | 0.910 | 0.915
SGC | 0.827 | 0.836 | 0.832 | 0.824
TAGCN | 0.901 | 0.907 | 0.910 | 0.904
Decision Tree | 0.700 | 0.737 | 0.735 | 0.722
Gradient Boosting | 0.877 | 0.900 | 0.868 | 0.871
Multi-layer Perceptron | 0.864 | 0.902 | 0.851 | 0.870
Naive Bayes | 0.531 | 0.533 | 0.444 | 0.425
Random Forest | 0.888 | 0.897 | 0.895 | 0.888

ChebNet, GraphSage, and TAGCN showed the best performance accuracy-wise, all
having roughly the same accuracy value of 0.901. AGNN came next with the
accuracy value of 0.877. GIN and SGC both had the approximate accuracy of
0.827. GCN came next with the accuracy value of 0.827. The worst accuracy was
demonstrated by GAT, with the score of 0.753. As compared to accuracy values
of the baseline models, three GNN models outperformed random forest
classifier, which had the highest accuracy. Additionally, the mean accuracy
score for GNNs was approximately 0.849 while the mean accuracy value for the
baseline models was roughly 0.772.

GraphSAGE had the highest precision value of 0.926. TAGCN and ChebNet came
next with the values of 0.907 and 902 respectively. AGNN and GCN had the
precision values of 0.890 and 0.855 respectively. SGC had the precision value
of 0.836 and GIN’s precision was 0.821. GAT had the precision value 0.726,
which was the lowest among the GNNs. Across all baseline models, multi-layer
perceptron had the highest precision value of 0.902 and gradient boosting came
close with the value of 0.900. GraphSAGE and TAGCN outperformed the best
baseline and ChebNet had the same approximate value. The mean precision for
GNNs was 0.858 and the mean precision for the baseline models was 0.794.

As for the recall, ChebNet showed the best performance among the GNNs, with
the value of 0.920. GraphSAGE and TAGCN both had the same recall value of
0.910. AGNN came next with the value of 0.889. SGC had the value of 0.832. GIN
and GCN had the recall values of 0.809 and 787 respectively. GAT had the value
of 0.684. The highest value for the baseline models was 0.895, with the random
forest classifier. ChebNet, GraphSAGE, and TAGCN had the higher recall values.
The mean recall values for GNN and baseline models were 0.843 and 0.759
respectively.

Across all GNNs, GraphSAGE had the highest F1 score, roughly equal to 0.915.
ChebNet and TAGCN came next with the values of 0.907 and 0.904 respectively.
AGNN had the F1 score of 0.855. SGC and GIN had the scores of 0.824 and 0.805
respectively. GCN had the score of 0.804. GAT had the worst F1 score of 0.697.
Random forest classifier had the highest F1 score among all baseline models,
with the value of 0.888. Among GNNs, GraphSAGE, ChebNet, and TAGCN had higher
scores. The mean F1 score for the GNN models was 0.843 and the mean F1 score
for the baseline models was 0.755.

Overall, GNN models, on average, showed better performance in all metrics than
the baseline models. The mean values for all four metrics (accuracy,
precision, recall and F1 score) were higher for GNN models. Among the GNN
models, ChebNet, GraphSage, and TAGCN had the best overall performance and GAT
showed the worst overall performance.

Table 2: Disease-specific evaluation of models based on F1 scores. Model | Lung | Prostate | Breast | Ovarian | Pancreas | Colon/Rectum | Liver
---|---|---|---|---|---|---|---
AGNN | 0.89 | 0.95 | 0.83 | 0.87 | 0.93 | 0.88 | 0.85
ChebNet | 0.89 | 0.90 | 0.93 | 0.86 | 1.00 | 0.92 | 0.85
GAT | 0.29 | 0.70 | 0.80 | 0.80 | 0.77 | 0.67 | 0.86
GCN | 0.67 | 0.90 | 0.77 | 0.81 | 1.00 | 0.62 | 0.86
GIN | 0.67 | 0.76 | 0.85 | 0.86 | 0.75 | 0.92 | 0.83
GraphSAGE | 1.00 | 0.90 | 0.93 | 0.87 | 0.93 | 0.92 | 0.85
SGC | 0.67 | 0.90 | 0.86 | 0.81 | 0.93 | 0.71 | 0.89
TAGCN | 0.89 | 0.90 | 0.93 | 0.87 | 0.93 | 0.96 | 0.85
Decision Tree | 0.80 | 0.67 | 0.71 | 0.63 | 0.86 | 0.67 | 0.75
Gradient Boosting | 0.89 | 0.86 | 0.89 | 0.86 | 0.77 | 0.91 | 0.92
Multi-layer Perceptron | 0.86 | 0.95 | 0.81 | 0.88 | 0.93 | 0.81 | 0.85
Naive Bayes | 0.00 | 0.35 | 0.65 | 0.67 | 0.62 | 0.44 | 0.25
Random Forest | 0.89 | 0.90 | 0.88 | 0.84 | 0.86 | 0.83 | 0.89

In the lung cancer classification task, GraphSAGE showed F1 score of roughly
1.00 and clearly outperformed any GNN or baseline model in this specific
disease classification task. AGNN, ChebNet, and TAGCN came next, all three
having the same the F1 score of 0.89. GCN, GIN, and SGC all had the F1 score
of 0.67. GAT showed the worst F1 score across the GNNs, approximately equal to
0.29. Among the baseline models, gradient boosting and random forest showed
the best F1 scores, both having 0.89. They showed performance similar to AGNN,
ChebNet, and TAGCN, yet could not outperform GraphSAGE.

As for the prostate cancer, AGNN had the highest score across the GNNS, with
the value of 0.95. Multi-layer perceptron baseline had the same F1 score.
ChebNet, GCN, SGC, GraphSAGE, SGC, and TAGCN all had the same F1 scores of
0.90 which was marginally worse than that of multi-layer perceptron baseline.
GIN had the lowest score among the GNNs, but performed well as compared with
baseline models. GAT had the worst performance among GNNs and second worst
overall performance, with the F1 score of 0.29.

In the breast cancer classification task, ChebNet, GraphSAGE, and TAGCN had
the same F1 score of 0.93. SGC and GIN had the F1 scores of 0.86 and 0.85
respectively. AGNN had the F1 score of 0.83 and GAT the F1 score of 0.80.
Among GNNs, GCN had the lowest F1 score, approximately equal to 0.77. The
highest score among baseline models was 0.89 with gradient boosting. ChebNet,
GraphSAGE, and TAGCN all had the higher F1 scores, but the rest of the models
showed a worse performance.

For the ovarian cancer, AGNN, GraphSAGE, and TAGCN all had the same F1 score,
equal to 0.87. ChebNet and GIN had marginally lower F1 scores, both equal to
0.86. GCN and SGC had the F1 scores of 0.81. GAT had the F1 score of 0.80.
Multi-layer perceptron was the baseline with the highest F1 score of 0.88. It
outperformed all GNN and baseline models.

In pancreas cancer classification, ChebNet and GCN had nearly perfect F1
scores, both approximately equal to 1.0. AGNN, GraphSAGE, SGC, and TAGCN came
next with the same score of 0.93. GAT and GIN had the F1 scores of 0.77 and
0.75 respectively. Among the baseline models, multi-layer perceptron showed
the highest score of 0.93, which was lower than all ChebNet and GCN, but on
par with AGNN, GraphSAGE, SGC, and TAGCN and higher that both GAT and GCN.

As for colon/rectum cancer classification, TAGCN had the highest F1 score
among the GNNs, equal to 0.96. ChebNet, GIN, and GraphSAGE had the same F1
score of 0.92. AGNN had the F1 score of 0.88. SGC, GCN, and GAT had the scores
of 0.71, 0.67, and 0.62 respectively. Gradient boosting had the highest F1
score across all baselines, with the value of 0.91.

In the liver cancer classification task, SGC had the highest F1 score, equal
to 0.89. GAT and GCN came next with the same F1 score of 0.86. AGNN, ChebNet,
GraphSAGE, and TAGCN all had marginally lower F1 score of 0.85. GIN had the
score of 0.83. Interestingly, gradient boosting outperformed all GNN as well
as baseline models with the F1 score of 0.92.

To recapitulate, out of the seven classes, in four (lung cancer, breast
cancer, pancreas cancer, and colon/rectum cancer), the highest F1 scores came
from GNNs. In one of the classes (prostate cancer), the highest overall score
was present among both GNNs and baseline models had the same F1 score. In
three of the classes (ovarian cancer and liver cancer), highest F1 scores came
from baseline models.

##  6 Conclusions

We harvested a cancer dataset comprising of 7 cancer disease classes and
applied 8 GNNs - AGNN, ChebNet, GAT, GCN, GIN, GraphSAGE, SGC, and TAGCN - for
disease classification purposes. We assessed the performance of the proposed
models using various metrics as well as compared them to the baseline models
including decision tree, gradient boosting, multi-layer perceptron, naive
bayes, and random forest. In this study, we discovered the potential usage of
applying deep learning algorithms over joint phenotypic and genetic features
for a better cancer type characterization. Moreover, we observed that ChebNet,
GraphSage, and TAGCN outperformed other GNNs and baseline models, showing
strong potential to be used for cancer type classification tasks.

## Acknowledgements

This work has been supported by the Gerstner Family Foundation, Center for
Individualized Medicine (CIM) of Mayo Clinic.

## References

  * [1] Yuchen Yuan, Yi Shi, Changyang Li, Jinman Kim, Weidong Cai, Zeguang Han, and David Dagan Feng.  Deepgene: an advanced cancer type classifier based on deep learning and somatic point mutations.  BMC bioinformatics, 17(17):476, 2016.
  * [2] Michael Feuerstein.  Defining cancer survivorship.  Journal of cancer survivorship: research and practice, 1(1):5–7, 2007.
  * [3] Lee Schwartzberg, Edward S Kim, David Liu, and Deborah Schrag.  Precision oncology: who, how, what, when, and when not?  American Society of Clinical Oncology Educational Book, 37:160–169, 2017.
  * [4] Forum on Neuroscience, Nervous System Disorders, et al.  Improving and accelerating therapeutic development for nervous system disorders: Workshop summary.  2014\.
  * [5] Enea Parimbelli, Simone Marini, Lucia Sacchi, and Riccardo Bellazzi.  Patient similarity for precision medicine: A systematic review.  Journal of biomedical informatics, 83:87–96, 2018.
  * [6] Nathan D Seligson, Jeremy L Warner, William S Dalton, David Martin, Robert S Miller, Debra Patt, Kenneth L Kehl, Matvey B Palchuk, Gil Alterovitz, Laura K Wiley, et al.  Recommendations for patient similarity classes: results of the amia 2019 workshop on defining patient similarity.  Journal of the American Medical Informatics Association, 2020.
  * [7] Brian D Lehmann, Joshua A Bauer, Xi Chen, Melinda E Sanders, A Bapsi Chakravarthy, Yu Shyr, Jennifer A Pietenpol, et al.  Identification of human triple-negative breast cancer subtypes and preclinical models for selection of targeted therapies.  The Journal of clinical investigation, 121(7):2750–2767, 2011.
  * [8] Hiroko Masuda, Keith A Baggerly, Ying Wang, Ya Zhang, Ana Maria Gonzalez-Angulo, Funda Meric-Bernstam, Vicente Valero, Brian D Lehmann, Jennifer A Pietenpol, Gabriel N Hortobagyi, et al.  Differential response to neoadjuvant chemotherapy among 7 triple-negative breast cancer molecular subtypes.  Clinical cancer research, 19(19):5533–5540, 2013.
  * [9] Guergana K Savova, Eugene Tseytlin, Sean Finan, Melissa Castine, Timothy Miller, Olga Medvedeva, David Harris, Harry Hochheiser, Chen Lin, Girish Chavan, et al.  Deepphe: a natural language processing system for extracting cancer phenotypes from clinical records.  Cancer research, 77(21):e115–e118, 2017.
  * [10] Alicia Beeghly-Fadiel, Jeremy L Warner, Sean Finan, James Masanz, Harry Hochheiser, and Guergana Savova.  Deep phenotype extraction to facilitate cancer research: Extending deepphe to ovarian cancer, 2019.
  * [11] Feichen Shen, Suyuan Peng, Yadan Fan, Andrew Wen, Sijia Liu, Yanshan Wang, Liwei Wang, and Hongfang Liu.  Hpo2vec+: Leveraging heterogeneous knowledge resources to enrich node embeddings for the human phenotype ontology.  Journal of biomedical informatics, 96:103246, 2019.
  * [12] Feichen Shen, Sijia Liu, Yanshan Wang, Liwei Wang, Andrew Wen, Andrew H Limper, and Hongfang Liu.  Constructing node embeddings for human phenotype ontology to assist phenotypic similarity measurement.  In 2018 IEEE International Conference on Healthcare Informatics Workshop (ICHI-W), pages 29–33. IEEE, 2018.
  * [13] Ian J Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial networks.  arXiv preprint arXiv:1406.2661, 2017.
  * [14] Ramon Maldonado, Meliha Yetisgen, and Sanda M Harabagiu.  Adversarial learning of knowledge embeddings for the unified medical language system.  AMIA Summits on Translational Science Proceedings, 2019:543, 2019\.
  * [15] Charles Sawyers.  Targeted cancer therapy.  Nature, 432(7015):294–297, 2004.
  * [16] Carmen Lai, Anjali D Zimmer, Robert O’Connor, Serra Kim, Ray Chan, Jeroen van den Akker, Alicia Y Zhou, Scott Topper, and Gilad Mishne.  Leap: Using machine learning to support variant classification in a clinical setting.  Human Mutation, 41(6):1079–1090, 2020.
  * [17] Giulia Federici and Silvia Soddu.  Variants of uncertain significance in the era of high-throughput genome sequencing: a lesson from breast and ovary cancers.  Journal of Experimental & Clinical Cancer Research, 39(1):1–12, 2020.
  * [18] Sue Richards, Nazneen Aziz, Sherri Bale, David Bick, Soma Das, Julie Gastier-Foster, Wayne W Grody, Madhuri Hegde, Elaine Lyon, Elaine Spector, et al.  Standards and guidelines for the interpretation of sequence variants: a joint consensus recommendation of the american college of medical genetics and genomics and the association for molecular pathology.  Genetics in medicine, 17(5):405–423, 2015.
  * [19] Heidi L Rehm, Sherri J Bale, Pinar Bayrak-Toydemir, Jonathan S Berg, Kerry K Brown, Joshua L Deignan, Michael J Friez, Birgit H Funke, Madhuri R Hegde, and Elaine Lyon.  Acmg clinical laboratory standards for next-generation sequencing.  Genetics in medicine, 15(9):733–747, 2013.
  * [20] Sunkyu Kim, Heewon Lee, Keonwoo Kim, and Jaewoo Kang.  Mut2vec: distributed representation of cancerous mutations.  BMC medical genomics, 11(2):33, 2018.
  * [21] Jie Zhou, Ganqu Cui, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changcheng Li, and Maosong Sun.  Graph neural networks: A review of methods and applications.  arXiv preprint arXiv:1812.08434, 2018.
  * [22] Andrew Wen, Sunyang Fu, Sungrim Moon, Mohamed El Wazir, Andrew Rosenbaum, Vinod C Kaggal, Sijia Liu, Sunghwan Sohn, Hongfang Liu, and Jungwei Fan.  Desiderata for delivering nlp to accelerate healthcare ai advancement and a mayo clinic nlp-as-a-service implementation.  NPJ digital medicine, 2(1):1–7, 2019.
  * [23] Olivier Bodenreider.  The unified medical language system (umls): integrating biomedical terminology.  Nucleic acids research, 32(suppl_1):D267–D270, 2004.
  * [24] Feichen Shen, Liwei Wang, and Hongfang Liu.  Phenotypic analysis of clinical narratives using human phenotype ontology.  Studies in health technology and informatics, 245:581, 2017.
  * [25] Peter N Robinson, Sebastian Köhler, Sebastian Bauer, Dominik Seelow, Denise Horn, and Stefan Mundlos.  The human phenotype ontology: a tool for annotating and analyzing human hereditary disease.  The American Journal of Human Genetics, 83(5):610–615, 2008.
  * [26] Wenguan Wang, Xiankai Lu, Jianbing Shen, David J. Crandall, and Ling Shao.  Zero-shot video object segmentation via attentive graph neural networks, 2019.
  * [27] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst.  Convolutional neural networks on graphs with fast localized spectral filtering, 2016.
  * [28] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks, 2016\.
  * [29] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio.  Graph attention networks, 2017.
  * [30] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?, 2018.
  * [31] William L. Hamilton, Rex Ying, and Jure Leskovec.  Inductive representation learning on large graphs, 2017.
  * [32] Felix Wu, Tianyi Zhang, Amauri Holanda de Souza Jr., Christopher Fifty, Tao Yu, and Kilian Q. Weinberger.  Simplifying graph convolutional networks, 2019.
  * [33] Jian Du, Shanghang Zhang, Guanhang Wu, Jose M. F. Moura, and Soummya Kar.  Topology adaptive graph convolutional networks, 2017.
  * [34] Jonathan Masci, Davide Boscaini, Michael M. Bronstein, and Pierre Vandergheynst.  Geodesic convolutional neural networks on riemannian manifolds, 2015.
  * [35] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding, 2018.
  * [36] Bryan Perozzi, Rami Al-Rfou, and Steven Skiena.  Deepwalk: Online learning of social representations, 2014.