microarray-analysis/de-tutorial.html at master · microarray/microarray-analysis

History

734 lines (709 loc) · 855 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

642

643

644

645

646

647

648

649

650

651

652

653

654

655

656

657

658

659

660

661

662

663

664

665

666

667

668

669

670

671

672

673

674

675

676

677

678

679

680

681

682

683

684

685

686

687

688

689

690

691

692

693

694

695

696

697

698

699

700

701

702

703

704

705

706

707

708

709

710

711

712

713

714

715

716

717

718

719

720

721

722

723

724

725

726

727

728

729

730

731

732

733

734

<!DOCTYPE html>

<head>

<title>Differential Expression Tutorial</title>

pre:not([class]) {

background-color: white;

}

</style>

if (window.hljs && document.readyState && document.readyState === "complete") {

window.setTimeout(function() {

hljs.initHighlighting();

}, 0);

}

</script>

</head>

<body>

.main-container {

max-width: 940px;

margin-left: auto;

margin-right: auto;

}

code {

color: inherit;

background-color: rgba(0, 0, 0, 0.04);

}

img {

max-width:100%;

height: auto;

}

</style>

<h1 class="title">Differential Expression Tutorial</h1>

<h4 class="author"><em>Mark Dunning; mark ‘dot’ dunning ‘at’ cruk.cam.ac.uk, Oscar Rueda; oscar ‘dot’ rueda ‘at’ cruk.cam.ac.uk</em></h4>

<h4 class="date"><em>Last modified: 22 Feb 2016</em></h4>

</div>

<h1>Importing the data</h1>

<p>The data from this experiment comprises nine paired tumor/normal colon tissues on Illumina HT12_v3 gene expression Beadchips. These data were generated to inform a comparison of technologies for microRNA profiling. However, we will only use the mRNA data here.</p>

<pre class="r"><code>library(GEOquery)

library(limma)

url <- "ftp://ftp.ncbi.nih.gov/pub/geo/DATA/SeriesMatrix/GSE33126/GSE33126_series_matrix.txt.gz"

filenm <- "data/GSE33126_series_matrix.txt.gz"

if(!file.exists(filenm)) download.file(url, destfile=filenm)

gse <- getGEO(filename=filenm)</code></pre>

<pre><code>## File stored at:

## /tmp/RtmpPKb33G/GPL6947.soft</code></pre>

<pre><code>## ExpressionSet (storageMode: lockedEnvironment)

## assayData: 48803 features, 18 samples

## element names: exprs

## protocolData: none

## phenoData

## sampleNames: GSM820516 GSM820517 ... GSM820533 (18 total)

## varLabels: title geo_accession ... data_row_count (31 total)

## varMetadata: labelDescription

## featureData

## featureNames: ILMN_1343291 ILMN_1343295 ... ILMN_2416019 (48803

## total)

## fvarLabels: ID nuID ... GB_ACC (30 total)

## fvarMetadata: Column Description labelDescription

## experimentData: use 'experimentData(object)'

## Annotation: GPL6947</code></pre>

<pre class="r"><code>head(exprs(gse))</code></pre>

<pre><code>## GSM820516 GSM820517 GSM820518 GSM820519 GSM820520

## ILMN_1343291 66665.3800 69404.6700 64128.0700 68943.9700 67827.2200

## ILMN_1343295 22040.1100 13046.3400 38678.9600 16641.8900 33719.8900

## ILMN_1651199 226.6081 205.4483 217.2475 229.0451 226.3029

## ILMN_1651209 278.5710 253.7044 211.8002 278.0423 259.8059

## ILMN_1651210 195.4914 195.9835 175.3356 193.9065 229.5674

## ILMN_1651221 217.2764 206.0723 219.5992 205.0462 194.7481

## GSM820521 GSM820522 GSM820523 GSM820524 GSM820525

## ILMN_1343291 71775.3000 62245.5900 69713.7000 69509.2700 68244.5900

## ILMN_1343295 18933.2900 26170.0400 9906.9130 17166.5200 12428.9500

## ILMN_1651199 203.8710 213.4431 210.4129 229.5394 212.7384

## ILMN_1651209 265.1900 321.2587 273.4458 253.6032 310.1582

## ILMN_1651210 164.0632 244.6696 190.9813 188.1039 199.3084

## ILMN_1651221 233.9729 233.2414 222.1413 244.1758 209.8247

## GSM820526 GSM820527 GSM820528 GSM820529 GSM820530

## ILMN_1343291 65427.4700 68436.5200 57608.6700 69959.7700 69509.2700

## ILMN_1343295 25297.5700 17535.6000 19749.1400 17854.2900 43670.6800

