Difference between revisions of "Tips kang"

From Crop Genomics Lab.
Jump to: navigation, search
Line 4: Line 4:
  
 
'''soap3-dp'''
 
'''soap3-dp'''
  /data/k821209/programs/soap3-dp/4.2-amazon/soap3-dp-builder adzuki.ref.fasta
+
  /data/k821209/programs/soap3-dp/4.2-amazon/soap3-dp-builder adzuki.ref.fasta # reference fasta file 에 컨티그나 스카폴드가 너무 많으면 에러남.. 몇개가 적정인지는 모르겠음
 
  /data/k821209/programs/soap3-dp/4.2-amazon/BGS-Build adzuki.ref.fasta.index
 
  /data/k821209/programs/soap3-dp/4.2-amazon/BGS-Build adzuki.ref.fasta.index
 
  /data/k821209/programs/soap3-dp/4.2-amazon/soap3-dp pair adzuki.ref.fasta.index test_1.fastq test_2.fastq -u 500 -v 200 -b 3
 
  /data/k821209/programs/soap3-dp/4.2-amazon/soap3-dp pair adzuki.ref.fasta.index test_1.fastq test_2.fastq -u 500 -v 200 -b 3

Revision as of 15:37, 19 July 2014

GBS

parallel -j 7 bcftools index {1} ::: *.bcf
parallel -j 7 --colsep='\t' python3 genotyping.py {1} {2} :::: sites.txt

soap3-dp

/data/k821209/programs/soap3-dp/4.2-amazon/soap3-dp-builder adzuki.ref.fasta # reference fasta file 에 컨티그나 스카폴드가 너무 많으면 에러남.. 몇개가 적정인지는 모르겠음
/data/k821209/programs/soap3-dp/4.2-amazon/BGS-Build adzuki.ref.fasta.index
/data/k821209/programs/soap3-dp/4.2-amazon/soap3-dp pair adzuki.ref.fasta.index test_1.fastq test_2.fastq -u 500 -v 200 -b 3
/data/k821209/programs/samtools-0.1.19/samtools cat -o V_nepal_to_redbean.bam V_nepal_1_1.fastq.gz.gout.1 V_nepal_1_1.fastq.gz.gout.2 V_nepal_1_1.fastq.gz.gout.3 V_nepal_1_1.fastq.gz.gout.4 V_nepal_1_1.fastq.gz.dpout.1


gene prediction length distribution

python3 0_length_dist.py Gmax_189_gene.gff3
paste -d , Gmax_189_gene.gff3_CDS.ld.out Gmax_189_gene.gff3_gene.ld.out Gmax_189_gene.gff3_intron.ld.out > gm.paste
python3 add_spcs.py gm.paste gm
모든 species에 대해서 같은 작업을 한뒤
*.Rin 중 하나만 첫번째 column을 Species로 이름을 정해주고 나머지들은 header 들을 지워준다.
cat gm.Rin 등등 > merged.Rin # 해더가 있는 애들 처음으로 둔다.
R
library(ggplot2)
Data <- read.csv("merged.Rin",header=T)
> ggplot(Data, aes(x=CDS,colour=Species, fill = Species)) + geom_bar(alpha = 0.3,binwidth=50,position="identity")+ xlim(0,2000)                                                
> ggplot(Data, aes(x=CDS,colour=Species, fill = Species)) + geom_bar(alpha = 0.2,binwidth=50,position="identity")+ xlim(0,2000)
> ggplot(Data, aes(x=intron,colour=Species, fill = Species)) + geom_bar(alpha = 0.2,binwidth=50,position="identity")+ xlim(0,2000)
> ggplot(Data, aes(x=gene,colour=Species, fill = Species)) + geom_bar(alpha = 0.2,binwidth=500,position="identity")+ xlim(0,20000) 
Density plot 가능
> ggplot(Data, aes(x=gene,colour=Species)) + geom_density(alpha = 0.4,position="identity")+ scale_x_continuous(limits = c(0, 2000)) 
문제는 NA가 있는 row를 다 제거한다는 문제가 있다.


qplot

http://ggplot2.org/book/qplot.pdf

Jatropha transcriptome

