Tag Archives: 突变

五 05

用VEP对vcf格式的突变数据进行注释

Posted on by

VEP是国际三大数据库之一的ENSEMBL提供的，也是非常主流和方便，但它是基于perl语言的，所以在模块方面可能会有点烦人。跟snpEFF一样，也是对遗传变异信息提供更具体的注释，而不仅仅是基于位点区域和基因。如果你熟悉外显子联盟这个数据库EXAC(ExAC.r0.3.sites.vep.vcf.gz)，你可以下载它所有的突变记录数据，看看它对每个变异位点到底注释了些什么，它就是典型的用VEP来注释的。 Continue reading →

五 05

用snpEFF对vcf格式的突变数据进行注释

Posted on by

这个软件比较重要，尤其是对做遗传变异相关研究的，很多人做完了snp-calling后喜欢用ANNOVAR来进行注释，但是那个注释还是相对比较简单，只能得到该突变位点在基因的哪个区域，那个基因这样的信息，如果想了解更具体一点，就需要更加功能化的软件了，snpEFF就是其中的佼佼者，而且是java平台软件，非常容易使用！而且它的手册写的非常详细：http://snpeff.sourceforge.net/SnpEff_manual.html Continue reading →

十二 29

用Mutation-Assessor软件来看突变位点对基因或者蛋白功能的影响

Posted on by
这是一个在线工具，非常好用，对snp位点进行注释，看看该突变是否影响蛋白功能，一定要收藏！！！

官网：http://mutationassessor.org/

也应该是有standalone版本，我没有去找，不过网页版就很好用了，只需要复制粘贴进去自己想分析的数据，按照一定的格式即可，比如：
clipboard
该软件就能分析出该突变位点发生在BRCA2这个基因上面，对氨基酸的改变也能写出来，对蛋白功能的改变等选项都是可以自由定制化的。
输入数据非常广泛：
The server accepts list of variants, one variant per line, plus optional text describing your variants,
in genomic coordinates, "+" strand assumed :
<genome build>,<chromosome>,<position>,<reference allele>,<substituted allele>
Genome build is optional (build 18 assumed), accepted values: 'hg18' and 'hg19'
Examples:

hg19,13,32912555,G,T   BRCA2
hg18,7,55178574,G,A   GBM
7,55178574,G,A   GBM

or in protein space: <protein ID> <variant> <text>, where <protein ID> can be :

1. Uniprot protein accession (i.e. EGFR_HUMAN)
2. NCBI Refseq protein ID (i.e. NP_005219)

examples:

EGFR_HUMAN R521K
EGFR_HUMAN R98Q Polymorphism
EGFR_HUMAN G719D disease
NP_000537 G356A
NP_000537 G360A dbSNP:rs35993958
NP_000537 S46A Abolishes phosphorylation

ID types can be mixed in one list in any way.

 

十 16

3000多份水稻全基因组测序数据共享-主要是突变数据

Posted on by

感觉最近接触的生物信息学知识越多，越对大数据时代的到来更有同感了。现在的研究者，其实很多都可以自己在家里做了，大量的数据基本都是公开的， 但是一个人闭门造车成就真的有限，与他人交流的思想碰撞还是蛮重要的。

这里面列出了3000多份水稻全基因组测序数据，都共享在亚马逊云上面，是全基因组的双端测序数据，共3,024个水稻数据，比对到了五种不同的水稻参考基因组上面，而且主要是用GATK来找差异基因的。
而且，数据收集者还给出了一个snp calling的标准流程
我以前也是用这样的流程
SNP Pipeline Commands

1. Index the reference genome using bwa index

   /software/bwa-0.7.10/bwa index /reference/japonica/reference.fa

2. Align the paired reads to reference genome using bwa mem. 
   Note: Specify the number of threads or processes to use using the -t parameter. The possible number of threads depends on the machine where the command will run.

   /software/bwa-0.7.10/bwa mem -M -t 8 /reference/japonica/reference.fa /reads/filename_1.fq.gz /reads/filename_2.fq.gz > /output/filename.sam

3. Sort SAM file and output as BAM file

   java -Xmx8g -jar /software/picard-tools-1.119/SortSam.jar INPUT=/output/filename.sam OUTPUT=/output/filename.sorted.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=TRUE

4. Fix mate information

   java -Xmx8g -jar /software/picard-tools-1.119/FixMateInformation.jar INPUT=/output/filename.sorted.bam OUTPUT=/output/filename.fxmt.bam SO=coordinate VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=TRUE