## ILMN_1651199 226.1345 232.2437 208.7316 229.2948 214.4033

## ILMN_1651209 275.0126 274.9519 250.6420 255.8540 219.5752

## ILMN_1651210 220.6229 213.3975 194.7746 173.7073 185.3380

## ILMN_1651221 227.8108 255.9940 233.5495 194.4097 218.6604

## GSM820531 GSM820532 GSM820533

## ILMN_1343291 70063.7700 69647.1700 70332.3400

## ILMN_1343295 22849.0800 23725.6600 28747.0100

## ILMN_1651199 216.6758 195.6539 252.1502

## ILMN_1651209 292.4965 253.3126 237.9844

## ILMN_1651210 174.6898 195.3534 191.9382

## ILMN_1651221 271.4561 222.7748 241.0521</code></pre>

<hr />

<h2>Q Do the data look to be normalised?</h2>

<hr />

<pre class="r"><code>exprs(gse) <- log2(exprs(gse))

boxplot(exprs(gse),outline=FALSE)</code></pre>

<p>Inspect the clinical variables</p>

<hr />

</div>

<h2>Q What column has the information about whether is sample is a tumour or normal?</h2>

</div>

<h2>Q How many do we have in each group?</h2>

<hr />

<pre><code>## title geo_accession status

## GSM820516 C1 (gene expression data) GSM820516 Public on Nov 26 2012

## GSM820517 C2 (gene expression data) GSM820517 Public on Nov 26 2012

## GSM820518 C3 (gene expression data) GSM820518 Public on Nov 26 2012

## GSM820519 C4 (gene expression data) GSM820519 Public on Nov 26 2012

## GSM820520 C5 (gene expression data) GSM820520 Public on Nov 26 2012

## submission_date last_update_date type channel_count

## GSM820516 Oct 20 2011 Nov 26 2012 RNA 1

## GSM820517 Oct 20 2011 Nov 26 2012 RNA 1

## GSM820518 Oct 20 2011 Nov 26 2012 RNA 1

## GSM820519 Oct 20 2011 Nov 26 2012 RNA 1

## GSM820520 Oct 20 2011 Nov 26 2012 RNA 1

## source_name_ch1

## GSM820516 tumor

## GSM820517 normal

## GSM820518 tumor

## GSM820519 normal

## GSM820520 tumor</code></pre>

<hr />

</div>

<h2>Q Create a design matrix for the experiment, on paper if you like</h2>

<hr />

</div>

<h1>Simple t-tests</h1>

<p>For completeness, we note that t-statistics and corresponding p-values can be generated quickly using the <code>genefilter</code> package.</p>

<pre class="r"><code>SampleGroup <- pd$source_name_ch1

library(genefilter)</code></pre>

<pre><code>##

## Attaching package: 'genefilter'

## The following object is masked from 'package:base':

## anyNA</code></pre>

<pre class="r"><code>destats <- rowttests(exprs(gse),fac=SampleGroup)

head(destats)</code></pre>

<pre><code>## statistic dm p.value

## ILMN_1343291 2.6776390 0.08354028 0.016510459

## ILMN_1343295 -3.0922257 -0.66722098 0.006992817

## ILMN_1651199 0.5595122 0.02324258 0.583560512

## ILMN_1651209 1.0935577 0.07758592 0.290340309

## ILMN_1651210 -1.5758446 -0.10374687 0.134625663

## ILMN_1651221 0.2322365 0.01441968 0.819298727</code></pre>

<ul>

<li>Can we remember the drawbacks of this approach?</li>

</ul>

<p>A useful function is the <code>model.matrix</code>, which understands the linear model syntax of R (i.e. using the tilde <code>~</code> symbol)</p>

<pre class="r"><code>design <- model.matrix(~0+SampleGroup)

colnames(design) <- c("Normal","Tumour")</code></pre>

<p>By-hand, we would do;</p>

<pre class="r"><code>design <- matrix(nrow=18,ncol=2)

design[,1] <- c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1)

design[,2] <- c(1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0)

colnames(design) <- c("Normal","Tumour")</code></pre>

<p>or a shortcut in this would use the <code>rep</code> function.</p>

<pre class="r"><code>design[,1] <- rep(c(0,1),by=9)

design[,2] <- rep(c(1,0),by=9)</code></pre>

<p>Prior to the differential expression analysis, we will <em>filter</em> the data so only the top 50% most-variable genes get analysed, which <a href="http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2906865/">has been shown</a> to increase our power to detect differential expression. We can do this using the <code>varFilter</code> function in the <code>genefilter</code> package.</p>

<pre class="r"><code>library(genefilter)

