I recently created a tutorial for the Harvard edX PH525.8x Case Study on DNA Methylation Analysis. This tutorial is a step-by-step guide of starting with raw Whole Genome Bisulfite Sequencing (WGBS) reads (SRA or fastq files), apply quality control filters, mapping mapping the filtered reads using Bismark and extracting the methylation calls for downstream analyses such as importing the methylation calls in R to find differently methylated CpGs or differentially methylated regions (see the bsseq R/Bioconductor package for more information).
As an example, we extracted six paired tumor-normal WGBS samples from Ziller et al. (2013) (PMID: 23925113). We have provided the coverage files (only chromosome 22) produced by Bismark in the colonCancerWGBS Github repository for others to use as a data example in R.
I hope others find the tutorial useful for the analysis of their own data!
Showing posts with label whole-genome sequencing. Show all posts
Showing posts with label whole-genome sequencing. Show all posts
Tuesday, May 5, 2015
Wednesday, June 19, 2013
Getting Started with SAMtools and Rsamtools
Have you heard of this thing called 'samtools', but not really sure what it is or how to get started? Maybe you are a postdoc in a wet-lab that is now starting to use whole-genome and whole-exome sequencing, but there isn't an in-house bioinformatician to help sift though the computational side of the project? If you're like me (a statistician by training and working with genomic data), learning about tools to handle next-generation sequencing data is essential. This tutorial is meant give a little background about samtools and then provide some technical (and hopefully helpful) information on the purpose of samtools, its relation to Rsamtools and examples of how to utilize both.
Beginning with SAMtools
What does SAM mean? How is SAM different from BAM?
The Sequence Alignment/Map (SAM) format was developed to store aligned nucleotide sequence reads produced from genome aligners (e.g. BWA) in a somewhat human readable tab-deliminted text file. Practically speaking, this format has become the standard alignment format widely used by groups such as 1000 Genomes Project. The Binary Alignment/Map (BAM) format is a binary version of SAM which contains the same information and simply compresses the data to save storage space (and BAM files are faster to manipulate). In the SAM format, there is a header section (all header rows begin with '@') and an alignment section (all alignment rows contain 11 mandatory fields + optional fields).
What is SAMtools and where did SAMtools come from?
The SAM/BAM format and SAMtools were published by Li et al. 2009. The purpose of samtools is provide a way to manipulate and interact with the aligned sequence reads in SAM/BAM format. Samtools does not work with the entire file of aligned reads, but rather works in a stream which allows Unix commands to be used in conjunction with the SAM files.
How do you install samtools?
The easiest way to do this to check out the most recent samtools source code from the github project page. If you're not familiar with git or github, I would suggest this wonderful tutorial by Karl Broman. The best description of github I've heard is this tutorial, Karl writes "Github is like facebook for programmers", which couldn't be more true. So, if you've set yourself up with a github account, then simply check out the lastest samtools source code
git clone git://github.com/samtools/samtools.git
cd samtools
make
make razip
This last line is to compile razip which is a tool to compress files and still make them randomly accessible by samtools (not normally possible with gzip files). Finally, the last step is to make samtools accessible to all the users on the machine which I extracted from this tutorial.
sudo cp ~/current/working/directory/samtools /usr/local/samtools-v0.1.19
sudo ln -s /usr/local/samtools-v0.1.19 /usr/local/samtools
The second line makes a symbolic link to the latest version directory. When a new version of samtools becomes available (e.g. samtools-v0.1.20), you can install it in the same directory and then just update the /usr/local/samtools symbolic link. You can check to make sure samtools was installed in your $PATH with the following command:
which samtools
If nothing happens, samtools was not properly put in your $PATH.
Examples of how to use SAMtools
There are many tutorials on samtools. The standard format of samtools requires a command and then allows for additional options.
samtools [command] <options>
If you type just 'samtools' in the terminal window, a list of available command (and their descriptions) will be provided. We begin by creating a folder which will contain our SAM file and reference genome in fasta format which we used to align the reads.
mkdir MyProject
cp ~/current/working/directory/test.sam MyProject
cp ~/current/working/directory/ref.fa MyProject
samtools view
To convert from SAM to BAM format (option -b) when header is available and input is SAM format (option -S)
samtools view -Sb test.sam > test.bam
To convert from SAM to BAM format when header is not available, but ref.fa was used to map the reads
samtools view -b ref.fa test.sam > test.bam
To filter out unmapped reads from BAM files (option -F)
samtools view -bF test.bam > mappedtest.bam
samtools faidx ref.fa
Examples of using SAM/BAM files
'samtools mpileup' and 'bcftools view'
To create a Binary Call Format (bcf) file (option -u) by mapping the bases using the indexed reference genome (option -f) and call genomic variants
samtools mpileup -u -f ref.fa test.sorted.bam > test.bcf
To convert the bcf file to a human readable Variant Call Format (vcf) file
bcftools view -v -g test.bcf > test.vcf
From here, you can also check out the Integrative Genomics Viewer (IGV) to visually look at the aligned reads or try importing the BAM files into tools like Rsamtools.
Examples of how to use Rsamtools
The purpose of Rsamtools is to provide an interface between R and BAM files produced by the tools samtools, bcftools, and tabix (not discussed here). To install Rsamtools in R, use
source("http://bioconductor.org/biocLite.R")
biocLite("Rsamtools")
The main purpose is to import (indexed!) BAM files into R using the scanBam() function. The vignette for Rsamtools can be found here. The entire BAM file should not be read into R, but rather just the portion of the genome of interest using the what and which arguments in ScanBamParam(which=which, what=what). There are other packages which can read in BAM files include ShortRead and GenomicRanges.
To extract the header information, use scanBamHeader(). Use filterBam() to filter reads from BAM file according to the criteria defined in ScanBamParam().
Another great use of Rsamtools is to access multiple BAM files using the BamViews class in Rsamtools. This allows you to obtain metadata by 'viewing' the BAM files rather than importing each BAM individually.
Beginning with SAMtools
What does SAM mean? How is SAM different from BAM?
The Sequence Alignment/Map (SAM) format was developed to store aligned nucleotide sequence reads produced from genome aligners (e.g. BWA) in a somewhat human readable tab-deliminted text file. Practically speaking, this format has become the standard alignment format widely used by groups such as 1000 Genomes Project. The Binary Alignment/Map (BAM) format is a binary version of SAM which contains the same information and simply compresses the data to save storage space (and BAM files are faster to manipulate). In the SAM format, there is a header section (all header rows begin with '@') and an alignment section (all alignment rows contain 11 mandatory fields + optional fields).
What is SAMtools and where did SAMtools come from?
The SAM/BAM format and SAMtools were published by Li et al. 2009. The purpose of samtools is provide a way to manipulate and interact with the aligned sequence reads in SAM/BAM format. Samtools does not work with the entire file of aligned reads, but rather works in a stream which allows Unix commands to be used in conjunction with the SAM files.
How do you install samtools?
The easiest way to do this to check out the most recent samtools source code from the github project page. If you're not familiar with git or github, I would suggest this wonderful tutorial by Karl Broman. The best description of github I've heard is this tutorial, Karl writes "Github is like facebook for programmers", which couldn't be more true. So, if you've set yourself up with a github account, then simply check out the lastest samtools source code
git clone git://github.com/samtools/samtools.git
cd samtools
make
make razip
This last line is to compile razip which is a tool to compress files and still make them randomly accessible by samtools (not normally possible with gzip files). Finally, the last step is to make samtools accessible to all the users on the machine which I extracted from this tutorial.
sudo cp ~/current/working/directory/samtools /usr/local/samtools-v0.1.19
sudo ln -s /usr/local/samtools-v0.1.19 /usr/local/samtools
The second line makes a symbolic link to the latest version directory. When a new version of samtools becomes available (e.g. samtools-v0.1.20), you can install it in the same directory and then just update the /usr/local/samtools symbolic link. You can check to make sure samtools was installed in your $PATH with the following command:
which samtools
If nothing happens, samtools was not properly put in your $PATH.
Examples of how to use SAMtools
There are many tutorials on samtools. The standard format of samtools requires a command and then allows for additional options.
samtools [command] <options>
If you type just 'samtools' in the terminal window, a list of available command (and their descriptions) will be provided. We begin by creating a folder which will contain our SAM file and reference genome in fasta format which we used to align the reads.
mkdir MyProject
cp ~/current/working/directory/test.sam MyProject
cp ~/current/working/directory/ref.fa MyProject
samtools view
To convert from SAM to BAM format (option -b) when header is available and input is SAM format (option -S)
samtools view -Sb test.sam > test.bam
To convert from SAM to BAM format when header is not available, but ref.fa was used to map the reads
samtools view -b ref.fa test.sam > test.bam
To filter out unmapped reads from BAM files (option -F)
samtools view -bF test.bam > mappedtest.bam
samtools sort
To sort BAM file by chromosomal coordinates (default) or read names (option -n) and create a test.sorted.bam file
samtools sort test.bam test.sorted
An example of using a Unix pipe to go from the SAM format to the sorted BAM format
samtools view -Sb test.sam | samtools sort test.sorted
samtools index
To create an indexed bam file (.bai) from the test.sorted.bam file
samtools index test.sorted.bam
samtools faidx
To create an indexed reference genome in FASTA format (creates a .fai file)samtools faidx ref.fa
Examples of using SAM/BAM files
'samtools mpileup' and 'bcftools view'
To create a Binary Call Format (bcf) file (option -u) by mapping the bases using the indexed reference genome (option -f) and call genomic variants
samtools mpileup -u -f ref.fa test.sorted.bam > test.bcf
To convert the bcf file to a human readable Variant Call Format (vcf) file
bcftools view -v -g test.bcf > test.vcf
From here, you can also check out the Integrative Genomics Viewer (IGV) to visually look at the aligned reads or try importing the BAM files into tools like Rsamtools.
Examples of how to use Rsamtools
The purpose of Rsamtools is to provide an interface between R and BAM files produced by the tools samtools, bcftools, and tabix (not discussed here). To install Rsamtools in R, use
source("http://bioconductor.org/biocLite.R")
biocLite("Rsamtools")
The main purpose is to import (indexed!) BAM files into R using the scanBam() function. The vignette for Rsamtools can be found here. The entire BAM file should not be read into R, but rather just the portion of the genome of interest using the what and which arguments in ScanBamParam(which=which, what=what). There are other packages which can read in BAM files include ShortRead and GenomicRanges.
To extract the header information, use scanBamHeader(). Use filterBam() to filter reads from BAM file according to the criteria defined in ScanBamParam().
Another great use of Rsamtools is to access multiple BAM files using the BamViews class in Rsamtools. This allows you to obtain metadata by 'viewing' the BAM files rather than importing each BAM individually.
Thursday, October 18, 2012
UCSC Genome Browser: A few useful tips
Most of my research revolves around analyzing next-generation sequencing data such as whole-exome sequencing data. As a statistician, I always appreciate finding useful bioinformatic tricks/tips from various tools that make me more efficient. Here are three examples of using the UCSC Genome Browser that I've found helpful.
Tip #1: How do you find a list of chromosome positions given a list of dbSNP identifiers? (Taken from the Guide to the UCSC Genome Browser FAQ by Nature)
Use the 'Variation and Repeats' group in Table Browser and the SNPs track of choice. Just specify the genome (e.g. Human) and assembly (e.g. hg19). For 'Region', if you want the chromosomal positions for a specific regions, click position and specify the region OR click genome and upload a list of dbSNP identifiers. Finally, choose your output format (e.g. GTF, BED) and click 'get output'.
Tip #2: How do you lift over a set of genomic coordinates from hg18 to hg19?
Use the Batch Coordinate Conversion (liftOver) tool in Utilities. Selected Original and New assemblies. Upload the original genomic coordinates (in a BED format) and submit.
Tip #3: How can you extract data from the UCSC browser and use it in R?
For this, we need to install the package rtracklayer. Here is an example on how to extract recombination rates:
library(rtracklayer)
my.session <- browserSession()
genome(my.session) <- "hg19"
recomb.rates <- getTable(ucscTableQuery(my.session, "recombRate"))
Tip #1: How do you find a list of chromosome positions given a list of dbSNP identifiers? (Taken from the Guide to the UCSC Genome Browser FAQ by Nature)
Use the 'Variation and Repeats' group in Table Browser and the SNPs track of choice. Just specify the genome (e.g. Human) and assembly (e.g. hg19). For 'Region', if you want the chromosomal positions for a specific regions, click position and specify the region OR click genome and upload a list of dbSNP identifiers. Finally, choose your output format (e.g. GTF, BED) and click 'get output'.
Tip #2: How do you lift over a set of genomic coordinates from hg18 to hg19?
Use the Batch Coordinate Conversion (liftOver) tool in Utilities. Selected Original and New assemblies. Upload the original genomic coordinates (in a BED format) and submit.
Tip #3: How can you extract data from the UCSC browser and use it in R?
For this, we need to install the package rtracklayer. Here is an example on how to extract recombination rates:
library(rtracklayer)
my.session <- browserSession()
genome(my.session) <- "hg19"
recomb.rates <- getTable(ucscTableQuery(my.session, "recombRate"))
Friday, August 31, 2012
biomaRt: Find gene name using chromosome number and position
In a previous post, I gave a few examples of using biomaRt in R. This is a continuation giving another useful example of using biomaRt: How to obtain gene names (e.g. HGNC) or really any information in listAttributes() function using only chromosome number and chromosome position. I was interested in obtaining the gene names for a set of mutations and decided to use biomaRt. I created a tab-delimited file called 'positions.txt' containing three columns. The first contained the chromosome number, followed by the start and end position (in this case they were the same). The following code identifies what gene the chromosome position is in and reports the HGNC gene symbol.
# Load the library
# Load the library
library(biomaRt)
# Define biomart object
mart <- useMart(biomart="ensembl", dataset="hsapiens_gene_ensembl")
# Gives a list of all possible annotations; Currently there are 1668 listed
listAttributes(mart)
# Gives a list of all filters or criteria to search by; Currently there are 333 listed
# I chose to filter by: chromosome_name, start, end
listFilters(mart)
# Read in tab-delimited file with three columns: chromosome number, start position and end position
positions <- read.table("positions.txt")
# Extract HGNC gene symbol
results <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "start_position", "end_position"), filters = c("chromosome_name", "start", "end"), values = list(positions[,1], positions[,2], positions[,3]), mart = mart)
Friday, July 13, 2012
War on Cancer
As a statistician applying my research to cancer genetics it is always encouraging to see stories like this that make me feel sometimes all this work isn't for nothing.
The article published last week in the NY Times describes the story of a woman who was diagnosed with a rare form of lymphoma in which white blood cells (T cells) become cancerous and move to the skin. Faced with a disease with no cure and no standard treatment, she was able to keep the cancer at bay using chemotherapy for five years. At that point, her health took a turn for the worse and her son who worked at Illumina, which produces some of the latest sequencing technology, decided to quit his job and help his mom full-time by sequencing her genome. After sequencing her normal and tumor DNA, they found 18K differences (or mutations) with no known significance for disease. In the analysis of her DNA, the researchers found two genes fused together forcing the growth signals in cancer cells to be reversed: the signal to stop was forcing the cells to grow and the signal to grow was forcing the cells to stop growth. The researchers decided to give her a new melanoma drug ipilimumab which forces normal T cells to grow, in hopes of stopping the growth of the tumor cells. The drug performed beautifully keeping the cancerous T cells at bay for 8 weeks before the cancer came back and ultimately she passed away a few weeks later.
41 years ago Richard Nixon signed the National Cancer Act of 1971 which started a "war on cancer". The goal was to find a "cure for cancer" by increasing the funding toward cancer research and find more effective treatments. As defined by PubMed "cancer is the uncontrolled growth of abnormal cells in the body". Cancer is not just one disease, in fact it is the word we have for hundreds of diseases we label as 'cancer'. When people talk about "curing cancer", it suggests that once a treatment/cure for one type of cancer is found, it should theoretically be applied to another form of cancer. In my experience this is definitely not the case. Every type of cancer is unique in its etiology and treatment. It is true that some drugs today developed for one type of cancer can be applied to another type of cancer because the target molecule may happen be the same for two different cancers, but this is a separate idea than the finding a "cure for cancer".
I work with several groups of researchers in the Texas Medical Center at University of Texas MD Anderson, Baylor College of Medicine and Texas Children's Hospital who are all trying to exactly this. They are sequencing the genomes of individuals affected by a particular disease and trying to find the causal mutation or reason for the disease. My small contribution comes in helping to analyze the results that come out of the sequencing. Bioinformatics attempts to take in the large amount of sequencing data and make sense of it. This is by no means easy and will take much longer to properly analyze the data than to just sequence it. But, even with these small steps I look forward to more success stories like these and seeing further progress being made in the war on cancer.
The article published last week in the NY Times describes the story of a woman who was diagnosed with a rare form of lymphoma in which white blood cells (T cells) become cancerous and move to the skin. Faced with a disease with no cure and no standard treatment, she was able to keep the cancer at bay using chemotherapy for five years. At that point, her health took a turn for the worse and her son who worked at Illumina, which produces some of the latest sequencing technology, decided to quit his job and help his mom full-time by sequencing her genome. After sequencing her normal and tumor DNA, they found 18K differences (or mutations) with no known significance for disease. In the analysis of her DNA, the researchers found two genes fused together forcing the growth signals in cancer cells to be reversed: the signal to stop was forcing the cells to grow and the signal to grow was forcing the cells to stop growth. The researchers decided to give her a new melanoma drug ipilimumab which forces normal T cells to grow, in hopes of stopping the growth of the tumor cells. The drug performed beautifully keeping the cancerous T cells at bay for 8 weeks before the cancer came back and ultimately she passed away a few weeks later.
41 years ago Richard Nixon signed the National Cancer Act of 1971 which started a "war on cancer". The goal was to find a "cure for cancer" by increasing the funding toward cancer research and find more effective treatments. As defined by PubMed "cancer is the uncontrolled growth of abnormal cells in the body". Cancer is not just one disease, in fact it is the word we have for hundreds of diseases we label as 'cancer'. When people talk about "curing cancer", it suggests that once a treatment/cure for one type of cancer is found, it should theoretically be applied to another form of cancer. In my experience this is definitely not the case. Every type of cancer is unique in its etiology and treatment. It is true that some drugs today developed for one type of cancer can be applied to another type of cancer because the target molecule may happen be the same for two different cancers, but this is a separate idea than the finding a "cure for cancer".
I work with several groups of researchers in the Texas Medical Center at University of Texas MD Anderson, Baylor College of Medicine and Texas Children's Hospital who are all trying to exactly this. They are sequencing the genomes of individuals affected by a particular disease and trying to find the causal mutation or reason for the disease. My small contribution comes in helping to analyze the results that come out of the sequencing. Bioinformatics attempts to take in the large amount of sequencing data and make sense of it. This is by no means easy and will take much longer to properly analyze the data than to just sequence it. But, even with these small steps I look forward to more success stories like these and seeing further progress being made in the war on cancer.
Friday, June 1, 2012
Private Genetics Company Announces First Mutation Patent
Last fall when I attended the American Society of Human Genetics in Montreal, there was a buzz about the company genetics company 23andMe. With the cost of exome sequencing falling, this company was offering the opportunity for anyone to have their exome sequenced quite inexpensively. On Monday 23andMe announced it's first patent titled "Polymorphisms Related to Parkinson's Disease". They discovered the G2019S variant in the SGK1 gene and suggest it may be protective against the disease. Interestingly, they say "our patent is an important step in ensuring that we've done all we can towards successful translation of this discovery" and they "want to those discoveries to move from the realm of academic publishing to the world of impacting lives by preventing, treating or curing disease". The response from the bloggers has been quite controversial. At ASHG the idea of patenting genes or mutations was discussed in an open forum which was equally controversial. Nature posted a blog about it yesterday which caught my eye.
The main problem with this announcement is the company is advertising the idea of patenting a gene or mutation as being beneficial because it will give individuals access to their genomes. In reality, these are two separate issues. The idea of patenting a gene or mutation is beneficial for the company to make money. Giving people access to their genome is a nice thought, but as a researcher who works in the field of interpreting the variants discovered in exomes or genomes, there is not a sufficient amount of information available yet to properly interpret the variants. We are just now discovering ways to efficiently sequence genomes and only at the very beginning of interpreting them. There is this idea of a $1,000 genome, but $1 million interpretation. I'm curious to see how researchers in both academia and industry respond to this new world of patenting genes and mutations.
The main problem with this announcement is the company is advertising the idea of patenting a gene or mutation as being beneficial because it will give individuals access to their genomes. In reality, these are two separate issues. The idea of patenting a gene or mutation is beneficial for the company to make money. Giving people access to their genome is a nice thought, but as a researcher who works in the field of interpreting the variants discovered in exomes or genomes, there is not a sufficient amount of information available yet to properly interpret the variants. We are just now discovering ways to efficiently sequence genomes and only at the very beginning of interpreting them. There is this idea of a $1,000 genome, but $1 million interpretation. I'm curious to see how researchers in both academia and industry respond to this new world of patenting genes and mutations.
Wednesday, April 11, 2012
Using BioMart and biomaRt
I've been wanting to explore the tools BioMart and the corresponding R package biomaRt which is a part of the bioconductor suite. I recently came across this blogpost explaining a bit more about the strength of the package.
After working through all this, I've quickly learned this is a very powerful tool in bioinformatics.
biomaRt is a package which interfaces with a large number of databases implemented by the BioMart suite. You don't need to know SQL, just R. Examples of BioMart database include Ensembl and HapMap. As the blogpost above says, "The concept is simple. You have a set of identifiers that describe a biological object, such as a gene. These are called filters. They have values – for example, HGNC symbols. You want to retrieve other identifiers – attributes – for your objects."
First, we must install the R library biomaRt.
source("http://bioconductor.org/biocLite.R")
biocLite("biomaRt")
biocLite("biomaRt")
We use the useMart() function to interface with a particular database. To see the available marts, use listMarts(). Within the database, we need to pick a particular dataset. You can see what datasets are available using the function listDatasets(). If we want to extract particular attributes from the database, we need to know what attributes are available. This can be found using the listAttributes() function. Finally, we use the getBM() function to actually extract the information.
Next, I will consider two examples.
Example 1: Randomly sample n = 500 HGNC gene IDs from the human genome
# Load library
library(biomaRt)
# Define biomart object
mart <- useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")
# listDatasets(mart)
# listAttributes(mart)
# Extract information from biomart
results <- getBM(attributes = c("hgnc_symbol"), mart = mart)
# Randomly sample the gene name list
N <- 500
sample.hgnc <- sample(results$hgnc_symbol,N)
sample.hgnc <- sample(results$hgnc_symbol,N)
Sample results
> head(sample.hgnc)
> head(sample.hgnc)
[1] "LINC00293" "C6orf223" "PRMT5-AS1" "SYT11" "FLNB" "SNORA49"
Example 2: Given a list of REFSEQ IDs, convert gene IDs to HGNC IDs or Uniprot Swissprot IDs
library(biomaRt)
# Define biomart object
# Define biomart object
mart <- useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")
# Read in file with gene names
genes <- read.csv("refseq.csv")
# Extract information from biomart
results <- getBM(attributes = c("refseq_mrna", "hgnc_symbol"), filters = "refseq_mrna", values = genes[,1], mart = mart)
results <- getBM(attributes = c("refseq_mrna", "uniprot_swissprot"), filters = "refseq_mrna", values = genes[,1], mart = mart)
# see uniqueRows = TRUE/FALSE to return unique list of IDs or not
# Match the RefSeq names with the Uniprot names
# Match the RefSeq names with the Uniprot names
matched <- match(genes[,1], results[,2])
cbind(genes,results[matched,2])
cbind(genes,results[matched,2])
Sample results
refseq Uniprot
1 NM_023018 O95544
2 NM_178545 Q8NDY8
3 NM_033467 Q495T6
4 NM_004402 O76075
5 NM_018198 Q9NVH1
6 NM_018198 Q9NVH1
Tuesday, April 3, 2012
Eliminating Highly Polymorphic Genes from Whole Exome Sequencing
Whole exome sequencing yields thousands of variants for each individual sequenced. Many of these variants are in genes that are highly polymorphic or in regions that do not sequence well and therefore would not be of interest when searching for a disease-causing variant. This is because any frequently mutated gene containing many deleterious variants will have a low probability of containing the disease-causing mutation. Also, if the gene is in a region that does not sequence well then a high number of variants will be often reported. When you find such variants, this would be considered a false positive. There is a big interest in detecting these false positive signals and eliminating these variants from the list reported back from whole exome sequencing.
A paper was just published this month in Human Mutation describing a way to do this. They published several lists of genes for researchers to use in their own projects. The researchers hypothesize these genes will not contain the disease-causing mutation. This list of genes could be an incredibly useful tool to filter out highly polymorphic genes or genes that simply do not sequence well.
A paper was just published this month in Human Mutation describing a way to do this. They published several lists of genes for researchers to use in their own projects. The researchers hypothesize these genes will not contain the disease-causing mutation. This list of genes could be an incredibly useful tool to filter out highly polymorphic genes or genes that simply do not sequence well.
Friday, March 30, 2012
Policy on 'Secondary findings' from Whole Genome Sequencing in Clinical Tests
The American College of Medical Genetics and Genomics (ACMG) is having their annual meeting this week in North Carolina. One of the major discussion points is: when a patient has their genome sequenced to look for disease-causing mutations for a specific disease in a clinical setting, what do you do with the 'secondary-findings' or other mutations unrelated to the disease in question that have been found? This is an incredibly difficult and convoluted question to answer.
For example in the clinical setting, say a patient's genome is being sequenced to test if their genome contains mutations related to high cholesterol, but in the process other mutations come back positive for Alzhimer's. Should a patient be informed of the 'secondary-finding' information? What if it were a child? Should the child know at a young age that they have a high predisposition to Alzhimer's?
In a research setting, there currently exist hundreds of large sequencing studies which sequence the genomes of many individuals suffering from a particular disease in an effort to study the etiology of that disease. When patients participate in these sequencing projects, thousands of mutations are often found which may or may not be related to a wide spectrum of diseases. When researchers find mutations related to other diseases, should the researchers be responsible of reporting the information to the patient? If the individual's genome is sequenced a second time at a later point in the future and mutations related to diseases that were not known before, but are now known are found, should the researcher be responsible of tracking down the individual to inform them? If a patient was informed at one point in time to have a deleterious mutation, but in the future that mutation is no longer considered to be deleterious, what should happen? At the American Society of Human Genetics (ASHG) annual meeting this fall in Montreal, I attended a similar forum that discussed many of these questions. The conversation can only be described as "intense and very heated". There were individuals who were adamantly in support of informing patients of secondary-findings and individuals who were adamantly against it in both the research and clinical-based setting.
The ACMG is releasing a policy statement which will be finalized this summer in support of reporting secondary findings to patients in the clinical setting. The policy says only disease-causing mutations with a high-prevelance for a treatable condition will be included for this clinical-based testing. Mutations for diseases with no known treatments will not be included in the list. I will be interested to see how we as a society decide to deal with all the other issues that will come out of this policy. A few of the issues include: How we will relay the information to the patients? Who is responsible to relay the information? Who will help the patients interpret these variants? Who is responsible for updating the patient on new information in the future? Of course there are also the legal issues related to the patient's privacy?
For example in the clinical setting, say a patient's genome is being sequenced to test if their genome contains mutations related to high cholesterol, but in the process other mutations come back positive for Alzhimer's. Should a patient be informed of the 'secondary-finding' information? What if it were a child? Should the child know at a young age that they have a high predisposition to Alzhimer's?
In a research setting, there currently exist hundreds of large sequencing studies which sequence the genomes of many individuals suffering from a particular disease in an effort to study the etiology of that disease. When patients participate in these sequencing projects, thousands of mutations are often found which may or may not be related to a wide spectrum of diseases. When researchers find mutations related to other diseases, should the researchers be responsible of reporting the information to the patient? If the individual's genome is sequenced a second time at a later point in the future and mutations related to diseases that were not known before, but are now known are found, should the researcher be responsible of tracking down the individual to inform them? If a patient was informed at one point in time to have a deleterious mutation, but in the future that mutation is no longer considered to be deleterious, what should happen? At the American Society of Human Genetics (ASHG) annual meeting this fall in Montreal, I attended a similar forum that discussed many of these questions. The conversation can only be described as "intense and very heated". There were individuals who were adamantly in support of informing patients of secondary-findings and individuals who were adamantly against it in both the research and clinical-based setting.
The ACMG is releasing a policy statement which will be finalized this summer in support of reporting secondary findings to patients in the clinical setting. The policy says only disease-causing mutations with a high-prevelance for a treatable condition will be included for this clinical-based testing. Mutations for diseases with no known treatments will not be included in the list. I will be interested to see how we as a society decide to deal with all the other issues that will come out of this policy. A few of the issues include: How we will relay the information to the patients? Who is responsible to relay the information? Who will help the patients interpret these variants? Who is responsible for updating the patient on new information in the future? Of course there are also the legal issues related to the patient's privacy?
Wednesday, March 28, 2012
Live Forum Tomorrow on 'Big Data' from the White House
Tomorrow afternoon the White House will be hosting a 90 minute forum on the Challenges and Opportunities in Big Data! Tune in at 2pm Eastern live on Thursday March 29th to see leaders in academia, industry, and the heads of these governmental agencies OSTP, NSF, NIH, DoE, DoD, DARPA and USGS. A blogpost from R-bloggers.com suggested even though it is of interest in how to store large data, what's more interesting is how to infer information from large data. This is exactly the question that genetics and genomics is asking with next-generation sequencing. We are now at a point that sequencing someone's genome is cheap. Interpreting the variants from someone's genome is the million, no billion dollar question. As a statistician, I'm curious to see the government's stance on analyzing not only large data coming out of companies such as Amazon, Google, Netflix, but also genome data (hopefully).
Friday, February 17, 2012
Loss-of-function mutations
A study by the Welcome Trust Sanger Institute and Yale University released in Science this month set out to determine on average how many genuine loss-of-function mutations do humans carry and how many genes are inactivated because of the mutations. Depending on what definition you use, humans carry ~20,000 genes. These loss-of-function mutations cause the protein to lose its structure or function. For example, one of the most common cancer genes, TP53, is called a tumor-suppressing gene because it controls the cell cycle. When TP53 is mutated, tumor cells can replicate uncontrollably because TP53 has lost its ability to control (or suppress) the cell cycle properly.
Using the three pilot phases of the 1000 Genomes Project, the researchers suggest humans carry ~100 loss-of-function (or deleterious) mutations and ~20 genes that have been inactivated (that's ~.1% of your genes)! This is such an interesting topic because up till now whenever researchers have found these loss-of-function mutations, they normally assumed it is somehow disease-causing. This is no longer the case. This news article from GenomeWeb states "as more and more apparently healthy individuals have their genomes and exomes sequenced, he added, investigators have unearthed a raft of apparent loss-of-function variants that are both intriguing and puzzling. " The article in Science is suggesting that we should expect humans to have a given number loss-of-function mutations (~100). What is still unclear is how to differentiate between the loss-of-function mutations that are disease-causing and the ones that are more benign. As personalized medicine is becoming an increasingly important topic, this type of research will be critical when whole-genome sequencing becomes cost effective.
Using the three pilot phases of the 1000 Genomes Project, the researchers suggest humans carry ~100 loss-of-function (or deleterious) mutations and ~20 genes that have been inactivated (that's ~.1% of your genes)! This is such an interesting topic because up till now whenever researchers have found these loss-of-function mutations, they normally assumed it is somehow disease-causing. This is no longer the case. This news article from GenomeWeb states "as more and more apparently healthy individuals have their genomes and exomes sequenced, he added, investigators have unearthed a raft of apparent loss-of-function variants that are both intriguing and puzzling. " The article in Science is suggesting that we should expect humans to have a given number loss-of-function mutations (~100). What is still unclear is how to differentiate between the loss-of-function mutations that are disease-causing and the ones that are more benign. As personalized medicine is becoming an increasingly important topic, this type of research will be critical when whole-genome sequencing becomes cost effective.
Subscribe to:
Posts (Atom)