5. Mark duplicate reads

   java -Xmx8g -jar /software/picard-tools-1.119/MarkDuplicates.jar INPUT=/output/filename.fxmt.bam OUTPUT=/output/filename.mkdup.bam METRICS_FILE=/output/filename.metrics VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=TRUE MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000

6. Add or replace read groups

   java -Xmx8g -jar /software/picard-tools-1.119/AddOrReplaceReadGroups.jar INPUT=/output/filename.mkdup.bam OUTPUT=/output/filename.addrep.bam RGID=readname PL=Illumina SM=readname CN=BGI VALIDATION_STRINGENCY=LENIENT SO=coordinate CREATE_INDEX=TRUE

7. Create index and dictionary for reference genome

   /software/samtools-1.0/samtools faidx /reference/japonica/reference.fa
   
   java -Xmx8g -jar /software/picard-tools-1.119/CreateSequenceDictionary.jar REFERENCE=/reference/japonica/reference.fa OUTPUT=/reference/reference.dict

8. Realign Target 

   java -Xmx8g -jar /software/GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar -T RealignerTargetCreator -I /output/filename.addrep.bam -R /reference/japonica/reference.fa -o /output/filename.intervals -fixMisencodedQuals -nt 8

9. Indel Realigner

   java -Xmx8g -jar /software/GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar -T IndelRealigner -fixMisencodedQuals -I /output/filename.addrep.bam -R /reference/japonica/reference.fa -targetIntervals /output/filename.intervals -o /output/filename.realn.bam 

10. Merge individual BAM files if there are multiple read pairs per sample

   /software/samtools-1.0/samtools merge /output/filename.merged.bam /output/*.realn.bam

11. Call SNPs using Unified Genotyper

   java -Xmx8g -jar /software/GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar -T UnifiedGenotyper -R /reference/japonica/reference.fa -I /output/filename.merged.bam -o filename.merged.vcf -glm BOTH -mbq 20 --genotyping_mode DISCOVERY -out_mode EMIT_ALL_SITES
十 09

对vcf突变数据与公开发表的进行比对

Posted on by
当我们对NGS数据call了snp之后一般会输出成vcf格式的数据，一行代表一个突变，例如
20      2451451 .       G       T       1939.77 .
AC=1;AF=0.500;AN=2;BaseQRankSum=-10.134;DP=239;Dels=0.00;FS=2.276;HaplotypeScore=0.0000;MLEAC=1;MLEAF=0.500;MQ=60.00;MQ0=0;MQRankSum=-0.258;QD=8.12;ReadPosRankSum=0.823;SOR=0.870
GT:AD:DP:GQ:PL  0/1:150,89:239:99:1968,0,3874
#前几列记录着该突变发生在第几号染色体以及该染色体的哪个坐标，我们的参考基因组在该位点是什么碱基，我们测到的突变成了什么碱基。
最后两列是测序深度以及正负测序深度，或者ref和allele的测序深度。
只有第8列是最复杂的，可以有高达几百个数据信息，取决于我们用什么样的软件来call的snp，以及call了snp之后用什么样的软件做的注释。
接下来我们还需要探究我们找到的突变是否在其它以及公开发表的数据库里面被找到过，所以可以下载非常多的公共数据库进行比对，我所见过的有一下一些，估计完全下载有0.5T
dbsnp144 （这个是ncbi提供的最权威的啦）
cgi69
ExAC.vcf.gz（这个是broadinstitute提供的外显子联盟）
Cosmic_v73.ann.vcf.gz （这个是癌症突变信息集）
finalTCGA.vcf.gz （TCGA计划也是癌症相关的）
icgc.vcf.gz
dbNSFP2.6vcf
SCLP.ann.vcf.gz
CCLE.ann.vcf.gz
ESP6500-SIv2.vcf.gz （Variants from the Exome Sequencing Project (ESP)）
adni-sum
safs-sum.indel.vcf.gz
gonl.vcf.gz
ssm.vcf.gz
ssi.vcf.gz
uk10k.vcf.gz
1000g-ph3v5.gff.gz  （千人基因组计划）
gwasCatalog.gff.gz  \
phewascatalog.gff.gz  \
dbgap-gwas.gff.gz  \
interproDomain.gff.gz \
clinvar.gff.gz \
RegulomeDB.gff.gz \
CancerGAMAdb.gff.gz \