gse.expFilt <- varFilter(gse)

gse.expFilt</code></pre>

<pre><code>## ExpressionSet (storageMode: lockedEnvironment)

## assayData: 24401 features, 18 samples

## element names: exprs

## protocolData: none

## phenoData

## sampleNames: GSM820516 GSM820517 ... GSM820533 (18 total)

## varLabels: title geo_accession ... data_row_count (31 total)

## varMetadata: labelDescription

## featureData

## featureNames: ILMN_1343295 ILMN_1651228 ... ILMN_2416019 (24401

## total)

## fvarLabels: ID nuID ... GB_ACC (30 total)

## fvarMetadata: Column Description labelDescription

## experimentData: use 'experimentData(object)'

## Annotation: GPL6947</code></pre>

<p>The <code>lmFit</code> funcion is used to fit the model to the data. The output has various components, most of which you probably don’t need to know about.</p>

<pre class="r"><code>fit <- lmFit(exprs(gse.expFilt), design)

names(fit)</code></pre>

<pre><code>## [1] "coefficients" "rank" "assign"

## [4] "qr" "df.residual" "sigma"

## [7] "cov.coefficients" "stdev.unscaled" "pivot"

## [10] "Amean" "method" "design"</code></pre>

<pre class="r"><code>dim(fit$coefficients)</code></pre>

<pre class="r"><code>head(fit$coefficients)</code></pre>

<pre><code>## Normal Tumour

## ILMN_1343295 14.032346 14.699567

## ILMN_1651228 13.984486 14.329732

## ILMN_1651229 9.661445 9.528996

## ILMN_1651232 8.089572 8.115662

## ILMN_1651237 9.214117 10.153798

## ILMN_1651238 7.710813 7.545338</code></pre>

<hr />

<h2>Q. What is the interpretation of the <code>coeff</code> item in the output</h2>

<hr />

<p>Now we define the contrast</p>

<pre class="r"><code>contrasts <- makeContrasts(Tumour - Normal, levels=design)

fit2 <- contrasts.fit(fit, contrasts)

head(fit2$coeff)</code></pre>

<pre><code>## Contrasts

## Tumour - Normal

## ILMN_1343295 0.66722098

## ILMN_1651228 0.34524651

## ILMN_1651229 -0.13244960

## ILMN_1651232 0.02609049

## ILMN_1651237 0.93968052

## ILMN_1651238 -0.16547496</code></pre>

<p>Finally, apply the <em>empirical Bayes’</em> step</p>

<pre class="r"><code>fit2 <- eBayes(fit2)

fit2</code></pre>

<pre><code>## An object of class "MArrayLM"

## $coefficients

## ILMN_1343295 ILMN_1651228 ILMN_1651229 ILMN_1651232 ILMN_1651237

## 0.66722098 0.34524651 -0.13244960 0.02609049 0.93968052

## 24396 more rows ...

## $rank

## [1] 2

## $assign

## NULL

## $qr

## Normal Tumour

## [1,] -3.0000000 0.0000000

## [2,] 0.3333333 3.0000000

## [3,] 0.0000000 -0.3333333

## [4,] 0.3333333 0.1111111

## [5,] 0.0000000 -0.3333333

## 13 more rows ...

## $qraux

## [1] 1.000000 1.111111

## $pivot

## [1] 1 2

## $tol

## [1] 1e-07

## $rank

## [1] 2

## $df.residual

## [1] 16 16 16 16 16

## 24396 more elements ...

## $sigma

## ILMN_1343295 ILMN_1651228 ILMN_1651229 ILMN_1651232 ILMN_1651237

## 0.4577251 0.2678995 0.4035902 0.1808880 0.7133376

## 24396 more elements ...

## $cov.coefficients

## Contrasts

## Contrasts Tumour - Normal

## Tumour - Normal 0.2222222

## $stdev.unscaled

## ILMN_1343295 ILMN_1651228 ILMN_1651229 ILMN_1651232 ILMN_1651237

## 0.4714045 0.4714045 0.4714045 0.4714045 0.4714045

## 24396 more rows ...

## $Amean

## ILMN_1343295 ILMN_1651228 ILMN_1651229 ILMN_1651232 ILMN_1651237

## 14.365957 14.157109 9.595220 8.102617 9.683957

## 24396 more elements ...

## $method

## [1] "ls"

## $design

## Normal Tumour

## [1,] 0 1

## [2,] 1 0

## [3,] 0 1

## [4,] 1 0

## [5,] 0 1

## 13 more rows ...

## $contrasts

## Contrasts

## Levels Tumour - Normal

## Normal -1

## Tumour 1

## $df.prior