244:/home/k821209/py/ # used code
transdecoder runMe.sh
244:$ cdhit -i Trinity.fasta.transdecoder.pep -o Trinity.fasta.transdecoder.pep.cdhit.fa
244:$ python3 pep_header_parser.py Trinity.fasta.transdecoder.pep.cdhit.fa JPO # JPO : 3글자 약어를 각 종 마다 지정해줌

div_WGD_gen.Ks.py

input : filename.colinearity.kaks # MCSCANX result
output : filename.colinearity.kaks.Rin
R
library(ggplot2)
Data <- read.table("Va2Gm.collinearity.kaks.Rin",header=F)
qplot(V1,data=Data,binwidth=0.01,xlim=c(0,3),ylim=c(0,150),xlab="Ks bin",ylab="Frequency")


CDHIT

est
63:$ /data/program/cd-hit-v4.6.1-2012-08-27/cd-hit-est -i merge.est.fa -o merge.est.fa.cdhit.fa -T 10 -M 10000

Phyml

193:$ mpirun -n 16 phyml-mpi -i concate.fa.phy -q -d aa -m JTT -b 500

Interpro

193:$ /data2/k821209/programs/interproscan-5.4-47.0/interproscan.sh --pathways --goterms -iprlookup -T /data2/k821209/Mungbean/maker4/TMP/ -i /data2/k821209/Redbean/maker/redbean.ref.maker.output/redbean.ref.fasta.pep.fa.noasterix.fa

jellyfish

$ sudo bash
$ unlimit -a 4096

파일을 동시에 많이 쓰기 때문에 열어줘야된다.

$ zcat ysp-2_?.fastq.gz | jellyfish count -C -m 22 -s 50G -t 10 -c 6 -o output /dev/fd/0
$ jellyfish merge ./output_* # 결과 파일이 많을 경우 
$ jellyfish histo -o hist.out output_0 
$ R
library(ggplot2)
library(quantmod)
Data<-read.table("hist.out",header=F)
p <- qplot(V1,V2,data=Data,xlim=c(3,200),ylim=c(1,20000000),xlab="Kmer frequency",ylab="Count")
p + annotate("text", label = sprintf("Peak : %s",findPeaks(Data$V2[2:200])[1]) , x = findPeaks(Data$V2[2:200])[1], y = Data$V2[findPeaks(Data$V2[2:200])[1]]+1000000, size = 8, colour = "red")

Species divergence time with fixed substitution rate

63:/data2/Mungbean_revision/speciation/Orthomcl/OnlyVigna/pairs

작업경로

63:/home/k821209/py/pairwiseKs/

저장경로

Orthomcl을 통해 얻은 Ortholog.txt가 필요
Orthomcl의 orthomclAdjustFasta로 만들어진 cds sequence set이 필요
$ python3 shared_ortholog.py
$ python3 get_table.py

BEASTv1.8

http://www.molecularevolution.org/resources/activities/beast_activity/viruses

tutorial

193:/data2/k821209/good_python_code/BEAST/

쓰이는 코드가 저장되어있는 경로

PREPARE
- MCSCANX -
Vr2Cc.collinearity.kaks    : V.radiata vs C. cajan
Vr2Gm.collinearity.kaks    : V.radiata vs G. max
Vr2Wd.collinearity.kaks    : V.radiata vs V.radiata var. sublobata (wild mungbean denovo assembly)
Vra2Vrr.collinearity.kaks  : V.radiata vs V.reflexo-pilosa var. glabra (allotetraploid denovo assembly)
GM2GM.collinearity.kaks    : G.max within genome synteny block
Vrg2Vrg.collinearity.kaks  : V.reflexo-pilosa var. glabra within genome synteny
- ORTHOMCL - 
$ python3 div_WGD_gen.Ks.py filename.kaks 

본 스크립트를 실행하면 0.01 Ks window로 Frequency를 눈으로 볼 수 있는데 peak 가 보이는 영역 중 Frequency가 10 이상인 영역을 확인하고 코드 안에서 targetMin = 0.04, targetMax = 0.12 값을 수정해주면 filename.kaks.target이라는 파일을 생성한다. 해당 영역의 평균, 최대, 최소 MYA를 계산한다. 고정값 상수를 사용하므로 염두에 둔다. 그리고 생성되는 filename.kaks.recentpeakWGD 에는 설정한 범위에 해당하는 synteny pair들이 저장된다. Within genome 의 경우는 A genome, B genome을 구분할 기본 데이터가 되고, Inter genome 의 경우에는 high confident ortholog가 된다.

