加载中…
访客
加载中…
个人资料
fanyucai
fanyucai
  • 博客等级:
  • 博客积分:0
  • 博客访问:347,706
  • 关注人气:340
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
博文
(2018-12-14 12:49)
分类: 文献推荐
1:在目标区域捕获中如下文献比较了几种常用方法,其中HaloPlex方法在比对完成后是不需要去除PCR冗余的。
Samorodnitsky E , Datta J , Jewell B M , et al. Comparison of Custom Capture for Targeted Next-Generation DNA Sequencing[J]. The Journal of Molecular Diagnostics, 2015, 17(1):64-75.


阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
(2018-12-11 21:56)
分类: 生物信息学
1:gtf转bed
在使用RSeQC(http://rseqc.sourceforge.net/#对RNAseq数据进行评估的时候,需要将gtf转换为bed文件,网上尽管有很多工具,但是使用起来转换过后都不在是标准的bed12格式,这里推荐使用QuickRNASeq(https://sourceforge.net/projects/quickrnaseq/files/)管道中的脚本就可以参考文献如下:
Zhao S, Xi L, Quan J, et al. QuickRNASeq lifts large-scale RNA-seq data analyses to the next level of automation and interactive visualization[J]. BMC genomics, 2016, 17(1): 39.
阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
(2018-12-11 21:50)
分类: 生物信息学
python的版本选择建议是2.7.9
建立index的时候是gff3的文件,该文件可以从gencode下载,默认是端测序,也接受双末端测序
# Compute Psi values for control sample and knockdown sample
miso --run mm9/pickled/SE data/control.bam --output-dir SE/control/ --read-len 35 --paired-end 250 15 --use-cluster
miso --run mm9/pickled/SE data/knockdown.bam --output-dir SE/knockdown/ --read-len 35 --paired-end 250 15 --use-cluster


## Summarize the output (only run this once --run finished!)This will create a "summary" directory in SE/control/ and in SE/knockdown/
summarize_miso --summarize-samples SE/control/ SE/control/
summarize_miso --summarize-samples SE/knockdown/ SE/knockdown/

## Detect differentially expressed isoforms between "control" and "knockdown"
## This will compute Bayes factors and delta Psi values between the samples
## and place the results in the directory SE/comparisons/control_vs_knockdown
compare_miso --compare-samples SE/control/ SE/knockdown/ SE/comparisons/

--paired-end后面跟的参数是插入片段大小及其方差可以通过一下命令进行计算
exon_utils --get-const-exons Mus_musculus.NCBIM37.65.gff --min-exon-size 1000 --output-dir exons/
pe_utils --compute-insert-len sample.bam Mus_musculus.NCBIM37.65.min_1000.const_exons.gff --output-dir insert-dist/
阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
分类: 生物信息学
比对并排序:
hisat2 -x GRCh38.p12.genome.fa -p 20 -1 sample1_R1.fq.gz -2 sample1_R2.fq.gz -S sample1.sam

=====================
当然比对也可以使用tophat:
tophat2 -T -G gencode.v29.annotation.gtf --transcriptome-index gencode.v29.transcripts -o out GRCh38.p12.genome sample1_R1.fq.gz sample1_R2.fq.gz
======================
samtools view -@ 20 -uS sample1.sam | samtools sort -@ 20 -o sample1.bam

统计获得比对结果:
samtools flagstat sample1.bam

############using stringtie to assemble RNA-Seq alignments into potential transcripts
stringtie sample1.bam -p 20 -G gencode.v29.annotation.gtf -l sample1 -F 0.01 -m 100 -o sample1.gtf
多个样本单独拼接完成后,你需要手动生产一个文本文件,该文件包含了生产的新的gtf文件的路径

##############Merged GTF and compare use cuffcompare to compare gtf file
stringtie --merge -G /data/Database/gencode_hsa/gencode.v29.annotation.gtf -p 20 -o AML20/stringtie_merge.gtf AML20.merge.list

gffcompare -R -r /data/Database/gencode_hsa/gencode.v29.annotation.gtf -o ALL21/strtcmp ALL21/stringtie_merge.gtf

###############Estimate transcript abundances and create table counts
htseq-count --format=bam --stranded=no --type=gene ERR315333.bam /data02/database/gencode/gencode.v29.annotation.gtf >DGE/ERR315333.counts
htseq-count --format=bam --stranded=no --type=gene ERR315486.bam /data02/database/gencode/gencode.v29.annotation.gtf >DGE/ERR315486.counts
htseq-count --format=bam --stranded=no --type=gene ALL21.bam /data02/database/gencode/gencode.v29.annotation.gtf >DGE/ALL21.counts
htseq-count --format=bam --stranded=no --type=gene AML20.bam /data02/database/gencode/gencode.v29.annotation.gtf >DGE/AML20.counts
###############join
join DGE/ALL21.counts DGE/ERR315333.counts | join - DGE/ERR315486.counts >ALL21_counts.tsv
echo "GeneID ALL21 ERR315333 ERR315486" > header.txt
cat header.txt ALL21_counts.tsv | grep -v "__" | perl -ne 'chomp $_; $_ =~ s/\s+/\t/g; print "$_\n"' > ALL21.counts.final.tsv

join DGE/AML20.counts DGE/ERR315333.counts | join - DGE/ERR315486.counts >AML20counts.tsv
echo "GeneID AML20 ERR315333 ERR315486" > header.txt
cat header.txt AML20counts.tsv | grep -v "__" | perl -ne 'chomp $_; $_ =~ s/\s+/\t/g; print "$_\n"' >AML20.counts.final.tsv

阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
(2018-11-26 11:50)
标签:

sv

分类: 生物信息学
1:关于结构变异分析的一篇综述性文件,感觉这篇文献最大的特点就是在末尾依据不同的分析软件对SV不同类型的检测分类:
Guan P, Sung W K. Structural variation detection using next-generation sequencing data: a comparative technical review[J]. Methods, 2016, 102: 36-49.






2:关于千人基因组分析结构变异的文章,极具参考价值,最起码也是使用了很多分析软件:
Sudmant P H, Rausch T, Gardner E J, et al. An integrated map of structural variation in 2,504 human genomes[J]. Nature, 2015, 526(7571): 75.

3:万人基因组发在PNAS上,这篇文章关于比较了7款主流的分析SV的分析软件for deletion: Pindel (8), DELLY (9), GenomeSTRiP (10), BreakDancer (11), LUMPY (12), MatchClip2 (13), Manta (14); for insertion: Pindel (8), Manta (14).,最终认为Manta效果比较好。对于CNV, we compared the performance of 5 software: cn.mops (16), CNVnator (17), GenomeSTRiP (10), MatchClip2 (13), Canvas (18).认为Canvas比较好
Telenti A, Pierce L C T, Biggs W H, et al. Deep sequencing of 10,000 human genomes[J]. Proceedings of the National Academy of Sciences, 2016, 113(42): 11901-11906.

4:个人建议目前分析结构变异单一软件很难解决问题,一般需要多个软件。或者某一软件只针对几种结构变异而不是全部,在大的方面结构变异可以尝试使用lumpy以及Delly,针对CNV可以使用cnvkit。

5:还有推荐的关于CNV的综述性文献里面列举了发现CNV的几种常规方法,一些常规算法不需要懂,但名字还是要知道,因为在分析过程中会产生假阳性,还要对数据进行过滤。
Pirooznia M, Goes F S, Zandi P P. Whole-genome CNV analysis: advances in computational approaches[J]. Frontiers in genetics, 2015, 6: 138.

6:分析参考文献;Dixon J R, Xu J, Dileep V, et al. Integrative detection and analysis of structural variation in cancer genomes[J]. Nature genetics, 2018, 50(10): 1388.
阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
标签:

vcf

lcr

low

complexity

region

分类: 生物信息学
参考文章链接:
http://bcb.io/2014/05/12/wgs-trio-variant-evaluation/
参考文献:
Li, H. Toward better understanding of artifacts in variant calling from high-coverage samples[J]. Bioinformatics, 2014, 30(20):2843-2851.

坐标转换使用liftOver,参考链接:http://www.bio-info-trainee.com/990.html,转换过程中需要https://hgdownload-test.gi.ucsc.edu/goldenPath/hg38/liftOver/hg38ToHg19.over.chain.gz文件

已经基于hg38或者hg37的bed文件:https://github.com/lh3/varcmp/tree/master/scripts 如果需要基于hg19版本的就需要你使用上面的工具进行转换了

有了bed文件可以使用GATK的VariantFiltration功能,vcftools也是可以的
阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
分类: Perl
#!/usr/bin/perl
use strict;
use warnings;
use Spreadsheet::WriteExcel;
my $file=shift;

open(IN,$file);
my %hash;
my @sample;
while()
{
    chomp;
    my @array=split(/\t/,$_);
    if($_=~/^#/)
    {
        for(my $i=1;$i<=$#array;$i++)
        {
           $sample[$i-1]=$array[$i]; 
        }
    }
    else
    {
        for(my $i=0;$i<=$#array;$i++)
        {
            
            if($i>0)
            {
                if($array[$i]>0)
                {
                    $hash{$sample[$i-1]}{$array[0]}=$array[$i];
                }
            }
        }
    }
}
my $workbook = Spreadsheet::WriteExcel->new('newexcel.xls');
foreach my $key(sort {$a cmp $b} keys %hash)
{
    my $worksheet = $workbook->add_worksheet($key);
    my $j=0;
    foreach my $key2(sort {$a cmp $b} keys %{$hash{$key}})
    {
        $j++;
        $worksheet->write("A$j",$key2);
        $worksheet->write("B$j",$hash{$key}{$key2});
    }
}
$workbook->close();                         
阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
分类: 医学相关
关于dbsnp数据库说明(https://www.ncbi.nlm.nih.gov/snp)相关文件说明参见:https://www.ncbi.nlm.nih.gov/variation/docs/human_variation_vcf/,该数据库可以分为两部分:
没有临床特征的VCF数据:ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/
有临床表征的VCF数据:ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar

截止到2018年4月23日已经更新到最新版本为151,兼容GATK版本的vcf文件:ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/GATK/00-All.vcf.gz。
注意该vcf文件不包含等位基因频率信息

Genome Aggregation Database(简称gnomAD)突变频率数据,库该数据库提供的数据集包括125,748个个体的全外显子组测序数据和15,708个个体的全基因组测序数据,这些数据来源于各种疾病研究项目及大型人群测序项目。数据下载链接:http://gnomad.broadinstitute.org/downloads 目前版本是2.1

ExAC(http://exac.broadinstitute.org)是一个外显子测序变异数据库,最新版本是2017.02.27升级的1.0版本,数据下载ftp地址:ftp://ftp.broadinstitute.org/pub/ExAC_release/release1/ExAC.r1.sites.vep.vcf.gz

COSMIC
千人基因组http://www.internationalgenome.org
* Pilot Analysis
   * A map of human genome variation from population-scale sequencing Nature 467, 1061–1073 (28 October 2010)
* Phase 1 Analysis
   * An integrated map of genetic variation from 1,092 human genomes Nature 491, 56–65 (01 November 2012)
* Phase 3 Analysis
   * A global reference for human genetic variation Nature 526, 68–74 (01 October 2015)
   * An integrated map of structural variation in 2,504 human genomes Nature 526, 75–81 (01 October 2015)

phase1和phase3的一些区别:http://www.internationalgenome.org/category/phase-3/

由于GATK在call变异的过程中会用到局部重比对(主要是出现indels附近往往会出现较多的snp)往往会用到已知的高质量的Indel文件例如:
* Mills_and_1000G_gold_standard.indels.b37.vcf
* 1000G_phase1.indels.b37.vcf
这个Mills是Mills R E, Luttig C T, Larkins C E, et al. An initial map of insertion and deletion (INDEL) variation in the human genome[J]. Genome research, 2006, 16(9): 1182-1190.这篇文章的作者以上文件的下载链接:

Omni是Illumina一款芯片(https://www.illumina.com/Documents/products/datasheets/datasheet_gwas_roadmap.pdf)包含基因型信息
一般可以在GATK的VariantRecalibrator步骤可以使用(https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantrecalibration_VariantRecalibrator.php)这一步可以使用基于hapmap的vcf数据

hapmap_3.3.b37.sites.vcf 
1000G_omni2.5.b37.sites.vcf

https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_mutect_Mutect2.php#--germline-resource
如何使用https://gatkforums.broadinstitute.org/gatk/discussion/11136/how-to-call-somatic-mutations-using-gatk4-mutect2

Mutect2基于haplotypes开发
Sentieon
自行设计的panel大小至少要有3M否则会影响TMB的检测
Buchhalter I, Rempel E, Endris V, et al. Size Matters: Dissecting Key Parameters for PanelBased Tumor Mutational Burden (TMB) Analysis[J]. International journal of cancer, 2018.

Microsatellite instability(微卫星不稳定性MSI)又称短串联重复序列(STR)是指基因组中小于10个核苷酸的简单重复序列,以2个核苷酸组成的重复序列最为丰富,主要位于基因组的非编码区。DNA错配修复(mutations of mismatch repair)是细胞复制后的一种修复机制,起维持DNA复制保真度,控制基因变异的作用。MMR系统由MMR基因编码的一系列的MMR蛋白组成。MMR基因的突变或者修饰(如甲基化)可以导致MMR蛋白的缺乏。

1997年美国国家癌症研究所工作组推荐通过检测基因组上的5个微卫星位点(BAT26, BAT25, D2S123, D5S346, D17S250)的不稳定性来判定微卫星不稳定性成都。MSI-H2个或2个以上位点的不稳定定为微卫星高度不稳定,1个位点不稳定定为微卫星低度不稳定MSI-L,无位点出现不稳定性MSS。如果是你自己设计的panel的话,位点不稳定性占你总panel数目的30%则认为MSI-H。
阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
标签:

somatic

snv

mutect2

gatk

umi

分类: 医学相关
Kalatskaya I, Trinh Q M, Spears M, et al. ISOWN: accurate somatic mutation identification in the absence of normal tissue controls[J]. Genome Medicine, 2017, 9(1):59.

变异检测可以分为三类:
single nucleotide variant (SNV), insertion and deletion (indel), and structural variant (SV, including copy number variation, duplication, translocation, etc.)典型的SNV和小的indel一般小等于10bp

比对:
Illumina(BWA)、TMAP (for Ion Torrent reads) for DNA reads
splice-aware aligners such as TopHat and STAR for RNA sequencing

关于在比对之前是否要做数据质控的问题,在这里做数据质控也只是去掉接头序列。因为很多变异检测的软件都是基于位点的检测策略,所以整条reads的质量情况不是那么重要,另外局部重比对也就是BQSR (base quality score recalibration)。基于PCR扩增的数据不需要在数据比对后去除PCR冗余。

tumor-normal变异检测模式

基于启发式算法的编译检测算法有VarScan2, qSNP, Shimmer, RADIA, SOAPsnv, and VarDict
加入genotype analysis的分析软件有SomaticSniper, FaSD- somatic, SAMtools, JointSNVMix2, Virmid, SNVSniffer, Seurat, and CaVEMan,这些软件一般使用在低覆盖的数据分析中(WGS, WES, or targeted sequencing with low depth),但是对低频突变不敏感
基于Haplotype-based strategy检测策略不需要局部重比对,因为该变异检测方法是基于reads组装后的结果进行编译检测的,这样的软件有Platypus, HapMuC, LocHap, FreeBayes, and MuTect2 。
基于机器学习方法的软件MutationSeq, SomaticSeq, SNooPer, and BAYSIC

如果是高测序覆盖深度的低频突变建议使用Strelka, MuTect, LoFreq, EBCall, deepSNV, LoLoPicker, and MuSE,启发式算法的软件对于发现低频突变也有较好的效果(1% variant calling with VarDict) and (< 5% variant calling with VarScan2)

Single-sample 变异检测模式

SNVMix2, Shearwater, SPLINTER, SNVer, OutLyzer, and Pisces这些软件都可以进行单样本变异检测但是不能区分somatic and germline
ISOWN, SomVarIUS, and SiNVICT可以提供单样本的变异检测但是同时也可同时区分somatic and germline,ISOWN软件是依赖于MuTect2,随后依赖somatic (COSMIC) and germline mutations (ExAC and dbSNP)来做进一步区分,OutLyzer, Pisces, ISOWN, SomVarIUS, SiNVICT已经被应用到靶向测序的应用

UMI-based variant calling
一般低频突变定义为((VAF ≤5%) )目标就是排出测序错误Illumina(0.01–0.1 ),目前给予UMI分析的软件有三款: DeepSNVMiner, MAGERI, and smCounter
其中 DeepSNVMiner, MAGERI, and smCounter输入都是原始数据而只有smCounter的输入是BAM格式,基于PGM平台已经有了一个处理UMI的插件TVC。此外还有一个开源的软件Fgbio。Illumina建议DNA输入量30ng 测序层数40000X 中值覆盖度可达到~2500X 敏感性变异检测为:0.4%

RNA-seq variant calling

基于RNA数据变异检测的软件有RADIA, Seurat, VarDict, VarScan2, SNPiR, and eSNVdetect,但是RADIA and Seurat 需要整合RADIA and Seurat 的DNA数据

2014年Genome in a Bottle Consortium简称(GIAB)通过整合多种测序科技和比对分析软件公布了NA12878 cell line 细胞系高质量可信的变异检测结果
Zook J M, Chapman B, Wang J, et al. Integrating human sequence data sets provides a resource of benchmark SNP and indel genotype calls[J]. Nature Biotechnology, 2014, 32(3):246-51.
阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
(2018-11-02 12:54)
分类: 生物信息学

最近在学习弄spark很感觉挫败,没有搞明白,现把一些学习经历分享一下。

1:在安装hadoop与spark的时候需注意版本搭配,GATK也给出了基于谷歌云的相关版本可参加如下信息:

2:在做人的vcf注释的时候可以下载的人注释文件下载网站:
UCSC:http://hgdownload.cse.ucsc.edu/gbdb/hg19/
GATK:https://software.broadinstitute.org/gatk/download/bundle

3:运行BwaSpark的时候需要建立全新的参考基因组index文件,生产的文件后缀名是img格式
gatk BwaMemIndexImageCreator -I ucsc.hg19.fasta -O ucsc.hg19.fasta.img

4:基于GATK的序列比对输入并不是fastq而是未排序的bam文件,可以通过FastqToSam程序将fastq转成未比对的bam或者sam
gatk FastqToSam -F1 /data/test/B-23_R1.fq.gz -F2 /data/test/B-23_R2.fq.gz -O B-23unaligned.bam -PL illumina -SM B-23 -RG B-23

5:spark对文件小于60G的数据文件加速效果并不明显,甚至还没有在本地运行的快(个人经验谨慎参考),运行spark的时候最关键的参数有2个,executor-cores(Number of cores for the executors)与executor-memory(Memory per Executor)还有1个参数--driver-memory 这个值的设置一般在4-8G之间为宜,过小会报错

6:对于普通的16cpu64G内存的节点,在设计参数的时候首先要避免满负荷运转,所以可用内存以及可用cpu都要适当减少。对于如下条件设定(http://blog.cloudera.com/blog/2015/03/how-to-tune-your-apache-spark-jobs-part-2/):
--num-executors 6 --executor-cores 15 --executor-memory 63G
不如设定为:
--num-executors 17 --executor-cores 5 --executor-memory 19G
以上设置是博客里提供的,但是在我的测试里还是貌似第一种方式比较好,原因未明
华为提供的配置建议:
--num-executors这个参数可以不设置
7:关于内存和cpu的分配可以参考如下网址:
https://blog.csdn.net/fansy1990/article/details/54314249/
https://github.com/broadinstitute/gatk/wiki/Troubleshooting-Spark
阅读  ┆ 评论  ┆ 转载 ┆ 收藏 
  

新浪BLOG意见反馈留言板 电话:4006900000 提示音后按1键(按当地市话标准计费) 欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

新浪公司 版权所有