## [1] 2.839956

## $s2.prior

## [1] 0.077334

## $var.prior

## [1] 11.67797

## $proportion

## [1] 0.01

## $s2.post

## ILMN_1343295 ILMN_1651228 ILMN_1651229 ILMN_1651232 ILMN_1651237

## 0.18958760 0.07260886 0.14998900 0.03944557 0.44380322

## 24396 more elements ...

## $t

## ILMN_1343295 ILMN_1651228 ILMN_1651229 ILMN_1651232 ILMN_1651237

## 3.2506558 2.7179456 -0.7254829 0.2786695 2.9922044

## 24396 more rows ...

## $df.total

## [1] 18.83996 18.83996 18.83996 18.83996 18.83996

## 24396 more elements ...

## $p.value

## ILMN_1343295 ILMN_1651228 ILMN_1651229 ILMN_1651232 ILMN_1651237

## 0.004242497 0.013718934 0.477071595 0.783534214 0.007539433

## 24396 more rows ...

## $lods

## ILMN_1343295 ILMN_1651228 ILMN_1651229 ILMN_1651232 ILMN_1651237

## -2.271990 -3.376114 -6.317279 -6.545394 -2.816065

## 24396 more rows ...

## $F

## [1] 10.56676330 7.38722816 0.52632548 0.07765668 8.95328702

## 24396 more elements ...

## $F.p.value

## [1] 0.004242497 0.013718934 0.477071595 0.783534214 0.007539433

## 24396 more elements ...</code></pre>

<p>We usually get our first look at the results by using the <code>topTable</code> command</p>

<pre class="r"><code>topTable(fit2)</code></pre>

<pre><code>## logFC AveExpr t P.Value adj.P.Val

## ILMN_1704294 3.269052 9.833508 20.44807 2.544269e-14 6.208271e-10

## ILMN_1813704 4.749574 10.782232 15.96861 2.102340e-12 2.564960e-08

## ILMN_1724686 3.627208 10.814662 15.16574 5.197074e-12 4.227127e-08

## ILMN_1811598 -2.380786 9.150986 -13.41658 4.338822e-11 2.646790e-07

## ILMN_1763749 -5.039637 11.894505 -12.42644 1.599366e-10 7.805228e-07

## ILMN_1652431 -5.998422 11.580616 -11.94769 3.096877e-10 1.259448e-06

## ILMN_1659990 1.259224 8.898440 11.72716 4.228175e-10 1.473881e-06

## ILMN_1702858 -1.511602 9.233874 -11.12718 1.009704e-09 3.079722e-06

## ILMN_1663866 2.038634 14.040020 10.82858 1.577919e-09 4.278090e-06

## ILMN_1657435 -3.354608 10.528543 -10.63195 2.127825e-09 5.111229e-06

## B

## ILMN_1704294 21.16212

## ILMN_1813704 17.73053

## ILMN_1724686 16.97912

## ILMN_1811598 15.16284

## ILMN_1763749 14.01219

## ILMN_1652431 13.42051

## ILMN_1659990 13.13976

## ILMN_1702858 12.34868

## ILMN_1663866 11.93958

## ILMN_1657435 11.66437</code></pre>

<p>The rows in this table are ordered according the <em>B-statistic</em>, or log-odds from the empirical Bayes’ analysis. You should be cautious when interpreting these values. Although a B-statistic of 0 corresponds to a 50/50 chance of a gene being DE, this is contingent on a certain percentage of genes in the experiment being DE. In experiments such as the one we are looked at here, there are a large number of DE genes. In this case genes with B-statistic < 0 could still be DE. Thus, you should use the B-statistic as a ranking criteria and not attempt to compare the values between experiments.</p>

</div>

<h1>Including annotation</h1>

<pre class="r"><code>anno <- fData(gse.expFilt)

head(anno)[,1:5]</code></pre>

<pre><code>## ID nuID Species Source

## ILMN_1343295 ILMN_1343295 9fQSYRUddRf4Z6p6T4 Homo sapiens RefSeq

## ILMN_1651228 ILMN_1651228 fZRGweI51.AwJv7v0I Homo sapiens RefSeq

## ILMN_1651229 ILMN_1651229 9EIpq4KI64KL0R57lo Homo sapiens RefSeq

## ILMN_1651232 ILMN_1651232 cqKdrRDdYiRJTTd0aw Homo sapiens RefSeq

## ILMN_1651237 ILMN_1651237 BiXUgqHpd5URfh9LP0 Homo sapiens RefSeq

## ILMN_1651238 ILMN_1651238 WKQZZ4QIE_WAr06inU Homo sapiens RefSeq

## Search_Key