file_fa = 'final.assembly.11.fasta.nonATGC.fa.cds.primary.fasta.pep.sh.fa'
file_vra2vrp = 'Vra2Vrr.collinearity.kaks.recentpeakWGD'
file_vrp2vrp = 'Vrg2Vrg.collinearity.kaks.recentpeakWGD'
$ python3 sep.py

Within genome synteny에서 A genome B genome 을 구분하는 코드. Within genome의 Recent peak을 A genome과 B genome관계라고 가정하고 V radiata와 가까운 Synteny pair를 A genome, 먼쪽을 B genome으로 저장한다. match.txt에는 A genome과 B genome내 유전자 matching 정보가 들어간다.

mcscanout_list = ['Vr2Cc.collinearity.kaks.recentpeakWGD','Vr2Gm.collinearity.kaks.recentpeakWGD','Vr2Wd.collinearity.kaks.recentpeakWGD','Vra2Vrr.collinearity.kaks.recentpeakWGD']
file_GMAA       = 'GMAA.fasta'
file_VRGA       = 'VRPA.fasta'
file_matchGM    = 'match.GM.txt'
file_matchVRG   = 'match.Vrg.txt'
file_fa         = 'all.fas'
$ python3 0_get_SYNB.py

V.radiata, C.cajan, G.max, V. radiata var. sublobata, V. reflexo-pilosa 간의 공유된 Synteny block을 GMAX A genome과 V.replexo A genome을 기준으로 가져온다. 본 연구에서 173개의 Locus가 나옴.

fa_list = glob.glob('*.prealn.fa')
orthologs = 'orthologs.txt'
dicAHD2seq = kang.Fasta2dic('all.fa')
$ python3 0_ortholog_add.py

ortholog.txt는 Vigna species만 사용해서 만든 orthomcl pair파일. *.prealn.fa 는 0_get_SYNB.py 의 결과파일. all.fa는 필요한 시퀀스들.. orthomcl의 compliantFasta/ 를 모으면 됨.

$ parallel -j 10 /data2/k821209/programs/prank/bin/prank -d={1} -o={1}.aln ::: *.good.fa

결과 파일 Alignment

$ python3 3_fastaaln2nexus.py
$ parallel -j 10 java -jar /data2/k821209/programs/prottest-3.4-20140123/prottest-3.4.jar -i {1} -S 0 -F -AIC -BIC -tc 0.5 -o {1}.prottest  ::: *.good.fa.aln.best.fas.nex

nexus format으로 변환. 꼭 안바꿔도 되는듯(?)

$ cat *.good.fa | grep '>' > prespecies
$ python3 species.py > species

beauti에 들어가는 species 파일 생성

BEAST 는 두번돌림. 50,000,000 change of length로 두번돌림

/data2/k821209/programs/BEASTv1.8.0/bin/logcombiner -trees -burnin 0 StarBEASTLog.species.trees StarBEASTLog.species.second.trees merged.trees

두번돌리고 Merge


Maker

$ /data2/k821209/programs/maker/bin/gff3_merge -d redbean.ref_master_datastore_index.log
$ cat redbean.ref.all.gff | awk '$2=="maker"' - > redbean.ref.all.gff.maker.gff3
$ python3 sort.py redbean.ref.all.gff.maker.gff3
$ python3 genename_change_nonAnchor.py
$ /data2/k821209/programs/maker/bin/map_gff_ids Vang.scaffold.map redbean.ref.all.gff
$ cat redbean.ref.all.gff | awk '$2=="maker"' - > redbean.ref.all.gff.final.maker.gff3

Orthomcl

