dreamfishes / rgaat Goto Github PK
View Code? Open in Web Editor NEWThis project forked from wushyer/rgaat
Reference based genome assembly and annotation tool for new genome and known genome upgrade
This project forked from wushyer/rgaat
Reference based genome assembly and annotation tool for new genome and known genome upgrade
NAME ==== RGAAT - Reference based genome assembly and annotation tool for new genome and known genome upgrade VERSION ======= Version 1.0 COPYRIGHT ========= Copyright 2016 - Wanfei Liu & Shuangyang Wu & Qiang Lin - Beijing Institute of Genomics, Chinese Academy of Sciences. This tool is a free package; you can redistribute it and/or modify it freely. INTRODUCTION ============ This program can assemble and/or annotate genome for new genome and known genome upgrade using sequence alignment file (SAM or BAM format), sequence variant file (VCF format or five coloum table (tab-delimited, including chromosome, position, id, reference allele and alternative allele)) or new genome sequence file (FASTA format) based on reference genome sequence file (FASTA format) and annotation file (TBL, GTF, GFF, GFF3 or BED format). Genome sequence, annotation and sequence variant files should be uncompressed files. The name of annotation file should be ended with suffix *.tbl, *.gtf, *.gff, *.gff2, *.gff3 or *.bed. If you want run program multiple times with different parameters, please use "-o" option to assign different prefix name for output files. To filter false positive sequence variants, we use 3 as the default minimum reads depth and 3 as the default minimum alt allele reads coverage, because most of false positive variants come from the low depth and low alt coverage records according to our previous study. Although we filtered genome comparative result (obtained by BLAT) for reference and new genome sequence, user can manually check the *.psl file for higher accuracy if possible and only remain the most possible direction results (plus or minus). Moreover, if comparing with different assembly version or different strain of same species, we can use the default minIdentity parameter (90) for BLAT, however, if comparing with different species, we recommend using 50 as minIdentity by option -m (-m 50) for BLAT. To convert BAM file to SAM file, RGAAT need samtools. RGAAT is freely available at github. PREREQUISITES ============= Before running, the following software or program are required: - Perl Environment (perl v5.18.4 or later (it can be lower)) - Perl module: Getopt::Long; Text::NSP; SVG - BLAT (BLAT v. 35 or above) - samtools (tested on samtools Version: 1.2) RGAAT has been tested on Linux. INSTALLATION ============ This package doesn't need installation. You can run it using absolute path after decompressed it. COMMAND LINE ============ Run RGAAT.pl without any parameter, it will print below usage on the screen. Author: Wanfei Liu Email: <[email protected]> Date: May 16, 2016 Version: 1.0 version Introduction This program can assemble and/or annotate genome for new genome and known genome upgrade using sequence alignment file (SAM or BAM format), sequence variant file (VCF format or five coloum table (tab-delimited, including chromosome, position, id, reference allele and alternative allele)) or new genome sequence file (FASTA format) based on reference genome sequence file (FASTA format) and/or annotation file (TBL, GTF, GFF, GFF3 or BED format). Usage: perl RGAAT.pl -g genome_sequence(FASTA) -a annotation(TBL, GTF, GFF, GFF3 or BED) [-s SAM_file | -b BAM_file | -v VCF_file | -n new_genome_file] -o prefix_of_output_file. -g: the absolute path of genome sequence file (FASTA format, required). -a: the absolute path of genome annotation file (TBL, GTF, GFF, GFF3 or BED format, required for annotation). -s: the absolute path of SAM file (sorted according to coordinate, it is mututally exclusive with -b, -v and -n). -b: the absolute path of BAM file (sorted according to coordinate, it is mututally exclusive with -s, -v and -n). -v: the absolute path of sequence variant file (VCF format or five coloum table (tab-delimited, including chromosome, position, id, reference allele and alternative allele), it is mututally exclusive with -s, -b and -n). -n: the absolute path of new genome sequence file (FASTA format, it is mututally exclusive with -s, -b and -v). -m: the minIdentity for BLAT comparative between reference and new genome (default value is 90). -q: quality threshold for reads (default value is 30). -l: read length threshold (default value is 30). -d: read depth threshold for sequence variant (default value is 3). -c: allele depth threshold for sequence variant (default value is 3). -p: allele proportion threshold for sequence variant (default value is 0.5). -fu: filter unpaired reads (yes or no, default value is no). -fm: filter multiple mapping reads (yes or no, default value is yes). -o: the prefix of output file (required). Note: Genome sequence, annotation and sequence variant files should be uncompressed files. The name of annotation file should be ended with suffix *.tbl, *.gtf, *.gff, *.gff2, *.gff3 or *.bed. If you want run program multiple times with different parameters, please use "-o" option to assign different prefix name for output files. To filter false positive sequence variants, we use 3 as the default minimum reads depth and 3 as the default minimum alt allele reads coverage, because most of false positive variants come from the low depth and low alt coverage records according to our previous study. Although we filtered genome comparative result (obtained by BLAT) for reference and new genome sequence, user can manually check the *.psl file for higher accuracy if possible and only remain the most possible direction results (plus or minus). Moreover, if comparing with different assembly version or different strain of same species, we can use the default minIdentity parameter (90) for BLAT, however, if comparing with different species, we recommend using 50 as minIdentity by option -m (-m 50) for BLAT. To convert BAM file to SAM file, RGAAT need samtools. INPUTFILES ========== RGAAT need at least two necessary input files: - Reference genome file in FASTA format. - Sequence alignment file (SAM or BAM format), sequence variant file (VCF format or five coloum table) or new genome sequence file (FASTA format). If you need to annoate upgrade or new genome, you must provide reference annotation file (TBL, GTF, GFF, GFF3 or BED format). We prefer to use TBL format file (feature table file, http://www.ncbi.nlm.nih.gov/projects/Sequin/table.html) as annotation file, which can be download from each genome sequence record by click “Send -> Complete Record -> File -> Format -> Feature Table -> Create File” easily in nucleotide database of NCBI. OUTPUTFILES =========== RGAAT mainly produces output files in simple tab-delimited table files. 1. *.var (variant calling result) The identified sequence variants in reference genome. ##fileformat=VCF like format ##INFO=<ID=RDP,Number=1,Type=Integer,Description=\"All mapped reads depth after reads filter\"> ##INFO=<ID=DP,Number=1,Type=Integer,Description=\"High quality reads depth; some reads and loci may have been filtered\"> ##INFO=<ID=AF,Number=A,Type=Float,Description=\"Allele Frequency, for each ALT allele, in the same order as listed\"> ##INFO=<ID=BQ,Number=1,Type=Integer,Description=\"Locus base Quality\"> ##INFO=<ID=MP,Number=1,Type=Integer,Description=\"Multiple mapped reads depth\"> ##INFO=<ID=IC,Number=1,Type=Integer,Description=\"Insertion reads count\"> ##INFO=<ID=DC,Number=1,Type=Integer,Description=\"Deletion reads count\"> ##INFO=<ID=PC,Number=1,Type=Integer,Description=\"Plus strand reads count\"> ##INFO=<ID=MC,Number=1,Type=Integer,Description=\"Minus strand reads count\"> ##INFO=<ID=SC,Number=1,Type=Integer,Description=\"Startpoint count\"> ##INFO=<ID=Ref,Number=1,Type=String,Description=\"Reference allele\"> ##INFO=<ID=Alt,Number=1,Type=String,Description=\"Alternative allele\"> #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Main fields are: - CHROM: chromosome name; - POS: nucleotide position in chromosome (1-based); - ID: sequence variant ID; - Ref: reference allele; - Alt: alternative allele; - QUAL: variant quality; - FILTER: variant filter; - INFO: detail information for variant as explained above; - FORMAT: ignored; 2. *.genome_update (main result) The updated genome sequence in FASTA format according to sequence alignment result or sequence variant result. 3. *.pos_change The genome coordinate conversion file between reference genome and updated genome/new genome. #Chr OldStart OldEnd NewStart NewEnd Main fields are: - Chr: chromosome name; - OldStart: start position in reference genome; - OldEnd: end position in reference genome; - NewStart: start position in updated genome/new genome; - NewEnd: end position in updated genome/new genome; - Strand: strand information ("+" or "-") for updated genome/new genome comparing with reference genome; 4. *.anno_update The annotation file for updated genome/new genome. The file format is as reference annotation file. 5. *.RefChr_NewChr.svg The genome comparison picture between reference genome and updated/new genome. Gene names and gene regions were plotted for reference and updated/new genome. The identical regions were connected by blue quadrilateral while the inconsistent regions were connected by purple quadrilateral. 6. *.psl The genome alignment result by BLAT. This file can be produced when assign "-n" parameter. 7. *.filter.psl The filtered genome alignment result. This file can be produced when assign "-n" parameter. Note: In default, the error messages or running messages are printed in screen. You can redirect the screen output into a file. TEST & EXAMPLE ============== RGAAT has been tested in three ways based on sequence alignment result, sequence variant result and genome comparison. 1 Run RGAAT using sequence alignment result For sequence alignment result, we first use bowtie2 (v2.2.4) to align DNA-seq data to reference genome with default parameters and obtained *.sam file. Then the *.sam file was converted to *.bam file and sorted according to the coodinate (*.sorted.bam) by samtools (Version 1.2). 1.1 For SAM file, the commond line and running messages are as following: Commond line: nohup perl RGAAT.pl -g coco_cp.v1.fa -a coco_cp.v1.tbl -s coco_cp.sam -o coco_cp.sam > coco_cp.sam.nohup & Running messages: Start time: Sun May 22 19:01:53 2016. RGAAT.pl -a coco_cp.v1.tbl -g coco_cp.v1.fa -s coco_cp.sam -o coco_cp.sam Start identify sequence variant: Sun May 22 19:01:53 2016. Processed 1000000 reads. Processed 2000000 reads. Processed 3000000 reads. Processed 4000000 reads. Processed 5000000 reads. Processed 6000000 reads. Processed 7000000 reads. Processed 8000000 reads. Processed 9000000 reads. #Total reads is 9402591. #Total filtered reads is 9294498 #The filter reads are as following: Unmap 61279 PCR_duplicate 0 Quality_control 46814 Multiple 0 Short_anchor 0 Low_quality 0 Sequence variant identified: Sun May 22 20:40:03 2016. Read genome file: Sun May 22 20:40:04 2016. Update genome: Sun May 22 20:40:04 2016. Genome comparison: Sun May 22 20:40:04 2016. End time: Sun May 22 20:40:07 2016. 1.2 For BAM file, the commond line and running messages are as following: Commond line: nohup perl RGAAT.pl -g coco_cp.v1.fa -a coco_cp.v1.tbl -s coco_cp.sorted.bam -o coco_cp.bam > coco_cp.bam.nohup & Running messages: Start time: Sun May 22 19:16:50 2016. RGAAT.pl -a coco_cp.v1.tbl -o coco_cp.bam -g coco_cp.v1.fa -b coco_cp.sorted.bam Start identify sequence variant: Sun May 22 19:16:50 2016. Processed 1000000 reads. Processed 2000000 reads. Processed 3000000 reads. Processed 4000000 reads. Processed 5000000 reads. Processed 6000000 reads. Processed 7000000 reads. Processed 8000000 reads. Processed 9000000 reads. #Total reads is 9402591. #Total filtered reads is 9294498 #The filter reads are as following: Unmap 61279 PCR_duplicate 0 Quality_control 46814 Multiple 0 Short_anchor 0 Low_quality 0 Sequence variant identified: Sun May 22 20:24:31 2016. Read genome file: Sun May 22 20:24:32 2016. Update genome: Sun May 22 20:24:32 2016. Genome comparison: Sun May 22 20:24:32 2016. End time: Sun May 22 20:24:35 2016. 2 Run RGAAT using sequence variant result To obtian the sequence variant result, we first use samtools and gatk to identify sequence variant. For samtools, we used the following commond lines (samtools: Version 1.2, bcftools: Version 1.2): samtools mpileup -d 100000 -L 100000 -m 3 -v -f coco_cp.correct.fa coco_cp.sorted.bam --output coco_cp.sorted.vcf& bcftools call -o coco_cp.sorted.call.vcf -O z -c coco_cp.sorted.vcf& bcftools view -O v -o coco_cp.samtools.vcf coco_cp.sorted.call.vcf& For gatk, we used the following commond lines (GATK: v3.3-0-g37228af, picard-tools: Version 1.119): java -jar CreateSequenceDictionary.jar R= coco_cp.correct.fa O= coco_cp.correct.dict & java -jar MarkDuplicates.jar INPUT=coco_cp.sorted.bam OUTPUT=coco_cp.sorted.markdup.bam METRICS_FILE=coco_cp.sorted.markdup.metrics ASSUME_SORTED=true & java -jar AddOrReplaceReadGroups.jar INPUT=coco_cp.sorted.markdup.bam OUTPUT=coco_cp.sorted.markdupnew.bam SORT_ORDER=coordinate RGPL=illumina RGLB=pairend180bp RGSM=coco_cp RGPU=NONE & java -jar GenomeAnalysisTK.jar -U ALL --filter_reads_with_N_cigar -T RealignerTargetCreator -R coco_cp.correct.fa -I coco_cp.sorted.markdupnew.bam -o coco_cp.sorted.markdupnew.intervals & java -jar GenomeAnalysisTK.jar -U ALL --filter_reads_with_N_cigar -T IndelRealigner -R coco_cp.correct.fa -I coco_cp.sorted.markdupnew.bam -targetIntervals coco_cp.sorted.markdupnew.intervals -o coco_cp.sorted.markdupnew2.bam & java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R coco_cp.correct.fa -I coco_cp.sorted.markdupnew2.bam -o coco_cp.gatk.vcf -stand_call_conf 20 -stand_emit_conf 20 & 2.1 For VCF file obtained by gatk, the commond line and running messages are as following: Commond line: perl RGAAT.pl -g coco_cp.v1.fa -a coco_cp.v1.tbl -v coco_cp.gatk.vcf -o coco_cp.gatk > coco_cp.gatk.nohup & Running messages: Start time: Sun May 22 19:20:12 2016. RGAAT.pl -v coco_cp.gatk.vcf -g coco_cp.v1.fa -o coco_cp.gatk -a coco_cp.v1.tbl Reads sequence variant file: Sun May 22 19:20:12 2016. Read genome file: Sun May 22 19:20:12 2016. Update genome: Sun May 22 19:20:12 2016. Genome comparison: Sun May 22 19:20:13 2016. End time: Sun May 22 19:20:15 2016. 2.2 For VCF file obtained by samtools, the commond line and running messages are as following: Commond line: nohup perl RGAAT.pl -g coco_cp.v1.fa -a coco_cp.v1.tbl -v coco_cp.samtools.vcf -o coco_cp.samtools > coco_cp.samtools.nohup & Running messages: Start time: Sun May 22 19:48:16 2016. RGAAT.pl -o coco_cp.samtools -a coco_cp.v1.tbl -v coco_cp.samtools.vcf -g coco_cp.v1.fa Reads sequence variant file: Sun May 22 19:48:16 2016. Read genome file: Sun May 22 19:48:19 2016. Update genome: Sun May 22 19:48:19 2016. Genome comparison: Sun May 22 19:48:19 2016. End time: Sun May 22 19:48:22 2016. 2.3 For five table file obtained by manually inspecting sequence alignment file, the commond line and running messages are as following: Commond line: perl RGAAT.pl -g coco_cp.v1.fa -a coco_cp.v1.tbl -v coco_cp.true.vcf -o coco_cp.true > coco_cp.true.nohup & Running messages: Start time: Sun May 29 13:06:01 2016. RGAAT.pl -o coco_cp.true -v coco_cp.true.vcf -a coco_cp.v1.tbl -g coco_cp.v1.fa Reads sequence variant file: Sun May 29 13:06:01 2016. Read genome file: Sun May 29 13:06:01 2016. Update genome: Sun May 29 13:06:01 2016. Genome comparison: Sun May 29 13:06:01 2016. End time: Sun May 29 13:06:04 2016. 3 Run RGAAT by genome comparison Commond line: perl RGAAT.pl -g coco_cp.v1.fa -a coco_cp.v1.tbl -n coco_cp.v2.fa -o coco_cp.blat > coco_cp.blat.nohup & Running messages: Start time: Sun May 22 20:46:29 2016. RGAAT.pl -n coco_cp.v2.fa -a coco_cp.v1.tbl -o coco_cp.blat -g coco_cp.v1.fa Genome comparison by BLAT: Sun May 22 20:46:29 2016. Loaded 154535 letters in 1 sequences Searched 154740 bases in 1 sequences BLAT result filter: Sun May 22 20:46:31 2016. Identify sequence variant: Sun May 22 20:46:31 2016. Annotation transfer: Sun May 22 20:46:31 2016. #Raw element: CDS 93 exon 54 gene 141 rRNA 8 tRNA 40 #Complete element: CDS 93 exon 54 gene 141 rRNA 8 tRNA 40 #Check element: Genome comparison: Sun May 22 20:46:33 2016. End time: Sun May 22 20:46:36 2016. 4. IGAAT result files (partial) 4.1 Figure - coco_cp.true.Cp_Cp.svg. genome comparison between reference and updated/new genome. Gene names and gene regions were plotted for reference and updated/new genome. The identical regions were connected by blue quadrilateral while the inconsistent regions were connected by purple quadrilateral. 4.2 Tables Table 1 - coco_cp.true.pos_change (partial) #Chr OldStart OldEnd NewStart NewEnd Strand Cp 1 186 1 186 NewGenome + Cp 187 464 188 465 NewGenome + Cp 465 520 467 522 NewGenome + Cp 521 3179 524 3182 NewGenome + Cp 3180 6977 3184 6981 NewGenome + SUBPROGRAM ========== We also provide a subprogram named sam2variant.pl for identification sequence variants using sequence alignment file (SAM format) and a subprogram named variant_filter.pl for filter sequence variants according to the attributes of sequence variant, such as base quality, allele frequency, multiple mapping reads percent, reads depth, and alternative allele depth. CONTACT ======= If you have any question or suggestion, please contact us. Wanfei Liu & Shuangyang Wu & Qiang Lin Email: <[email protected]> & <[email protected]> & <[email protected]>
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.