## ILMN_1343295 GAPDH

## ILMN_1651228 ILMN_992

## ILMN_1651229 ILMN_10123

## ILMN_1651232 ILMN_34533

## ILMN_1651237 ILMN_18895

## ILMN_1651238 ILMN_18299</code></pre>

<pre class="r"><code>anno <- anno[,c("Symbol","Entrez_Gene_ID","Chromosome","Cytoband")]

fit2$genes <- anno

topTable(fit2)</code></pre>

<pre><code>## Symbol Entrez_Gene_ID Chromosome Cytoband logFC

## ILMN_1704294 CDH3 1001 16 16q22.1c 3.269052

## ILMN_1813704 KIAA1199 57214 15 15q25.1b 4.749574

## ILMN_1724686 CLDN1 9076 3 3q28c 3.627208

## ILMN_1811598 ADH1B 125 4 4q23b -2.380786

## ILMN_1763749 GUCA2A 2980 1 1p34.2b -5.039637

## ILMN_1652431 CA1 759 8 8q21.2b -5.998422

## ILMN_1659990 HIG2 29923 7 7q32.1a 1.259224

## ILMN_1702858 ADHFE1 137872 8 8q13.1b -1.511602

## ILMN_1663866 TGFBI 7045 5 5q31.1f-q31.2a 2.038634

## ILMN_1657435 MT1M 4499 16 16q13b -3.354608

## AveExpr t P.Value adj.P.Val B

## ILMN_1704294 9.833508 20.44807 2.544269e-14 6.208271e-10 21.16212

## ILMN_1813704 10.782232 15.96861 2.102340e-12 2.564960e-08 17.73053

## ILMN_1724686 10.814662 15.16574 5.197074e-12 4.227127e-08 16.97912

## ILMN_1811598 9.150986 -13.41658 4.338822e-11 2.646790e-07 15.16284

## ILMN_1763749 11.894505 -12.42644 1.599366e-10 7.805228e-07 14.01219

## ILMN_1652431 11.580616 -11.94769 3.096877e-10 1.259448e-06 13.42051

## ILMN_1659990 8.898440 11.72716 4.228175e-10 1.473881e-06 13.13976

## ILMN_1702858 9.233874 -11.12718 1.009704e-09 3.079722e-06 12.34868

## ILMN_1663866 14.040020 10.82858 1.577919e-09 4.278090e-06 11.93958

## ILMN_1657435 10.528543 -10.63195 2.127825e-09 5.111229e-06 11.66437</code></pre>

</div>

<h1>Diagnostic checks</h1>

<p>At this point it is useful to take stock and do some diagnostic checks</p>

<pre class="r"><code>dotchart(exprs(gse)["ILMN_1704294",])</code></pre>

<pre class="r"><code>boxplot(exprs(gse)["ILMN_1704294",]~SampleGroup)</code></pre>

<p>To find how many genes are DE we can use the <code>decideTests</code> function, which also applies multiple-test correction.</p>

<p><em>Make a note of how many genes are DE, you will need this later</em></p>

<pre class="r"><code>decideTests(fit2)</code></pre>

<pre><code>## TestResults matrix

## ILMN_1343295 ILMN_1651228 ILMN_1651229 ILMN_1651232 ILMN_1651237

## 1 0 0 0 0

## 24396 more rows ...</code></pre>

<pre class="r"><code>table(decideTests(fit2))</code></pre>

<pre><code>##

## -1 0 1

## 1438 21057 1906</code></pre>

<pre class="r"><code>sum(abs(decideTests(fit2))==1)</code></pre>

<p>The “<em>Volcano Plot</em>” function is a common way of visualising the results of a DE analysis. The <span class="math">$x$</span> axis shows the log-fold change and the <span class="math">$y$</span> axis is some measure of statistical significance, which in this case is the log-odds, or “B” statistic. A characteristic “volcano” shape should be seen.</p>

<pre class="r"><code>volcanoplot(fit2)</code></pre>

<p>The function also has options to highlight the postions of the top N genes from the test.</p>

<pre class="r"><code>volcanoplot(fit2,highlight=10,names = fit2$genes$"Symbol")</code></pre>

</div>

<h1>Exporting the results</h1>

<pre class="r"><code>write.fit(fit2, file = "de-results.txt",adjust="BH")</code></pre>

</div>

<h1>Getting the results for a particular gene</h1>

<p>The <code>topTable</code> function by default only prints the top 10 DE genes. However, by changing the arguments we can make it print the statistics for all the genes that we test.</p>

<pre class="r"><code>testResults <- topTable(fit2, number=nrow(fit2))</code></pre>

