https://www.biostars.org/p/125444/
Puget Peak Quad Xeon Tower offers 512GB, 1TB, or even 2TB RAM. This can be configured at: https://www.pugetsystems.
http://www.penguincomputing.com/products/workstations/tempest-4400/
Showing posts with label NGS. Show all posts
Showing posts with label NGS. Show all posts
Wednesday, January 4, 2017
NGS workstation configurations
https://www.biostars.org/p/125444/
Puget Peak Quad Xeon Tower offers 512GB, 1TB, or even 2TB RAM. This can be configured at: https://www.pugetsystems.
http://www.penguincomputing.com/products/workstations/tempest-4400/
Monday, December 28, 2015
yeast ngs data resource, Waern and Snyder 2013 G3.
Extensive Transcript Diversity and Novel Upstream
Open Reading Frame Regulation in Yeast
Karl Waern*,† and Michael Snyder*,†,1
*Department of Genetics, Stanford University School of Medicine, Stanford, California 94305, and †Department of Molecular, Cellular, and Developmental Biology, Yale University, New Haven, Connecticut 06520
http://downloads.yeastgenome.org/published_datasets/Waern_2013_PMID_23390610/fastq/
http://www.g3journal.org/content/3/2/343.full.pdf+html
Some discussion on problems during using these data
http://seqanswers.com/forums/archive/index.php/t-45228.html
"Data is single-end. I just checked the protocol accompanying the data- it confirms that the reads are indeed strand-specific."
Karl Waern*,† and Michael Snyder*,†,1
*Department of Genetics, Stanford University School of Medicine, Stanford, California 94305, and †Department of Molecular, Cellular, and Developmental Biology, Yale University, New Haven, Connecticut 06520
http://downloads.yeastgenome.org/published_datasets/Waern_2013_PMID_23390610/fastq/
http://www.g3journal.org/content/3/2/343.full.pdf+html
Some discussion on problems during using these data
http://seqanswers.com/forums/archive/index.php/t-45228.html
"Data is single-end. I just checked the protocol accompanying the data- it confirms that the reads are indeed strand-specific."
Saturday, July 18, 2015
NGS and SAGE
realized today that NGS and SAGE gene expression analysis share similar aspects.
my previous SAGE modeling did not considering isoforms, though
http://samplingprocess.org/2015/01/19/measuring-differential-expression-in-the-context-of-isoforms/
my previous SAGE modeling did not considering isoforms, though
http://samplingprocess.org/2015/01/19/measuring-differential-expression-in-the-context-of-isoforms/
Thursday, December 18, 2014
NGS method, RNA seq
Method in Lei 2013, Gene, Diminishing returns in next-generation sequencing (NGS)
transcriptome data.
The sequencing files downloaded from NCBI SRA database were initially
converted from SRA format to FASTQ format using SRA toolkit
(http://www.ncbi.nlm.nih.gov/Traces/sra/?view=software). Then, the
raw data were filtered using the following criteria: (1) the number of
unknown bases (N) was no more than two for each read; and (2) the
fraction of low quality sites (Q b 5) was no more than 50% for each
read. The data that passed this quality control were then used to map
back to their respective genome sequences using bowtie2 (Langmead
and Salzberg, 2012). Only uniquely mapped reads with no more than
two mismatches were retained for further analysis. After mapping, the
counts for each gene were summarized using HTSeq (http://wwwhuber.
embl.de/users/anders/HTSeq/doc/overview.html). In the simulation,
a predetermined-sized subset of reads was randomly selected
from the original file. Using the samemapping procedure as mentioned
above, the RPKM for each gene and depth of coverage were calculated
and comparedwith those fromoriginal data. In-house Perl and R scripts
were developed for data analysis and graphing (available upon request).
Sunday, March 9, 2014
Thursday, March 6, 2014
git, 2011scengs, shell operations, aborted and then GitHub piecemeal fashion.
I had trouble to push the large data using Gui version. So, I tried to use the shell commands.
git config --global user.name "hongqin"
biolaptop17:2011scengs hongqin$ sudo curl -s -O \
http://github-media-downloads.s3.amazonaws.com/osx/git-credential-osxkeychain
biolaptop17:bin hongqin$ ll git-credential-osxkeychain
-rw-r--r-- 1 root wheel 15K Mar 6 11:11 git-credential-osxkeychain
sudo chmod u+x git-credential-osxkeychain
biolaptop17:2011scengs hongqin$ pwd
/Users/hongqin/github/2011scengs
$ git add src/*
git config --global user.name "hongqin"
biolaptop17:2011scengs hongqin$ sudo curl -s -O \
http://github-media-downloads.s3.amazonaws.com/osx/git-credential-osxkeychain
biolaptop17:bin hongqin$ ll git-credential-osxkeychain
-rw-r--r-- 1 root wheel 15K Mar 6 11:11 git-credential-osxkeychain
sudo chmod u+x git-credential-osxkeychain
biolaptop17:bin hongqin$ sudo git config --global credential.helper osxkeychain
/Users/hongqin/github/2011scengs
$ git commit -m "src"
biolaptop17:2011scengs hongqin$ git remote add origin https://github.com/hongqin/2011scengs.git
fatal: remote origin already exists.
git push --force origin master
#fatal error
#fatal error
I deleted the .git directory. Re-cloned Gibhut directory, and then move small-sized folders into the github project. Synch with Github Gui worked this time.
It turns out that Gibhub limize files size to 100M. How did I uploaded 1.6G tarfile into the Bacillus coat repository?
Check Github show that 2011scengs was not updated.
It turns out that Gibhub limize files size to 100M. How did I uploaded 1.6G tarfile into the Bacillus coat repository?
Saturday, March 1, 2014
vcf format, va
From
http://www.broadinstitute.org/igv/viewing_vcf_files
VCF stands for Variant Call Format, and this file format is used by the 1000 Genomes project to encode SNPs and other structural genetic variants. The format is further described on the 1000 Genomes project Web site. VCF calls are available at EBI / NCBI.
http://www.broadinstitute.org/software/igv/?q=download
Wednesday, March 20, 2013
NGS services (in progress)
Beckman Coulter
Genomics, $899, based on a batch of 32 samples including library construction, 2x75bp sequencing on Illumina HiSeq 2000 and primary data analysis. This means 32x900 ~ 30K per order.
Monday, November 12, 2012
botseq_hong_variant_calling.mk
# Run freebayes on bam files
#
.DELETE_ON_ERROR:
ROOT_DIR?=/ndata/botlabseq/seqshare
include ${ROOT_DIR}/src/config.mk
include ${ROOT_DIR}/src/hong_analysis/botseq_hong_mapping.mk
REFERENCE_INDEX?=${ROOT_DIR}/reference_genomes/saccharomyces_cerevisiae_s288c_saccer3/saccharomyces_cerevisiae_s288c_saccer3.fasta
HONG_VCF_DIR=${ROOT_DIR}/hong_analysis/vcf
HONG_VCF_FILES=$(addprefix ${HONG_VCF_DIR}/, $(patsubst %.bam,%.vcf,$(notdir ${HONG_SORTEDBAM_FILES})))
HONG_VCF_GZ_FILES=$(addsuffix .gz, ${HONG_VCF_FILES})
HONG_VCF_INDEX_FILES=$(addsuffix .tbi, ${HONG_VCF_GZ_FILES})
.DEFAULT_GOAL:=hongvcall
hongvcall: ${HONG_VCF_INDEX_FILES}
${HONG_VCF_DIR}:
mkdir ${HONG_VCF_DIR}
.SECONDARY: ${HONG_VCF_FILES}
${HONG_VCF_FILES}: ${HONG_VCF_DIR}/%.vcf: ${HONG_MAPPED_BAM_DIR}/%.bam | ${HONG_VCF_DIR}
${FREEBAYES} --fasta ${REFERENCE_INDEX} --pooled $< > $@
${HONG_VCF_GZ_FILES}: ${HONG_VCF_DIR}/%.gz: ${HONG_VCF_DIR}/%
${BGZIP} $<
${HONG_VCF_INDEX_FILES}: ${HONG_VCF_DIR}/%.tbi: ${HONG_VCF_DIR}/%
${TABIX} $<
hongvctest:
@echo ${HONG_VCF_INDEX_FILES}
#
.DELETE_ON_ERROR:
ROOT_DIR?=/ndata/botlabseq/seqshare
include ${ROOT_DIR}/src/config.mk
include ${ROOT_DIR}/src/hong_analysis/botseq_hong_mapping.mk
REFERENCE_INDEX?=${ROOT_DIR}/reference_genomes/saccharomyces_cerevisiae_s288c_saccer3/saccharomyces_cerevisiae_s288c_saccer3.fasta
HONG_VCF_DIR=${ROOT_DIR}/hong_analysis/vcf
HONG_VCF_FILES=$(addprefix ${HONG_VCF_DIR}/, $(patsubst %.bam,%.vcf,$(notdir ${HONG_SORTEDBAM_FILES})))
HONG_VCF_GZ_FILES=$(addsuffix .gz, ${HONG_VCF_FILES})
HONG_VCF_INDEX_FILES=$(addsuffix .tbi, ${HONG_VCF_GZ_FILES})
.DEFAULT_GOAL:=hongvcall
hongvcall: ${HONG_VCF_INDEX_FILES}
${HONG_VCF_DIR}:
mkdir ${HONG_VCF_DIR}
.SECONDARY: ${HONG_VCF_FILES}
${HONG_VCF_FILES}: ${HONG_VCF_DIR}/%.vcf: ${HONG_MAPPED_BAM_DIR}/%.bam | ${HONG_VCF_DIR}
${FREEBAYES} --fasta ${REFERENCE_INDEX} --pooled $< > $@
${HONG_VCF_GZ_FILES}: ${HONG_VCF_DIR}/%.gz: ${HONG_VCF_DIR}/%
${BGZIP} $<
${HONG_VCF_INDEX_FILES}: ${HONG_VCF_DIR}/%.tbi: ${HONG_VCF_DIR}/%
${TABIX} $<
hongvctest:
@echo ${HONG_VCF_INDEX_FILES}
Friday, November 9, 2012
Parse genetic variants from NGS data, sequenced qin-lab natural isolates
These shell script were written by Lance at the Princeton genome center.
batch code are in src/
mkdir -p /Volumes/BotlabSeq/hong_analysis/mapped_bam
bwa aln -t 2 /Volumes/BotlabSeq/reference_genomes/saccharomyces_cerevisiae_s288c_saccer3/saccharomyces_cerevisiae_s288c_saccer3.fasta /Volumes/BotlabSeq/fastq/3-001_M1-2_lib1_read-1.fastq.gz > /Volumes/BotlabSeq/hong_analysis/mapped_bam/3-001_M1-2_lib1_read-1_aln_sa.sai
... ...
#!/bin/bash
DIRECTORY=$(cd `dirname $0` && pwd)
ROOT_DIR=$(cd "${DIRECTORY}/../../../" && pwd)
STRAINS_FILE=${ROOT_DIR}/Sequencing_Consortium_Strains-MASTER.txt
FASTQ_DIR=${ROOT_DIR}/fastq
ORIGINAL_BAM_DIR=${ROOT_DIR}/original_bam
function process_batch {
HONG_FASTQ_FILES=""
ORIGINAL_BAM_FILES=""
for strain in "${strain_batch[@]}"; do
HONG_FASTQ_FILES+=`find "${FASTQ_DIR}" -name "$strain*.fastq.gz" | xargs echo`
HONG_FASTQ_FILES+=" "
HONG_ORIGINAL_BAM_FILES+=`find "${ORIGINAL_BAM_DIR}" -name "$strain*.bam" | xargs echo`
HONG_ORIGINAL_BAM_FILES+=" "
done
export HONG_FASTQ_FILES
echo $HONG_FASTQ_FILES
export HONG_ORIGINAL_BAM_FILES
echo $HONG_ORIGINAL_BAM_FILES
echo "make -n -f ${DIRECTORY}/botseq_hong_mapping.mk -j 5"
echo ""
}
strain_list=()
while IFS=$'\t' read -r -a fields
do
genus_species=${fields[1]}
individual_name=${fields[4]}
number=${fields[11]}
if [[ $number == 3-* ]]; then
echo "$number $individual_name $genus_species - SELECTED"
strain_list+=("$number")
fi
done < $STRAINS_FILE
strain_batch=()
i=0
for strain in "${strain_list[@]}"; do
strain_batch+=("$strain")
i=$i+1
if [[ $i -eq 10 ]]; then
echo ${strain_batch[@]}
process_batch
strain_batch=()
i=0
fi
done
echo ${strain_batch[@]}
process_batch
exit
#${DIRECTORY}/botseq_lte_merge_bam.py > ${DIRECTORY}/botseq_combine_bam.mk
#make -f ${DIRECTORY}/botseq_combine_bam.mk -j 4
The mapped variant are in "vcf" files.
References
http://www.broadinstitute.org/igv/book/export/html/16
File formats
http://vcftools.sourceforge.net/perl_module.html
VCF perl tools
http://sourceforge.net/projects/samtools/files/tabix/
Saturday, September 10, 2011
2011 NGS discussion on mitbbs
There are old posts from 2011
标 题: 简单介绍 Bioinformatics Tools for NGS 分析
发信站: BBS 未名空间站 (Sat Sep 10 19:30:42 2011, 美东)
本人在wet lab里面做纯数据分析，for NGS data analysis, 简单介绍一些自己接触过
，并且觉得挺有用的工具，说的有点杂，权作抛砖引玉，还请不吝赐教。
Next-Gen sequencing(NGS)和现在正在发展的3rd-gen sequencing将会在生物学研究中
被越来越广泛应用。不管你信不信，反正我信了。一是基于实验成本的降低（$1k
whole-genome sequencing is coming），越来越多的实验室可以操作；二是可以提供
相对low throughput experiment多的多的数据和信息，可以看到很多从前看不到的东
西；三是sequencer本身对测序的准确性正在逐渐提高，所以实验固有错误率降低；四
是各种算法的成熟应用，这使得很多由于实验产生的误差在出数据后通过对数据的分析
得以过滤。按照library preparation来分，NGS主要有DNA-seq和RNA-seq
DNA-seq is usually used as ChIP-seq to study transcription factor(TF)-DNA
binding, epigenetic histome modification (such as H3K4me3, H3K36me3 etc),
DNA methylation etc. and Exon-seq or whole-DNA sequencing to study genetic
variation (SNPs/Indels etc). 从bioinformatics的角度来说，DNA-seq相对比较容易
处理，简单的说，第一步先做alignment，就是把millions of short reads from the
sequencer map back to the genome。很多genome已经有相对完整的reference，这些
genome的whole sequence可以从UCSC genome browser(http://genome.ucsc.edu/)下载。先前已经有朋友贴过现有的alignment tools（ http://en.wikipedia.org/wiki/List_of_sequence_alignment_software#Short- Read_Sequence_Alignment），目前比较常用的有bowtie，bwa，maq，soap，从他们各自的网站，应该都可以下载到 executable binary或者source code。我们比较常用的是bowtie（ http://bowtie-bio.sourceforge.net/index.shtml ）和bwa（ http://bio-bwa.sourceforge.net/），简单方便速度快（paralyzed）；output files例如SAM，BAM format比较standardized，方便后续处理；他们的网站有非常详细的manual。
如果是做ChIP-seq，通常第二步是call peaks。TF的peak一般比较peaky，MACS(http://liulab.dfci.harvard.edu/MACS /00README.html)是个不错的peak caller. TF ChIP-seq，一般会在call出peak之后assign peak to gene, in order to find TF regulated genes，这个比较arbitary，overlap of the called peak with gene promoter, enhancer, or use GREAT(http://great.stanford.edu/) to assign genes. 通常TF-DNA binding都有motif，meme（http://meme.sdsc.edu/meme/intro.html）是个不错的motif caller，use meme你可以看看这些TF enriched sites是否有significant motif, you can also check whether a known motif is enriched in a certain genomic region you interested.也就是说，既可以从genomic region找motif，也可以反向从motif找potential binding sites。
Histome modification，例如H3K4me3, K3K36me3,H3K27me3,H3K27ac, etc通常不具有
motif（not expected to），usually people are interested in their position,
downstream genes,how the intensity of peak changes (from average diagram),
whether there is bivalent domain (H3K4me3, H3K27me3) etc. 整体来说，对于ChIP
-seq,可以利用open source tool自己整workflow；可以利用galaxy(http://main.g2.bx.psu.edu/)已有的workflow 和file tools（上面也有详细的tutorial）；也可以用Partek Genomic Suite （not free,licence required）.后两个适合bench worker使用，非常容易上手。
Exon-seq and whole-DNA-seq asks very different question, usually people are
interested to know genetic variations, such as SNPs, Indels, Copy number
variation and looking for reoccurrence in tumor samples. 第一步alignment几乎
没有什么区别，找genetic variation, SeqGene(http://sourceforge.net/apps/mediawiki/seqgene /index.php?title=SeqGene)是一个比较简单且实用的tool。SNPs/Indels的prediction，基本上是比较准确 的，可以说那些high quality的prediction，95%是可以通过sangar sequencing validated。当然在找出SNPs/Indels之后，除了reoccurrence以外（这个需要大量sample），如何去ranking这些 SNPs/Indels的重要性，例如处在那些基因，有多少是missense mutation，在一个pathway中有多少基因被mutated，mutation所处的位点是否对some sexy gene的expression有影响，如何linkage with expression data等等，尚处在研究阶段。因此做为marker，reoccurrence mutation的detection有重要意义，但同时其后续的functional study比较难，因为在找到的mutation之中，很难区分哪个是driver，哪个是receiver。除此之外，由于cancer tumor的 heterogeneity，要想发现那些真正在tumor cell里面的mutation，asks for new bioinformatics algorithms。
RNA-seq通常又被称为transcriptome sequencing。RNA-seq can be used to study
RNA expression variation and gene isoforms structure variation between
samples, novel RNA(such as long-intergenic-noncoding RNA)detection, microRNA
targets etc. Sequencing depth对于RNA-seq来讲非常重要，很多，例如novel
transcripts没有deep到一定程度是很难detect到的。现在比较流行的有Illumina
highseq 2000，和Roche 454，前者可以sequence up to 200M reads（pair-end）后者
号称可以detect longer reads up to 800bp（我没用过，不做评论）。从
bioiformatics的角度讲，与DNA-seq相比，RNA-seq要稍微复杂一些，当然也更有意思
。主要需要处理的就是alternative splicing。Tophat（ http://tophat.cbcb.umd.edu/）是目前比较流行的RNA-seq aligner。其部分也是基于bowtie，只对bowtie unalignable reads处理的时候，才去考虑splicing。后续对transciptome的 prediction，scripture(http://www.broadinstitute.org/software/scripture/)和 cufflinks(http://cufflinks.cbcb.umd.edu/manual.html)都是比较常用的reference based transcriptome assembler. 他们的主要作用就是把tophat align好的reads（当然也包括tophat predicted splicing junctions）组装成transcriptome。通过对已知gene annotation（例如RefSeq genes，UCSC genes）等的筛选，可以找到那些处于intergenic 并且有信号的transcripts，而通常这些是novel transcripts。这两个assembler共同的缺点是，在同一个genomic loci都会predict非常多的isoforms，故而FP比较高，但好处是，可以帮助reseracher非常迅速的找到novel transcripts可能处于的位点。除此之外，cufflinks还可以用来检测gene/transcripts expression difference(there are better ways to do this)，比较两个assembled transcriptome的不同，merge multiple transcriptome etc.就annotate transcriptome来讲，RNA-seq is much better than previously used tilling-array。
总体来说，对RNA-seq的数据分析尚未被标准化，相对于对microarry的分析，即便是最
基本的call significant genes in RNA-seq依然处于发展阶段，在microarry，大家都
知道，有RMA normalization，数据被normalize 之后，可以用ttest，SAM等手段找
significant genes，但是RNA-seq，单就normalization仍尚无定论。不同之处在于
microarry更类似于模拟信号，sequence 更像数字信号,比microarry 用的probe要更准
确，精度更高。此前业内标准是用RPKM（read per kilo-base of exon per million
reads）来代表gene expression level，其实是比较粗糙的，因为在RNA-seq 的
library中，有些RNA expression非常高，而有些去非常低，variation非常大，如果仅
仅用number of mappable reads来normalize，并不是非常好。相对而言，DESeq（R
package）和next-gen SAM algorithm是比较好的tool，好像他们都是用了quantile
normalization(就是除去两头，用中间的50%做normalization)。
RNA-seq也有被用来做structure variation，我尚未用过一个比较reliable的tool。
RNA-seq也可以用来检测在RNA pool里面是否有bacterial和viral infection，PathSeq
(http://www.broadinstitute.org/software/pathseq/index.html)给出了一个比较好的pipeline，但是他们的software并不是work非常well，至少不是我想象中那么便于操作，我们自己写了一个。
俗话说，眼见为实，信号的visulazation对很多ongoing project有非常重要的意义，
有的甚至可以改变问题的问法。我们lab应该说，都是UCSC genome browser的忠实用户
，理由：1.UCSC 提供很多已经做好的track，例如各种gene annotation，Encode/
Gencode里面TF和histome marker的ChIP-seq/RNA-seq信号等等，只要active就可以看
到；2.可以无限量upload custom tracks from your own data, and share the link
to whoever you want to. 在做自己的track的时候，通常会遇到很多file format
conversion，UCSC自己有一套convert tools（http://hgdownload.cse.ucsc.edu/admin/exe/macOSX.i386/）, BedTools(http://code.google.com/p/bedtools/)也是比较常用的。具体来讲，就是通过这些tool，把 align好的reads（usually in SAM or BAM format）变成UCSC tracks(UCSC has a number of accepted track file formats: http://genome.ucsc.edu/goldenPath/help/customTrack.html)。UCSC 还support HTTP and FTP for tracks, usually if a file is greater than 300M, we will put it on Amazon cloud and give the link to UCSC, and we will see the track，and it’s actually very fast.
综上，简单的介绍了一些据我所知的NGS前期数据分析方法，对于每一个project而言，
其问题的提法不同，所需要用到的后续工具也不同，很多情况下需要自己写script去解
决。随着实验技术和软件技术的发展，我相信越来越多的非常有趣的问题可以通过一些
整合的方法得以解决，或者可以通过整合的方法得到解决这些问题的hint。
I forgot to talk about de novo assembly, and thank you for making this point.
For de novo assembly, you can try SOAPdenovo" http://soap.genomics.org.cn/soapdenovo.html", BGI actually has a lot SOAP series, I've never tried and compared to the others, but will soon, I guess.
I did try trinity from Broad "http://trinityrnaseq.sourceforge.net/". very similar algorithms, but people from Broad definitely have more connections than BGI, and I guess that was the main reason why their story gets published on Nature. Biotech.(too far away from the topic)
Usually de novo assembly takes much more RAM than reference based assembly,
so it should work for small genomes, for large genome like human/mouse, what
people usually do is locate a loci, extract reads mapped to this loci, and
do local de novo assembly. This local de novo assembly was also used for SNP
detection in some algorithms.
标 题: 简单介绍 Bioinformatics Tools for NGS 分析
发信站: BBS 未名空间站 (Sat Sep 10 19:30:42 2011, 美东)
本人在wet lab里面做纯数据分析，for NGS data analysis, 简单介绍一些自己接触过
，并且觉得挺有用的工具，说的有点杂，权作抛砖引玉，还请不吝赐教。
Next-Gen sequencing(NGS)和现在正在发展的3rd-gen sequencing将会在生物学研究中
被越来越广泛应用。不管你信不信，反正我信了。一是基于实验成本的降低（$1k
whole-genome sequencing is coming），越来越多的实验室可以操作；二是可以提供
相对low throughput experiment多的多的数据和信息，可以看到很多从前看不到的东
西；三是sequencer本身对测序的准确性正在逐渐提高，所以实验固有错误率降低；四
是各种算法的成熟应用，这使得很多由于实验产生的误差在出数据后通过对数据的分析
得以过滤。按照library preparation来分，NGS主要有DNA-seq和RNA-seq
DNA-seq is usually used as ChIP-seq to study transcription factor(TF)-DNA
binding, epigenetic histome modification (such as H3K4me3, H3K36me3 etc),
DNA methylation etc. and Exon-seq or whole-DNA sequencing to study genetic
variation (SNPs/Indels etc). 从bioinformatics的角度来说，DNA-seq相对比较容易
处理，简单的说，第一步先做alignment，就是把millions of short reads from the
sequencer map back to the genome。很多genome已经有相对完整的reference，这些
genome的whole sequence可以从UCSC genome browser(http://genome.ucsc.edu/)下载。先前已经有朋友贴过现有的alignment tools（ http://en.wikipedia.org/wiki/List_of_sequence_alignment_software#Short- Read_Sequence_Alignment），目前比较常用的有bowtie，bwa，maq，soap，从他们各自的网站，应该都可以下载到 executable binary或者source code。我们比较常用的是bowtie（ http://bowtie-bio.sourceforge.net/index.shtml ）和bwa（ http://bio-bwa.sourceforge.net/），简单方便速度快（paralyzed）；output files例如SAM，BAM format比较standardized，方便后续处理；他们的网站有非常详细的manual。
如果是做ChIP-seq，通常第二步是call peaks。TF的peak一般比较peaky，MACS(http://liulab.dfci.harvard.edu/MACS /00README.html)是个不错的peak caller. TF ChIP-seq，一般会在call出peak之后assign peak to gene, in order to find TF regulated genes，这个比较arbitary，overlap of the called peak with gene promoter, enhancer, or use GREAT(http://great.stanford.edu/) to assign genes. 通常TF-DNA binding都有motif，meme（http://meme.sdsc.edu/meme/intro.html）是个不错的motif caller，use meme你可以看看这些TF enriched sites是否有significant motif, you can also check whether a known motif is enriched in a certain genomic region you interested.也就是说，既可以从genomic region找motif，也可以反向从motif找potential binding sites。
Histome modification，例如H3K4me3, K3K36me3,H3K27me3,H3K27ac, etc通常不具有
motif（not expected to），usually people are interested in their position,
downstream genes,how the intensity of peak changes (from average diagram),
whether there is bivalent domain (H3K4me3, H3K27me3) etc. 整体来说，对于ChIP
-seq,可以利用open source tool自己整workflow；可以利用galaxy(http://main.g2.bx.psu.edu/)已有的workflow 和file tools（上面也有详细的tutorial）；也可以用Partek Genomic Suite （not free,licence required）.后两个适合bench worker使用，非常容易上手。
Exon-seq and whole-DNA-seq asks very different question, usually people are
interested to know genetic variations, such as SNPs, Indels, Copy number
variation and looking for reoccurrence in tumor samples. 第一步alignment几乎
没有什么区别，找genetic variation, SeqGene(http://sourceforge.net/apps/mediawiki/seqgene /index.php?title=SeqGene)是一个比较简单且实用的tool。SNPs/Indels的prediction，基本上是比较准确 的，可以说那些high quality的prediction，95%是可以通过sangar sequencing validated。当然在找出SNPs/Indels之后，除了reoccurrence以外（这个需要大量sample），如何去ranking这些 SNPs/Indels的重要性，例如处在那些基因，有多少是missense mutation，在一个pathway中有多少基因被mutated，mutation所处的位点是否对some sexy gene的expression有影响，如何linkage with expression data等等，尚处在研究阶段。因此做为marker，reoccurrence mutation的detection有重要意义，但同时其后续的functional study比较难，因为在找到的mutation之中，很难区分哪个是driver，哪个是receiver。除此之外，由于cancer tumor的 heterogeneity，要想发现那些真正在tumor cell里面的mutation，asks for new bioinformatics algorithms。
RNA-seq通常又被称为transcriptome sequencing。RNA-seq can be used to study
RNA expression variation and gene isoforms structure variation between
samples, novel RNA(such as long-intergenic-noncoding RNA)detection, microRNA
targets etc. Sequencing depth对于RNA-seq来讲非常重要，很多，例如novel
transcripts没有deep到一定程度是很难detect到的。现在比较流行的有Illumina
highseq 2000，和Roche 454，前者可以sequence up to 200M reads（pair-end）后者
号称可以detect longer reads up to 800bp（我没用过，不做评论）。从
bioiformatics的角度讲，与DNA-seq相比，RNA-seq要稍微复杂一些，当然也更有意思
。主要需要处理的就是alternative splicing。Tophat（ http://tophat.cbcb.umd.edu/）是目前比较流行的RNA-seq aligner。其部分也是基于bowtie，只对bowtie unalignable reads处理的时候，才去考虑splicing。后续对transciptome的 prediction，scripture(http://www.broadinstitute.org/software/scripture/)和 cufflinks(http://cufflinks.cbcb.umd.edu/manual.html)都是比较常用的reference based transcriptome assembler. 他们的主要作用就是把tophat align好的reads（当然也包括tophat predicted splicing junctions）组装成transcriptome。通过对已知gene annotation（例如RefSeq genes，UCSC genes）等的筛选，可以找到那些处于intergenic 并且有信号的transcripts，而通常这些是novel transcripts。这两个assembler共同的缺点是，在同一个genomic loci都会predict非常多的isoforms，故而FP比较高，但好处是，可以帮助reseracher非常迅速的找到novel transcripts可能处于的位点。除此之外，cufflinks还可以用来检测gene/transcripts expression difference(there are better ways to do this)，比较两个assembled transcriptome的不同，merge multiple transcriptome etc.就annotate transcriptome来讲，RNA-seq is much better than previously used tilling-array。
总体来说，对RNA-seq的数据分析尚未被标准化，相对于对microarry的分析，即便是最
基本的call significant genes in RNA-seq依然处于发展阶段，在microarry，大家都
知道，有RMA normalization，数据被normalize 之后，可以用ttest，SAM等手段找
significant genes，但是RNA-seq，单就normalization仍尚无定论。不同之处在于
microarry更类似于模拟信号，sequence 更像数字信号,比microarry 用的probe要更准
确，精度更高。此前业内标准是用RPKM（read per kilo-base of exon per million
reads）来代表gene expression level，其实是比较粗糙的，因为在RNA-seq 的
library中，有些RNA expression非常高，而有些去非常低，variation非常大，如果仅
仅用number of mappable reads来normalize，并不是非常好。相对而言，DESeq（R
package）和next-gen SAM algorithm是比较好的tool，好像他们都是用了quantile
normalization(就是除去两头，用中间的50%做normalization)。
RNA-seq也有被用来做structure variation，我尚未用过一个比较reliable的tool。
RNA-seq也可以用来检测在RNA pool里面是否有bacterial和viral infection，PathSeq
(http://www.broadinstitute.org/software/pathseq/index.html)给出了一个比较好的pipeline，但是他们的software并不是work非常well，至少不是我想象中那么便于操作，我们自己写了一个。
俗话说，眼见为实，信号的visulazation对很多ongoing project有非常重要的意义，
有的甚至可以改变问题的问法。我们lab应该说，都是UCSC genome browser的忠实用户
，理由：1.UCSC 提供很多已经做好的track，例如各种gene annotation，Encode/
Gencode里面TF和histome marker的ChIP-seq/RNA-seq信号等等，只要active就可以看
到；2.可以无限量upload custom tracks from your own data, and share the link
to whoever you want to. 在做自己的track的时候，通常会遇到很多file format
conversion，UCSC自己有一套convert tools（http://hgdownload.cse.ucsc.edu/admin/exe/macOSX.i386/）, BedTools(http://code.google.com/p/bedtools/)也是比较常用的。具体来讲，就是通过这些tool，把 align好的reads（usually in SAM or BAM format）变成UCSC tracks(UCSC has a number of accepted track file formats: http://genome.ucsc.edu/goldenPath/help/customTrack.html)。UCSC 还support HTTP and FTP for tracks, usually if a file is greater than 300M, we will put it on Amazon cloud and give the link to UCSC, and we will see the track，and it’s actually very fast.
综上，简单的介绍了一些据我所知的NGS前期数据分析方法，对于每一个project而言，
其问题的提法不同，所需要用到的后续工具也不同，很多情况下需要自己写script去解
决。随着实验技术和软件技术的发展，我相信越来越多的非常有趣的问题可以通过一些
整合的方法得以解决，或者可以通过整合的方法得到解决这些问题的hint。
I forgot to talk about de novo assembly, and thank you for making this point.
For de novo assembly, you can try SOAPdenovo" http://soap.genomics.org.cn/soapdenovo.html", BGI actually has a lot SOAP series, I've never tried and compared to the others, but will soon, I guess.
I did try trinity from Broad "http://trinityrnaseq.sourceforge.net/". very similar algorithms, but people from Broad definitely have more connections than BGI, and I guess that was the main reason why their story gets published on Nature. Biotech.(too far away from the topic)
Usually de novo assembly takes much more RAM than reference based assembly,
so it should work for small genomes, for large genome like human/mouse, what
people usually do is locate a loci, extract reads mapped to this loci, and
do local de novo assembly. This local de novo assembly was also used for SNP
detection in some algorithms.
Subscribe to:
Posts (Atom)