multiple myeloma is a malignant b - cell disorder characterized by proliferation of atypical plasma cells in bone marrow   with or without the presence of monoclonal immunoglobulin protein in serum and/or urine .
multiple myeloma has correlation with plasmacytoma , which is a mass of plasma cells found outside of bone marrow   that needs medical intervention with radiotherapy   or chemotherapy .
while multiple myeloma frequently accompanies plasmacytoma at the time of diagnosis , plasmacytoma precedes multiple myeloma in some cases .
the disease entity called solitary plasmacytoma exists in 4% of plasma cell tumors [ 5 , 6 ] , and approximately 4050% of patients with solitary plasmacytoma will develop multiple myeloma .
hence , plasmacytoma is an early form or an accompanying disease of myeloma , and the data regarding the  clinical behavior  of plasmacytoma are quite accumulated .
however , not much is known about the cellular biology of plasmacytoma per se . fortunately , we previously succeeded in establishing two cell lines with different tropism from a single patient .
briefly , 14 weeks after injecting mononuclear cells obtained from multiple myeloma patient 's bone marrow in a nrg / scid mouse via tail vein , tumor developed at subcutis of the mouse .
the engraftment of myeloma cells into mouse bone marrow ( bm ) was also observed .
after separation and culturing cells from subcutis and bm , we succeeded in the establishing of two cell lines ( snu_mm1393_sc and snu_mm1393_bm ) from a single patient with different tropism .
snu_mm1393_bm and snu_mm1393_sc showed a high degree of resistance against bortezomib compared to u266 cell line .
snu_mm1393_bm had the greater lethality compared to snu_mm1393_sc . from genetic perspective , we believe that unique somatic mutations may allow tumor cells to adapt and survive in tumor microenvironment .
in other words , there may be specific genetic changes that contribute to tumor tropism .
hence , in this study , we tried to characterize genomic profile specific for snu_mm1393_bm and snu_mm1393_sc to understand genetic background of difference in these cell lines .
we were particularly interested in genetic changes specific for snu_mm1393_sc , because we thought these snu_mm1393_sc specific genetic changes would reveal genetic background for plasmacytoma which exhibit tropism for extramedullary space . to find these genetic changes , we performed whole exome sequencing ( wes ) using dna of both cell lines .
as it is well known , wes allows comprehensive characterization of genomic changes in individual tumors .
japan ) was used to extract dna according to the manufacturer 's recommendations . for wes , we sequenced exome using the solexa sequencing technology platform ( hiseq2000 , illumina , san diego , ca , usa ) following the manufacturer 's instructions .
we randomly sheared 3  ug of genomic dna using covaris system to generate about 150  bp inserts .
the fragmented dna was end - repaired using t4 dna polymerase and klenow polymerase , and illumina paired - end adaptor - oligonucleotides were ligated to the sticky ends .
we analyzed the ligation mixture by electrophoresis on an agarose gel , sliced and purified fragments with 200250  bp sizes .
purified dna library was hybridized with sureselect human all exon v4 probes set ( agilent , santa clara , ca , usa ) to capture 50 mb targeted exons following manufacturer 's instruction .
we prepared the hiseq2000 paired - end flow cell to the manufacturer 's protocol using captured exome library .
clusters of pcr colonies were then sequenced on the hiseq2000 platform using recommended protocols from the manufacturer .
fastq files were aligned to human reference ( human_g1k_v 37.fasta ) by using the burrows - wheeler aligner ( bwa-0.7.5 )   to make sam file .
sortsam in picard - tools-1.68 was used to convert to bam file and sort with chromosome and went through a pcr duplicate marking process , which enables the genome analysis toolkit ( gatk-1.6.5 ) to ignore duplicates in subsequent processing .
we performed a local realignment prior to recalibration , which gives the most accurate quality scores for each sample .
for single nucleotide variant ( snv ) and small indel calling , we used varscan2 ( http://genome.wustl.edu/ ) . because germline control was absent , we used pilup2snp command for single sample calling using human_g1k_v 37 as reference genome .
we set the following options : ( 1 ) minimum coverage above value of 20 , ( 2 ) minimum variant allele frequency above value of 0.1 , and ( 3 ) default p value below 0.05 other option value set as default values . to select unique mutation , we performed comparison between two calling results . for functional annotation and prediction of variant effect , we used annovar   with polyphen   database version 2.2.2 . for comparing public data with results in this study
, we used datasets from tcga ( https://tcga-data.nci.nih.gov/tcga/tcgahome2.jsp ) , cbioportal for cancer genomics ( http://www.cbioportal.org/public-portal/ ) , and kegg database for pathway analysis ( http://david.abcc.ncifcrf.gov/ ) .
the raw data size of snu_1393mm_bm and snu_1393mm_sc was 9,090  mb and 8,979  mb , respectively .
approximately 99.00% of the targeted reads ( 165483843 reads ) were covered sufficiently to pass our thresholds for calling variants ( mapq > 20 by ngs qc toolkitv2.3 ) .
mapq distribution following that above 30 was 98.2% ( 164088367 ) , above 20 was 0.8% ( 1395476 ) , and below 20 of mapq was under 10% . for snu_1393mm_sc ,
mapq distribution following that above 30 was 98.1% ( 159871347 ) , above 20 was 0.8% ( 154084 ) , and below 20 of mapq was around 10% .
when snv calling was performed using varscan , a total of 18573 snvs were found in snu_mm1393_sc .
their distribution according to the functional consequences was as follows : 8595 ( 46.2% ) nonsynonymous , 9575 ( 51.5% ) synonymous , 68 ( 0.003% ) stop - gain , and 6 ( 0.0003% ) stop - loss . in snu_mm1393_bm , a total of 18781 snvs were found and their distribution was as follows : 8694 ( 46.2% ) nonsynonymous , 9667 ( 51.5% ) synonymous , 75 ( 0.004% ) stop - gain , and 5 ( 0.0003% ) stop - loss . as for nonsynonymous snvs
, we found 8595 nonsynonymous snvs in 4901 genes for snu_mm1393_sc , while 8694 nonsynonymous snvs in 4969 genes were found in snu_mm1393_bm .
there was overlapping of 8344 nonsynonymous snvs , and 251 nonsynonymous snvs and 350 nonsynonymous snvs were unique for snu_mm1393_sc and snu_mm1393_bm , respectively ( figures 1(a)1(c ) ) .
the rate of transversion and transition in the coding region was different between the two cell lines . while transversion was dominant event in snu_mm1393_bm cell line , transition was dominant event in snu_mm1393_sc .
absolute transversion rate was much higher in snu_mm1393_bm ( 65.5% ) than snu_mm1393_sc ( 34.0% ) ( figure 1(d ) ) .
after calling of snvs , we compared genomic signatures of snu_mm1393_sc and snu_mm1393_bm with those of tumors in public database . for this comparison
, we selected 12 nonsynonymous snvs that is unique for snu_mm1393_bm and 11 nonsynonymous snvs that is unique for snu_mm1393_sc .
these snvs were selected according to the criteria below : with the assumption that two cell lines consisted of single cell population , we selected genes with variant allele frequency between 0.4 and 0.6 . first , the frequencies of these snvs were investigated in open source data of multiple myeloma ( multiple myeloma research consortium )   using cbioportal for cancer genomics ( http://www.cbioportal.org ) . around half of snvs found in our cell lines were found with low frequency ( 0.52% ) in open source database of multiple myeloma ( table 1 ) .
then , the frequencies of these snvs were investigated in open source data of various cancers ( http://www.cbioportal.org ) .
when this analysis was performed , snvs which are unique for snu_mm1393_sc were frequently detected in melanoma ( 52.8% ) and uterine cancer ( 49.0% ) .
on the other hand , snvs those which are unique for snu_mm1393_bm were frequently found in ovarian cancer ( 74.7% ) and bladder cancer ( 72.8% ) ( figure 2 ) .
using monte carlo simulation ( ratio of dns / ds ) , we examined distribution statistics of snvs found in two cell lines , respectively .
it is believed that nonrandomness of snv distribution is related to the functional importance of those snvs . when this analysis was performed , ratio for snu_mm1393_bm was 2.8 ( p = 0.14 ) , while it was 1.1 for snu_mm1393_sc ( p = 0.07 ) . hence , snv distribution in both cell lines was random with cut - off p value of 0.05 .
our results indicated that unique nonsynonymous mutations of snu_mm1393_sc seemed biologically more neutral than those of snu_mm1393_bm although they were statistically insignificant . in kegg pathway analysis of unique somatic mutation from both cell lines ,
chemokine signaling pathway and chemokine - chemokine interaction pathway were showed in two cell lines at the same time .
while tight junction and adherens junction were shown to be involved only in snu_mm1393_sc , snu_mm1393_bm cell line showed diverse disorder in pathway such as cell cycle , wnt signaling pathway , and mapk signaling pathway ( figure 3 ) .
in this report , we analyzed genomic signature of two cell lines derived from single patient .
these cell lines have different tropism , and snu_mm1393_sc is plasmacytoma which has tropism for skin .
thus , we thought this study would reveal genetic changes that are related to skin tropism and plasmacytoma .
in fact , recent study reported the process of segregation of genetic changes in tumor cells during clonal expansion . as far as we know , this is the first wes study comparing the genomics of myeloma cells in bone marrow and plasmacytoma .
varscan2 which was used in our study for somatic variant calling is one of the most commonly used tools to detect somatic mutation along with mutect [ 12 , 17 ] .
in fact , many papers already published in the field of cancer genomics already used varscan2 [ 1820 ] . as expected , most of snvs ( 96.5% ) overlapped between the cell lines , suggesting their common originality from a single patient . and the most coding sequence snvs in both snu_mm1393_bm and snu_mm1393_sc were neutral with respect to adaptation and cancer cell growth .
this is also supported by the outcome of monte carlo simulation , where distribution of snvs did not show significant nonrandomness .
it has been suggested that bone marrow microenvironment strictly regulates the growth of cell via dynamic interplay among hematopoietic cells . and
our results indicate that only a small subset of nonsynonymous snvs in cancer are affected by selection , making it possible to interpret snv trends as reflection of underlying mutational processes . the most interesting finding in our data
we conjectured in the planning of this study that genetic change unique for snu_mm1393_sc would be related to skin tropism and formation of plasmacytoma . and ,
as expected , genomic signature that is unique to snu_mm1393_sc had more than 50% of overlapping with melanoma which is a primary skin tumor .
this finding highly coincides with our conjecture and snvs such as kiaa1199 , fry , ap3b2 , and optc may be the very gene related to skin tropism of cancers .
these 4 are the genes that are mutated in snu_mm1393_sc and are frequently found in melanoma samples .
in fact , it has been known that there are common genetics between melanoma and plasmacytoma such as cdkn2a germline mutation in prenext generation sequencing ( ngs ) era .
one more noticeable finding in our study is that the number of snvs was higher in snu_mm1393_bm than in snu_mm1393_sc and transversion rate was higher in snu_mm1393_bm than in snu_mm1393_sc .
also dynamic gene to gene interaction in snu_mm1393_bm was more complex than that in snu_mm1393_sc .
along with this phenomenon and from the previous report that the number and pattern of somatic snvs determine pathway underling cancer , we think that it would be biologically simpler to form a tumor in subcutis than to form a tumor in bone marrow in animal xenograft model used in our experiment .
we think this is related to the fact that growth and differentiation of cancer cells in bone marrow are strictly regulated by dynamic interplay among various hematopoietic cells , compared to subcutis .
in fact , we had difficulties in the analysis of wes data due to the lack of germline reference dna in this study . because we used public germline database as reference in the calling of snvs , the number of snvs was very large compared to the previous multiple myeloma genomic studies .
moreover , it is well known that the majority of mutations observed in cancer sequencing studies are believed to be passenger mutations having little impact on the cancer cell . to overcome this issue
rather , we analyzed the  pattern  of genetic changes found in our study and compared them with public database so as to find out genetic clues underlying tropism for skin and plasmacytoma . and
snvs for both snu_mm1393_bm and snu_mm1393_sc instead of whole genomic picture of snu_mm1393_bm and snu_mm1393_sc .
 
 background . 
 previously we established two cell lines ( snu_mm1393_bm and snu_mm1393_sc ) from different tissues ( bone marrow and subcutis ) of mice which were injected with single patient 's myeloma sample . 
 we tried to define genetic changes specific for each cell line using whole exome sequencing ( wes ) . 
 materials and methods . 
 we extracted dna from snu_mm1393_bm and snu_mm1393_sc and performed wes . 
 for single nucleotide variants ( snv ) calling , we used varscan2 . 
 annotation of mutation was performed using annovar . 
 results . 
 when calling of somatic mutations was performed , 68 genes were nonsynonymously mutated only in snu_mm1393_sc , while 136 genes were nonsynonymously mutated only in snu_mm1393_bm . 
 kiaa1199 , fry , ap3b2 , and optc were representative genes specifically mutated in snu_mm1393_sc . when comparison analysis was performed using tcga data , mutational pattern of snu_mm1393_sc resembled that of melanoma mostly . 
 pathway analysis using kegg database showed that mutated genes specific of snu_mm1393_bm were related to differentiation , while those of snu_mm1393_sc were related to tumorigenesis . 
 conclusion . 
 we found out genetic changes that underlie tropism of myeloma cells using wes . 
 genetic signature of cutaneous plasmacytoma shares that of melanoma implying common mechanism for skin tropism . 
 kiaa1199 , fry , ap3b2 , and optc are candidate genes for skin tropism of cancers . 
1. Introduction
2. Materials and Methods
3. Results
4. Discussion