<p>This gives us the opportunity to look-up the results for particular genes of interest. Lets say we’re interested in “OCIAD2” (for no particular reason…)</p>

<pre class="r"><code>testResults[which(testResults$Symbol == "OCIAD2"),]</code></pre>

<pre><code>## Symbol Entrez_Gene_ID Chromosome Cytoband logFC AveExpr

## ILMN_1700306 OCIAD2 132299 4 4p12a 0.7480189 13.40736

## ILMN_1772286 OCIAD2 132299 4 4p12a 0.7877357 10.92200

## t P.Value adj.P.Val B

## ILMN_1700306 4.190794 0.0005035995 0.009607765 -0.2227105

## ILMN_1772286 3.538375 0.0022183900 0.024230499 -1.6530080</code></pre>

<p>If we have a list of genes, the <code>%in%</code> function can be used.</p>

<pre class="r"><code>mylist <- c("LOC441066","ARF3","FMNL3","CSF1R","XLKD1","TTRAP","DMWD","SYNPO2L","PILRB","LAMP3")

testResults[which(testResults$Symbol %in% mylist),]</code></pre>

<pre><code>## Symbol Entrez_Gene_ID Chromosome Cytoband logFC

## ILMN_2210753 XLKD1 10894 11 11p15.4a -1.45372426

## ILMN_1785821 TTRAP 51567 6 6p22.2b -0.43214642

## ILMN_1768754 PILRB 29990 7 7q22.1c -0.37606793

## ILMN_1686623 CSF1R 1436 5 5q33.1c -0.39403423

## ILMN_2183216 TTRAP 51567 6 6p22.2b -0.34506377

## ILMN_1657955 FMNL3 91010 12 12q13.12c 0.10533704

## ILMN_1685534 PILRB 29990 7 7q22.1c 0.07702480

## ILMN_1690253 SYNPO2L 79933 10 10q22.2a 0.06322535

## ILMN_1807712 PILRB 29990 7 7q22.1c 0.08345081

## ILMN_1682938 ARF3 377 12 12q13.12a -0.07912146

## ILMN_1714352 DMWD 1762 19 19q13.32a 0.11488061

## ILMN_2170813 LAMP3 27074 3 3q27.1a 0.04895846

## ILMN_1714808 LOC441066 441066 5 5q21.1b 0.04097845

## ILMN_1723984 PILRB 29990 7 7q22.1c 0.04075150

## ILMN_2395214 FMNL3 91010 12 12q13.12c 0.02334380

## ILMN_2170814 LAMP3 27074 3 3q27.1a -0.02860983

## AveExpr t P.Value adj.P.Val B

## ILMN_2210753 8.469606 -6.23603683 5.672890e-06 0.0006207362 4.147419

## ILMN_1785821 8.700032 -1.88536236 7.489913e-02 0.2298027154 -4.905780

## ILMN_1768754 11.081650 -1.35504089 1.914305e-01 0.4149090657 -5.680934

## ILMN_1686623 11.593908 -1.25905686 2.233870e-01 0.4549592337 -5.799590

## ILMN_2183216 9.003461 -1.21319562 2.400473e-01 0.4741293394 -5.853696

## ILMN_1657955 8.010168 1.02756048 3.171695e-01 0.5584683031 -6.054871

## ILMN_1685534 8.116290 0.76377698 4.544597e-01 0.6807561543 -6.288670

## ILMN_1690253 7.719182 0.74547643 4.651835e-01 0.6894119309 -6.302515

## ILMN_1807712 8.638967 0.55453559 5.857359e-01 0.7787153453 -6.427850

## ILMN_1682938 11.506627 -0.48960363 6.300665e-01 0.8077682634 -6.462371

## ILMN_1714352 9.153734 0.43380919 6.693539e-01 0.8312757025 -6.488688

## ILMN_2170813 8.521685 0.33382435 7.421991e-01 0.8744761426 -6.528028

## ILMN_1714808 8.273997 0.25762280 7.994910e-01 0.9070871664 -6.551204

## ILMN_1723984 7.887657 0.23765047 8.147175e-01 0.9141899453 -6.556298

## ILMN_2395214 8.693602 0.12395878 9.026609e-01 0.9591042109 -6.577500

## ILMN_2170814 9.633219 -0.07762576 9.389442e-01 0.9745572234 -6.582323</code></pre>

<h2>Further filtering and ordering</h2>

<hr />

</div>

<h2>Q. Create a data frame containing genes with i) adjsuted p-value < 0.05 ii) log-fold change > 1. How many genes does this contain?</h2>

</div>

<h2>Q. Now order this data frame according to chromosome</h2>

<hr />