orthomcl 는 알아서 이름을 바꾸면된다.
Make config file : orthomcl.config.template
dbVendor=mysql
dbConnectString=dbi:mysql:orthomcl:mysql_local_infile=1:localhost:3306
dbLogin=orthomcl
dbPassword=8804555
similarSequencesTable=SimilarSequences
orthologTable=Ortholog
inParalogTable=InParalog
coOrthologTable=CoOrtholog
interTaxonMatchView=InterTaxonMatch
percentMatchCutoff=50
evalueExponentCutoff=-5
oracleIndexTblSpc=NONE
$ mysql -u root -p
mysql> CREATE DATABASE orthomcl;
mysql> GRANT SELECT,INSERT,UPDATE,DELETE,CREATE VIEW,CREATE, INDEX, DROP on orthomcl.* TO orthomcl@localhost;
mysql> set password for orthomcl@localhost = password('8804555');
$ orthomclInstallSchema orthomcl.config.template
$ orthomclAdjustFasta VRA final.assembly.pep.fasta.primary.SH.fa 1
$ mkdir ./compliantFasta/
$ mv VRA.fasta ./compliantFasta/ # same for the ATH.fasta etc.
$ orthomclFilterFasta ./compliantFasta/ 10 20
$ blastall -m 8 -a 10 -p blastp -F 'm S' -v 100000 -b 100000 -z 259701 -e 1e-5 -d goodProteins.fasta -i goodProteins.fasta -o goodProteins.fasta.self.out # 27개 할 때 보름 이상걸림
$ parallel -j 10 "blastall -m 8 -p blastp -F 'm S' -v 100000 -b 100000 -z 259701 -e 1e-5 -d goodProteins.fasta -i {1} -o {1}.out" ::: ???.fasta
$ orthomclBlastParser goodProteins.fasta.self.out compliantFasta/ >> similarSequences.txt # 27개 accession할 때 이 스텝만 6시간 정도 걸림
$ orthomclLoadBlast orthomcl.config.template similarSequences.txt # 용량 많이 먹음 충분히 확보
$ orthomclPairs orthomcl.config.template orthomcl_pairs.log cleanup=no
$ orthomclDumpPairsFiles orthomcl.config.template
$ mcl mclInput --abc -I 1.5 -o mclOutput
$ orthomclMclToGroups CLST 1000 < mclOutput > groups.txt

R
source("http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/overLapper.R")
sets <- read.delim("groups.txt.Rin")
setlistImp <- lapply(colnames(sets), function(x) as.character(sets[sets[,x]!="", x]))
names(setlistImp) <- colnames(sets)
setlist5 <- setlistImp[1:5]; OLlist5 <- overLapper(setlist=setlist5, sep="_", type="vennsets")
counts <- sapply(OLlist5$Venn_List, length)
vennPlot(counts=counts, ccol=c(rep(1,30),2), lcex=1.5, ccex=c(rep(1.5,5), rep(0.6,25),1.5))


Errors
execute failed: Can't create/write to file /tmp 어쩌고..
 이것은 mysql 버전에러라고 하는군. my.cnf 에서 tmp 디렉토리를 data 디렉토리랑 이름을 같게 만들어주면 돌아감. 
Duplicate entry error
 $ sort -k3,4 similarSequences.txt.primary > similarSequences.txt.primary.sort
 $ python3 duplicates_remove_ver1.py # 63:/data2/k821209/good_python_code/ 
 현재 테스트중

Protest

