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/


Monday, November 5, 2012

Wednesday, October 31, 2012

bio386, paper presentations

-----------
2012 Oct 31, Wed
Costanzo 2012 Science, genetic network -> Jasmine, 2012 Oct 31 PDF document
zotenko 2008 plos, essential hub genes -> Jessica, 2012 Oct 31 presentation

Monday, October 29, 2012

bio386, Dhami11, Gasch2000

------------
2012 Oct 29. Monday
Dhami 11 paper, Ella and Teneisha, 2012 Oct 29.

Gasch 2000 paper, Khayla and Asha, 2012 Oct 29
 Clustering exercise
 CV calculation

Saturday, September 1, 2012

bio386 fall 2012 log

------------
2012 Oct 29. Monday
Dhami 11 paper, Ella and Teneisha, 2012 Oct 29.

Gasch 2000 paper, Khayla and Asha, 2012 Oct 29
 Clustering exercise
 CV calculation

-------
2012 Oct 24 Wed.
inclass self-paced moodle quiz on qin's research proposal.
------------
2012 Oct 22 Monday
 review exam
Read DNA microarray wikipedia entry:
  http://en.wikipedia.org/wiki/DNA_microarray

read Qin's proposal
read Gasch paper, presentation
what is robustness? how to measure them?

--------------
2012 Oct 17 Wed, midterm exam, collaborative part
em: permutaiton on gIN

--------
3 days, intro to cellular aging and research
---------------------
2012 Oct 10.
 yeast aging presentation to m0 and G formula.

 take-home exam, due Oct 17 Wed before class. 1) permutation on genetic network 2) compare yeast gNet and PIN: are there correlation between genetic interactions and protein interactions?

 review genetic network assignment
 student presentations (Dhami paper, Newman paper)

--------
2012 Oct 8.
 intro to yeast aging.
    (aging, gompertz model, dS/dt = slope, non-aging)
    RLS, CLS
 figure 3B.

-----
2012 Oct 3, Wed
  match slides

  Figure 2
  partial regression demo, age, shoesize, readingability

  Permutation
  homework on figure 3B.
-----

quiz on Wed to see whether Fraser paper data can be loaded to R and first few lines can be ran.
set up dropbox accounts:

-----
2012 Sep 26 Wed
summary key points
  Ka, Ks, poistive, negative, neutral selection
  pairwise interaction -> network
  yeast ORF names
  competive growth fitness measures (using bar codes)
  use table to counter protein interactions
  data.frames

  match

For review and quiz, load genetic interaction data into R.
convert ORF to letters.

----------
2012 Sep 24, Mon
 discussion of Fraser Science paper

Some key ponts:

Protein interaction network
For Fig 1:
connectivity = number of interactions per gene
evolution rate =K
linear regression
For Fig 2: causal relationships (daycare example), multiple regression
For Fig 3: null distribution and p-value; permutation

Key messages:
 science paper is not a big deal
 how to read a scientific paper
 how to interpret figures
 p-value and null distribution

----------
2012 Sep 19, Wed
 ***collaborative exam, student demo and exercise on screen. This method works well!
-----------
2012 Sep 17, Mon
 exercises and reviews
 go over exam1
-------------
2012 Sep 12 Wed
  2.exericse
  function on make solution
  give home work
  input out,
 summary, leave homework, quiz on Monday.

---------------
=>2012 Sep 10, Mon
  quiz on daycare score again
  basic programming concepts: loops, conditions, functions.
---------------
=>2012 Sep 5, Wed
  Quiz on salary.R
  Homework on fraser paper
  salary.R 
  review new learnings. regression, t.test(), p-value.
  homework is requested.
//
-------------------
=>2012 August 29, Wed
  quiz: a list of values (0.1, 0.01, 0.001), take log
  irb signature: ella
  salary example. socratice method, working on it, the ask questions.
  give review guide.

skills and functions to cover:
str, read.table(), pick columns, pick rows, pick rows and columns,
hist
table() #how many female ‘Arts’ faculty
#how many female ‘Assistant’, ‘Associate’, and ‘Full’ professors
//

-------------------
=>2012 August 27, Mon
  welcome new comers
  simple.R takes 2 hours
//

-------------------
=> 2012 Aug 22, Wed
x irb signature
x go over syllabus, writing milestones, previous grades,
x pre-survey
* foucs on why computing for biology majors (four students just come to class for no reasons without computing backgrounds)
 googleDoc to work on list of reasons for biology majors to learn programming and computing concepts and skills
using cards to explain computational thinking and search space

Reasons why biology major should learn programming and computing
    Learning programming skills and computational thinking offers a key set of problem solving skills.
    genomics medicine
    highlights on resumes when applying for medical school or graduate schools
    more competitive for jobs
    Another useful language
    Computing holds the key to understand complicated biology phenomena.
    Further understanding of genomic dynamics of emerging pathogens and diseases
    To use computer storage and other programs to study the biological data and maps of the human genome sequence
    Updated on new technology
    Understanding programming and computing can aid in drug discovery
    Aids in conducting research projects

How to find out a missing card from a stack? (A key concept is search space)
Plan A: count card one by one. Search space is 52 possibilities.
Plan B: separate into sets and colors
Question: What are the key differences between these two search plans?
Answers: Random data versus organized (structured) data

* I should ask students to introduce themselves.

Tuesday, August 21, 2012

20120821 DHE DHR BY4743 H2O2

During Clibur measure, we measured blanks, DHE, DHR, and DHEDHR, and "on" for overnight DHR signals. The overnight DHR sample were basically left in the dark at room temperature (as I recall).





Wednesday, February 22, 2012

Tuesday, January 10, 2012

bio386, spring 2012 log

-> exam question: yeast gene in FASTA format, calculate its length, find out
 mean ORF length.



partion H2O2-LOH data to bio320 classes for assignment

Exam 2-> deg factor on salary


-> bayesian network modeling?
-> bio320 projects: story paper, morphology, expression; wynter gwa study; realibility simulation and threshold trait;


=> 2012 Feb 21
 Dropbox account creation for projects
 Levy's presentation and explanation of data;
 Go over research projects;
 Pubmed and Google scholar search to find the most recent protein interaction network data in S. cerevisiae.


=>2012 Jan 26, 3.zip. The exercises take the entire class,

=>2012 Jan 23Tue discuss my MRI proposal

=>2012 Jan 19Thu
 proofread an NSF proposal. 5 students, working 3 groups. 5 pages per group.
//

=> Jan 17, 2012 Tues
work through salary example. R2 and p-value take a long time.
One students ask for open-ended qustions. I asked them study
'deg' effect on 'salary', 'gender', and 'rank'.


=> Jan 12, 2012 Thursday
survey
hp laptops usages
Moodle site
go over syllabus, writing milestones
previous student projects, and posters
install R, Rstudio, notepad++,
worked on simple.R
assignment: watch youtube videos, 1st ppt, salary.R

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.