<pre><code>## Symbol Entrez_Gene_ID Chromosome Cytoband logFC

## ILMN_1811598 ADH1B 125 4 4q23b -2.380786

## ILMN_1763749 GUCA2A 2980 1 1p34.2b -5.039637

## ILMN_1652431 CA1 759 8 8q21.2b -5.998422

## ILMN_1702858 ADHFE1 137872 8 8q13.1b -1.511602

## ILMN_1657435 MT1M 4499 16 16q13b -3.354608

## ILMN_1743445 FAM107A 11170 3 3p14.2d -2.544128

## AveExpr t P.Value adj.P.Val B

## ILMN_1811598 9.150986 -13.41658 4.338822e-11 2.646790e-07 15.16284

## ILMN_1763749 11.894505 -12.42644 1.599366e-10 7.805228e-07 14.01219

## ILMN_1652431 11.580616 -11.94769 3.096877e-10 1.259448e-06 13.42051

## ILMN_1702858 9.233874 -11.12718 1.009704e-09 3.079722e-06 12.34868

## ILMN_1657435 10.528543 -10.63195 2.127825e-09 5.111229e-06 11.66437

## ILMN_1743445 10.102938 -10.57083 2.337004e-09 5.111229e-06 11.57787</code></pre>

<pre><code>## Symbol Entrez_Gene_ID Chromosome Cytoband logFC

## ILMN_2363634 ADHFE1 137872 8q13.1b -0.9144158

## ILMN_1668629 LOC401115 401115 4p16.3b 1.5394078

## ILMN_1849117 -0.7814479

## ILMN_1692198 LOC641765 641765 7q11.21e 0.8078182

## ILMN_1740609 CCL15 6359 17q12b -1.7684470

## ILMN_1683740 TOP1MT 116447 8q24.3f 1.0122906

## AveExpr t P.Value adj.P.Val B

## ILMN_2363634 8.438927 -7.157882 8.816662e-07 0.0001846922 5.953946

## ILMN_1668629 12.952781 7.141844 9.098006e-07 0.0001850004 5.923572

## ILMN_1849117 8.269265 -7.115754 9.575669e-07 0.0001915212 5.874083

## ILMN_1692198 8.328211 6.824160 1.707148e-06 0.0002913015 5.314140

## ILMN_1740609 9.892766 -6.694756 2.214637e-06 0.0003509049 5.061669

## ILMN_1683740 8.753496 6.593593 2.718609e-06 0.0004094863 4.862611</code></pre>

</div>

<h2>Improved analysis: Doing a paired analysis</h2>

<hr />

</div>

<h2>Q. What would the design matrix look like to incorporate patient ID into the model</h2>

</div>

<h2>Q. Create this model using the <code>model.matrix</code> function</h2>

<hr />

</div>

<h2>Q. Proceed to fit your model to the data</h2>

</div>

<h2>Q. How many genes are DE? How does this compare to the previous model?</h2>

<hr />

<pre class="r"><code>fit.better <- lmFit(exprs(gse.expFilt), design.better)

fit2.better <- eBayes(fit.better)

fit2.better$genes <- anno

topTable(fit2.better,coef=10)</code></pre>

<pre><code>## Symbol Entrez_Gene_ID Chromosome Cytoband logFC

## ILMN_1704294 CDH3 1001 16 16q22.1c 3.269052

## ILMN_1652431 CA1 759 8 8q21.2b -5.998422

## ILMN_1659990 HIG2 29923 7 7q32.1a 1.259224

## ILMN_1724686 CLDN1 9076 3 3q28c 3.627208

## ILMN_1813704 KIAA1199 57214 15 15q25.1b 4.749574

## ILMN_1740917 SCNN1B 6338 16 16p12.1c -3.081947

## ILMN_1772521 MTHFD1L 25902 6 6q25.1b-q25.1c 1.672602

## ILMN_1663866 TGFBI 7045 5 5q31.1f-q31.2a 2.038634

## ILMN_1811598 ADH1B 125 4 4q23b -2.380786

## ILMN_1735578 GUCA2B 2981 1 1p34.2b -4.018168

## AveExpr t P.Value adj.P.Val B

## ILMN_1704294 9.833508 22.75596 2.032905e-10 4.960491e-06 12.948698

## ILMN_1652431 11.580616 -18.76750 1.530057e-09 1.866747e-05 11.592372

## ILMN_1659990 8.898440 16.95215 4.405367e-09 3.555311e-05 10.806904

## ILMN_1724686 10.814662 16.33058 6.484600e-09 3.555311e-05 10.507833

## ILMN_1813704 10.782232 16.14761 7.285175e-09 3.555311e-05 10.416580