193server$ cd /data2/k821209/programs/prottest-3.4-20140123/  # 이 디렉토리로 이동해서 실행. 라이브러리 위치 때문에 따로 지정하기 귀찮기 때문에.. 
193server$ parallel -j 10 java -jar prottest-3.4.jar -i {1} -o {1}.protest.out -all-distributions -F -AIC -BIC -tc 0.5  :::  /data2/k821209/Mungbean/Mungbean_assemble/ver6/mcscanx/beast_prepare/all/others_ver2/manual/*.nex

jModeltest for DNA alignment

63server $ parallel -j 10 java -jar /data/program/jmodeltest-2.1.5/jModelTest.jar -d {1} -g 4 -i -f -AIC -BIC -a -o {1}.jmodel ::: *.nuc.fas

bamUtil

bam2fastq

bam bam2FastQ --in filename.bam

Picard

sam2fastq

java -jar SamToFastq.jar INPUT=[file_bam] INCLUDE_NON_PRIMARY_ALIGNMENTS=True FASTQ=Name_1.txt SECOND_END_FASTQ=Name_2.txt VALIDATION_STRINGENCY=LENIENT

BEAST

beauti

/data2/k821209/programs/BEAST/bin/beauti -template /data2/k821209/programs/BEAST/templates/StarBeast.xml

beast

'single locus'

java -Xms256m -Xmx202400m -Djava.library.path=/data/program/BEAST/lib/ -jar /data/program/BEAST/lib/beast.jar -beagle_SSE -beagle_instances 10 ver7.xml

'multi loci'

java -Xms256m -Xmx202400m -Djava.library.path=/data/program/BEAST/lib/ -jar /data/program/BEAST/lib/beast.jar -beagle_SSE ver7.xml
java -Xms256m -Xmx202400m -Djava.library.path=/data/program/BEAST/lib/ -jar /data/program/BEAST/lib/beast.jar -beagle_SSE 557.xml # 557 loci : Total calculation time: 373863.472 seconds 

GPU with multi loci

C:\Users\jun\Desktop\kangyangjae\BEAST.v2.1.2\BEAST\557>java -Xmx15g -jar ..\lib\beast.jar -beagle_order 0,0,0,0,0,0,0,1

0,0,0,0,0,0,0,1 의 의미는 "7개 loci정보 cpu에 던지고 8번째는 gpu에 던져라" loci가 많을 경우는 반복으로 할당. gpu에 다 던지면 좋겠지만 gpu메모리 한계 때문에 이렇게 해야함.

18개 loci, 27 species alignment 에 GPU memory 770 Mb 정도 소요.

NCBI에 WGS 올리기 팁

꼼수이다. 참고로. 일단 그냥 Whole genome seqeunce fasta file을 그냥 WGS submission portal로 넣어버리고 하루를 기다린다. 뭐가 잘못되었다고 구구절절 온다. 이렇게 하라고도 구구절절온다. 그렇게 해서 보내준다. 끝.

https://submit.ncbi.nlm.nih.gov/subs/wgs/

Mutation rate

Legumes

"If the older duplication is assumed to have occurred around 58Myr ago, then the calculated rate of silent mutations extending back to the duplication would be 5.17*10-3, similar to previous estimates of 5.2*10–3" [1]

Python

Fisher's exact test

from scipy import stats

oddsratio, pvalue = stats.fisher_exact([[A,B], [C, D]]) [2]

63:/home/k821209/py/NGS/vcfq2fa.py : vcfutil로 만들어진 fq 파일을 fa로 변환

Excel

=TEXT(2.2323,"(0.00)")

(2.23)

Softwares

SNP 반영된 fasta 만들기

  1. java -jar /data/program/picard-tools-1.91/CreateSequenceDictionary.jar R=Gmax_189.fa O=Gmax_189.dict
    • 인덱싱 작업인듯
  2. java -Xmx2g -jar /data/program/GenomeAnalysisTK-3.1-1/GenomeAnalysisTK.jar -R /data/ref/Gmax_189.fa -T FastaAlternateReferenceMaker -o JM14.bam.sort.bam.bcg.bcf.d2D50.vcf.fa --variant JM14.bam.sort.bam.bcg.bcf.d2D50.vcf
    • bcftools view [filename.bcf] | /data1/KimSue/Gmax_189/ver1/bcf/bcf2fa.py prefix # GATK의 FastaAlternateReferenceMaker 는 calling이 되지 않는 지역을 reference 서열로 가져오는 문제가 있음.

GATK pipe http://www.broadinstitute.org/gatk//events/2038/

  1. bwa mem -M -t 10 Va.ref.fa ysp-2_1.fastq.gz ysp-2_2.fastq.gz | /data/program/samtools-0.1.19/samtools view -Sb - | /data/program/samtools-0.1.19/samtools sort - ysp.bwamem.Va.ref.fa.sort # GATK pipe는 -M 옵션이 필요
  2. /data/program/jdk1.7.0_25/bin/java -jar /data/program/picard-tools-1.91/MarkDuplicates.jar INPUT=ysp.bwamem.Va.ref.fa.sort.bam OUTPUT=ysp.bwamem.Va.ref.fa.sort.bam.dedup.bam METRICS_FILE=metrics.txt MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000 TMP_DIR=./tmp/
    • MarkDuplicates 는 PCR duplicates를 제거하기 위한 작업. 왜 제거해야하는지는 자세히는 모르지만 PCR bias 때문에 Variant calling statistical model 에 문제가 된다고 함. 특히 PCR중에 polymerase가 실수 했을 경우 실수한 자리가 PCR bias로 대량 생산되었을 경우 마치 SNP처럼 calling되버림. 제거해야한다는 것이 결론. 그러나 특별히 제작된 라이브러리에는 적용해서는 안된다. 예를들면 특정 사이트가 잘리게 만들어놓은 라이브러리? 류들은 PCR bias로 오해받을 가능성? 이 있는듯..
    • MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000 값을 높여도 되는데 ulimit -n 값을 올려야함. 재부팅 해야된다해서 그냥 낮춰 쓰는중. 그렇게 느린줄 모르겠음. 높여서 안써봐서 그런가..
    • TMP_DIR=./tmp/ 이거 빼먹으면 용량 없음 크리를 먹는다.

Maker

  1. /data2/k821209/programs/maker/bin/gff3_merge -d Va.ref_master_datastore_index.log
  2. /data2/k821209/programs/maker/bin/maker_map_ids --prefix=Vang --iterate=1 --suffix=. Va.ref.all.gff > id_map.txt # 복잡하게 나오는 maker의 유전자 이름들을 심플하게 바꾸는 툴
    • cat final.assembly_40k_GSFLX_pseudo5k.fasta.nonATGC.all.gff | awk '$2=="maker"' - > final.assembly_40k_GSFLX_pseudo5k.fasta.nonATGC.all.gff.maker.gff
    • python3 gff_parse.py final.assembly_40k_GSFLX_pseudo5k.fasta.nonATGC.all.gff.maker.gff
    • python3 /data2/k821209/Redbean/maker_pseudo/Va.ref.maker.output/genename_change_nonAnchor.py # 이름이 맘에 안들게 바뀌어서 개인적으로 만든 툴
  3. /data2/k821209/programs/maker/bin/map_gff_ids id_map.txt Va.ref.all.gff # 그리 나온 이름들을 gff 반영하는 툴
  4. python3 /data2/k821209/Redbean/maker_pseudo/Va.ref.maker.output/header_change.py Va.ref.all.maker.proteins.fasta Vang.scaffold.map # 그리 나온 이름들을 fasta에 반영하는 툴
  5. /data2/k821209/programs/maker/bin/iprscan2gff3 Va.ref.all.maker.proteins.fasta.tsv.hc.tsv Va.ref.all.gff > Va.ref.all.gff.ipr.gff # interpro result를 jbrowser에 들어가는 모양으로 만들어주는 툴

Deconseq [3]

  1. Illumina read의 contamination을 확인한다.
  2. 63:/data/program/deconseq-standalone-0.4.3
  3. /usr/bin/perl deconseq.pl -keep_tmp_files -f 800_both.fq -dbs bact,vir,arch -dbs_retain gmax

ePCR

Re-PCR

  1. $ famap -tN -b genome.famap ./superscaf.ver1.fa
  2. $ fahash -b genome.hash -w 12 -f3 ${PWD}/genome.famap
  3. Work> /data/program/e-PCR-2.3.12/re-PCR -S genome.hash -n1 -g1 SSR.sts -o SSR.sts.mapped

SSR.sts

Mungbean_SSR_ID_1 CAAAAACATGAGTTGCACACAA TCATAACGCAGAACAGCGAA

Mungbean_SSR_ID_2 ATGTGTGTGAGCACCTCGAC TTTGGCCATGCAAGATGTAA

Mungbean_SSR_ID_4 GCGGTTCACCTAGCCATAAA GGACCCTTCTGTGCGTGTAT

Mungbean_SSR_ID_5 GTTTGTGCTGCGGATTCTTT TTGGCAATTTGGACTAAGGC

Mungbean_SSR_ID_7 TTGACCCAAAACTTACCAATTT GCTAAGGACTGGGGGTCTTC

Mummer, alignment draft genome to finished genome

$nucmer --prefix=ref_qry ref.fasta qry.fasta

$show-coords -rcl ref_qry.delta > ref_qry.coords

$show-aligns ref_qry.delta refname qryname > ref_qry.aligns

$show-tiling ref_qry.delta > ref_qry.tiling

추천 논문

Paterson, A.H., M. Freeling, H. Tang and X. Wang. 2010. Insights from the comparison of plant genome sequences. Annual Review of Plant Biology 61: 349-372. [4]

  1. doi:10.1038/nature08670
  2. scipy, fisher's exact
  3. deconseq
  4. Paterson, A.H., M. Freeling, H. Tang and X. Wang. 2010. Insights from the comparison of plant genome sequences. Annual Review of Plant Biology 61: 349-372.