TCGA突变数据的下载、整理和可视化

TCGA突变数据的下载、整理和可视化今天是生信星球陪你的第520天大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！花花写于2020-01-17小年快乐呀大家！今天上班最后一天，接下来开启长达半个月的年假（距离下一场讲课还有20天，期待）。我今年第一年在广东过年，不回家了，我父母会从山东老家过来和我一起过年，后天就到咯。。。你呢思维导图走起啦1.数据下载1.1突变数据TCGA的突变数据有4个软件得到的不同版本：这个可以在gdc的官网上找到，case选择KIRC，文件类型选择maf即可获得。选择mutect，就一个文件，直接点进去，download就行，下载下来只有一个tar.gz文件，解压放在工作目录下。tar-xzvffile.tar.gz解压，即可得到一个maf.gz文件。同样的筛选条件，参考/p/559d9604fcdf下载临床信息数据并整理。1mkdir

kirc-clinical2

./gdc-client.exe

download

-m

gdc_manifest.2020-01-17\

\(1\).txt

-d

kirc-clinical2.数据读取2.1突变数据使用maftools读取。1rm(list=ls())2options(stringsAsFactors

=

F)

3require(maftools)

4require(dplyr)5project='TCGA_KIRC'6laml

=

read.maf(maf

=

'TCGA.KIRC.mutect.somatic.maf.gz')7#>

-Reading8#>

-Validating9#>

-Silent

variants:

8383

10#>

-Summarizing11#>

--Mutiple

centers

found12#>

BCM;BI--Possible

FLAGS

among

top

ten

genes:13#>

TTN14#>

MUC1615#>

HMCN116#>

-Processing

clinical

data17#>

--Missing

clinical

data18#>

-Finished

in

3.750s

elapsed

(3.430s

cpu)19laml

20#>

An

object

of

class

MAF

21#>

ID

summary

Mean

Median22#>

1:

NCBI_Build

GRCh38

NA

NA23#>

2:

Center

BCM;BI

NA

NA24#>

3:

Samples

336

NA

NA25#>

4:

nGenes

9444

NA

NA26#>

5:

Frame_Shift_Del

1732

5.155

427#>

6:

Frame_Shift_Ins

1201

3.574

128#>

7:

In_Frame_Del

238

0.708

029#>

8:

In_Frame_Ins

350

1.042

030#>

9:

Missense_Mutation

12997

38.682

3631#>

10:

Nonsense_Mutation

1259

3.747

232#>

11:

Nonstop_Mutation

18

0.054

033#>

12:

Splice_Site

490

1.458

134#>

13:

Translation_Start_Site

25

0.074

035#>

14:

total

18310

54.494

4736maf_df

=

laml@data37save(laml,maf_df,file

=

'maf.Rdata')38length(unique(maf_df$Tumor_Sample_Barcode))39#>

[1]

33640length(unique(maf_df$Hugo_Symbol))41#>

[1]

9444因此，有336个病人,9444个突变基因信息。了解maf还可以用下面的几个函数：1getSampleSummary(laml)

2getGeneSummary(laml)

3getFields(laml)

2.2.临床信息将下载好的临床信息xml文件整理成一个数据框。1xmls

=

dir('kirc-clinical/',pattern

=

'*.xml$',recursive

=

T)2library(XML)3td

=

function(x){4

result

<-

xmlParse(file.path('kirc-clinical/',x))5

rootnode

<-

xmlRoot(result)6

xmldataframe

<-

xmlToDataFrame(rootnode[2])7

return(t(xmldataframe))8}910cl

=

lapply(xmls,td)11cl_df

<-

as.data.frame(t(do.call(cbind,cl)))12cl_df[1:3,1:3]13#>

additional_studies

tumor_tissue_site

histological_type14#>

1

Kidney

Kidney

Clear

Cell

Renal

Carcinoma15#>

2

Kidney

Kidney

Clear

Cell

Renal

Carcinoma16#>

3

Kidney

Kidney

Clear

Cell

Renal

Carcinoma17save(cl_df,file

=

'clinical.Rdata')3.突变数据的可视化3.1plotmafSummarymaftools自带可视化函数plotmafSummary，可以比较直观的统计maf文件的数据。1dev.off()2#>

null

device

3#>

14plotmafSummary(maf

=

laml,

rmOutlier

=

TRUE,showBarcodes

=

FALSE,5

addStat

=

'median',

dashboard

=

TRUE,

titvRaw

=

FALSE)就是将maf_df数据框做了统计，用barplot和boxplot做了可视化。3.2突变频谱图代码其实就一句！1oncoplot(maf

=

laml,

top

=

30,

fontSize

=

1)下面展开一下这个图的解读主体热图一行是一个基因，总共是9444个基因，从中截取了top30；一列是一个样本，总共是336个样本。不同颜色代表不同类型的突变。右侧条形图右侧的条形图是每个基因的突变样本数、突变类型和比例验证一下突变样本数1count(maf_df,Hugo_Symbol,sort

=

T)2#>

#

A

tibble:

9,444

x

23#>

Hugo_Symbol

n4#>

<chr>

<int>5#>

1

VHL

1696#>

2

PBRM1

1487#>

3

TTN

778#>

4

SETD2

469#>

5

BAP1

3710#>

6

MUC16

2811#>

7

MTOR

2312#>

8

KDM5C

2113#>

9

HMCN1

2014#>

10

ATM

1915#>

#

…

with

9,434

more

rows结果显示VHL在169样本中突变，样本总数336，所以是49%，以此类推条形图的颜色是突变类型，以VHL基因为例，他的突变类型分别是：1maf_df

%>%

filter(Hugo_Symbol=='VHL')

%>%2

count(Variant_Classification,sort

=

T)3#>

#

A

tibble:

7

x

24#>

Variant_Classification

n5#>

<fct>

<int>6#>

1

Missense_Mutation

607#>

2

Frame_Shift_Del

418#>

3

Nonsense_Mutation

279#>

4

Frame_Shift_Ins

2210#>

5

Splice_Site

1611#>

6

In_Frame_Del

212#>

7

Nonstop_Mutation

1顶部条形图显示每个样本里突变的基因个数，可以看到最高的是那个一枝独秀的1600多。1laml@variants.per.sample

%>%

head()2#>

Tumor_Sample_Barcode

Variants3#>

1:

TCGA-B8-4143-01A-01D-1806-10