## ILMN_1740917 10.237096 -14.41607 2.338530e-08 9.510411e-05 9.473201

## ILMN_1772521 10.610849 14.16790 2.792895e-08 9.735632e-05 9.325141

## ILMN_1663866 14.040020 13.69802 3.939725e-08 1.201665e-04 9.035112

## ILMN_1811598 9.150986 -12.58647 9.291885e-08 2.305872e-04 8.294591

## ILMN_1735578 9.916240 -12.56545 9.449906e-08 2.305872e-04 8.279803</code></pre>

<pre class="r"><code>table(decideTests(fit2.better[,10]))</code></pre>

<pre><code>##

## -1 0 1

## 1657 20431 2313</code></pre>

<pre class="r"><code>sum(abs(decideTests(fit2.better[,10]))==1)</code></pre>

</div>

<h1>The estrogen dataset</h1>

<p>Recall the estrogen dataset from yesterday</p>

<pre class="r"><code>library(affy)

targetsFile <- "estrogen/estrogen.txt"

pd <- read.AnnotatedDataFrame(targetsFile,header=TRUE,sep="",row.names=1)

pData(pd)</code></pre>

<pre><code>## estrogen time.h

## low10-1.cel absent 10

## low10-2.cel absent 10

## high10-1.cel present 10

## high10-2.cel present 10

## low48-1.cel absent 48

## low48-2.cel absent 48

## high48-1.cel present 48

## high48-2.cel present 48</code></pre>

<pre class="r"><code>ER <- pData(pd)$estrogen

Time <- factor(pData(pd)$time.h)

design <- model.matrix(~ER+Time)

design</code></pre>

<pre><code>## (Intercept) ERpresent Time48

## 1 1 0 0

## 2 1 0 0

## 3 1 1 0

## 4 1 1 0

## 5 1 0 1

## 6 1 0 1

## 7 1 1 1

## 8 1 1 1

## attr(,"assign")

## [1] 0 1 2

## attr(,"contrasts")

## attr(,"contrasts")$ER

## [1] "contr.treatment"

## attr(,"contrasts")$Time

## [1] "contr.treatment"</code></pre>

<pre class="r"><code>design2 <- model.matrix(~ER*Time)

design2</code></pre>

<pre><code>## (Intercept) ERpresent Time48 ERpresent:Time48

## 1 1 0 0 0

## 2 1 0 0 0

## 3 1 1 0 0

## 4 1 1 0 0

## 5 1 0 1 0

## 6 1 0 1 0

## 7 1 1 1 1

## 8 1 1 1 1

## attr(,"assign")

## [1] 0 1 2 3

## attr(,"contrasts")

## attr(,"contrasts")$ER

## [1] "contr.treatment"

## attr(,"contrasts")$Time

## [1] "contr.treatment"</code></pre>

<pre class="r"><code>raw <-ReadAffy(celfile.path = "estrogen", filenames=rownames(pData(pd)),phenoData = pd)

raw</code></pre>

<pre><code>## Warning: replacing previous import by 'utils::tail' when loading

## 'hgu95av2cdf'</code></pre>

<pre><code>## Warning: replacing previous import by 'utils::head' when loading

## 'hgu95av2cdf'</code></pre>

<pre><code>## AffyBatch object

## size of arrays=640x640 features (19 kb)

## cdf=HG_U95Av2 (12625 affyids)

## number of samples=8

## number of genes=12625

## annotation=hgu95av2

## notes=</code></pre>

<pre><code>## Background correcting

## Normalizing

## Calculating Expression</code></pre>

<pre><code>## ExpressionSet (storageMode: lockedEnvironment)

## assayData: 12625 features, 8 samples

## element names: exprs

## protocolData

## sampleNames: low10-1.cel low10-2.cel ... high48-2.cel (8 total)

## varLabels: ScanDate

## varMetadata: labelDescription

## phenoData

## sampleNames: low10-1.cel low10-2.cel ... high48-2.cel (8 total)

## varLabels: estrogen time.h

## varMetadata: labelDescription

## featureData: none

## experimentData: use 'experimentData(object)'

## Annotation: hgu95av2</code></pre>

</div>

// add bootstrap table styles to pandoc tables

$(document).ready(function () {

$('tr.header').parent('thead').parent('table').addClass('table table-condensed');

});

</script>

(function () {

var script = document.createElement("script");

script.type = "text/javascript";

script.src = "https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";

document.getElementsByTagName("head")[0].appendChild(script);

})();

</script>

</body>

</html>

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

FilesExpand file tree

de-tutorial.html

Latest commit

History

de-tutorial.html

File metadata and controls