项目一：使用二代测序数据进行基因组组装（局部组装）

项目数据：

  • kongyu_131_PCRfree_.CCAAT_L006_R1_001.fastq.gz （100X）（19G）
    kongyu_131_PCRfree_.CCAAT_L006_R2_001.fastq.gz （100X）（20G）
  • Y255_PCRfree_.TCCGC_L005_R1_001.fastq.gz （30X）（5.4G）
    Y255_PCRfree_.TCCGC_L005_R2_001.fastq.gz （30X）（6.0G）
  • all.chrs.con.fasta （364M） 日本晴  ( /public/genome/rice/msu/version_7.0/all.dir/all.chrs.con.fasta)

工具：

策略：

  • 将测序的二代reads使用BWA比对到参考基因组，分成不同的窗口，按窗口进行局部组装，然后合并。

预备知识：

  • 能用熟练使用 Perl 和 shell 写脚本
  • 会熟练使用 PBS 提交任务
  • BWA使用方法
  • IGV使用方法
  • SOAPdenovo使用方法

具体步骤:

1.前期准备

NP=`cat $PBS_NODEFILE | wc -l`
NN=`cat $PBS_NODEFILE | sort -u | wc -l`

cd $PBS_O_WORKDIR
export LD_LIBRARY_PATH=\$LD_LIBRARY_PATH:/public/software/htslib-1.3/lib
date

samtoolsdir=/public/software/samtools-1.3/bin
bwadir=/public/software/bwa-0.7.12-intel
dir=/public/scripts/liyan/scripts2016

sample=Y255
out=/public/home/zxli/Project/Project_Sliced_Assembly

ref=/public/genome/rice/msu/version_7.0/all.dir/all.chrs.con.fasta
fq1=/public/home/zxli/data/rice/Y255_PCRfree_.TCCGC_L005_R1_001.fastq.gz
fq2=/public/home/zxli/data/rice/Y255_PCRfree_.TCCGC_L005_R2_001.fastq.gz

2.比对

/public/software/bwa-0.7.12-intel/bwa index $ref
$bwadir/bwa mem -t $NP -f -M -R "@RG\tID:$sample\tLB:$sample\tSM:$sample\tPL:illumina\tPU:$sample" $ref $fq1 $fq2 > $out/bwamem_$sample.sam
date

3.可视化的前处理

samtools view -@ 10 -bS -F 4 ./contigs_sequence_align_to_public_genome.sam > ./contigs_sequence_align_to_public_genome.bam
samtools sort -@ 10 ./contigs_sequence_align_to_public_genome.bam contigs_sequence_align_to_public_genome.sorted
samtools index contigs_sequence_align_to_public_genome.sorted.bam contigs_sequence_align_to_public_genome.sorted.bai
samtools depth contigs_sequence_align_to_public_genome.sorted.bam >depth_reads.txt
wc -l depth_reads.txt > Coverage-aln_reads.txt
 
$samtoolsdir/samtools view -@ $NP -Sb $out/bwamem_$sample.sam -o $out/bwamem_$sample.bam
$samtoolsdir/samtools sort -@ $NP $out/bwamem_$sample.bam -o $out/bwamem_$sample.sorted.bam
$samtoolsdir/samtools index $out/bwamem_$sample.sorted.bam

4.按窗口分类reads

写perl脚本，提取SAM中的reads名称，去fastQ里提取原始reads，按窗口分类，为下一步的局部组装做准备。

 

5.局部组装

 

 

 

局部组装的问题：

已经有两批人没组出来了，局部组装大多不可能组装出完整的100K窗口，因为二代序列reads太短，重复序列太多，重复序列会导致连接中断，一个窗口会出现很多片段，而且也没有方法将其继续连接起来，所以他们都半途而废了。

后续可能会遇到的情况，必须借助后期的分析手段，将诸多片段连接成完整的序列。

杜发的文章，完全是在无参考基因组的情况下，denovo组装，利用多种手段，才将零碎的序列组装成完整的基因组。

 

其他:

PCRfree,就是DNA样品不是通过PCR进行扩增的,防止PCR中的错误造成影响.

map,就是确定基因在染色体中的位置.

众多软件都可以利用 fastq.gz 文件, less也可以查看此类型的文件.

 

问题:

简介:fastQ是二代测序下机数据的格式, 它存储了核酸序列和相应质量评价的信息,该格式包含四行:

第一行: @HWUSI-EAS100R:6:73:941:1973#0/1 ; 以@开头，后面是reads的ID以及其他信息，例如上例中 HWUSI-EAS100R代表Illmina设备名称，6代表flowcell中的第六个lane，73代表第六个lane中的第73个 tile，941:1973代表该read在该tile中的x：y坐标信息；#0，若为多样本的混合作为输入样本，则该标志代表样本的编号，用来区分个样本中的reads；/1代表paired end中的前一个read。

第二行: GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTT ; read序列

第三行: +HWUSI-EAS100R:6:73:941:1973#0/1 ; 以“+”开头，跟随者该read的名称（一般于@后面的内容相同），但有时可以省略，但“+”一定不能省。

第四行: !''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC6 ; 以“+”开头，跟随者该read的名称（一般于@后面的内容相同），但有时可以省略，但“+”一定不能省。

本项目中的原始reads格式如下: (每条read均为150 bp)

@ST-E00126:176:HL7H5CCXX:1:1101:23368:6319 1:N:0:GTCCGC
TATCGCTTGCTCAAACGCTGGGCAACTAGTCTCTAGTCTTGGTCGAGTGTGTGGTGGACTTGGTCGTGGACATGCTTGGTTCTTAGATCGTGTTTTGTATTCAGGGTTGCTTGTACCCTTTCTTTCTTGCACTGTCCATATTCTAATGCA
+
AAAAFKKFKKKKKKKKKKKKKKKKKKKKKKKFFKKKKKKKKKKKKFAFKK,,FKK,A<F<KAFKKKFKKKFKKKKKKKFKF<,,,AKKKKFK,FFKKKKKKFA<KK7<,<<,,7<AFFFKKKAFFKKKKKKK77FFA,,<FKKKKAAFAF

补充：双端测序时，fastq文件中，R1 和 R2 两个文件中的reads是如何排列的？

R1
@ST-E00126:176:HL7H5CCXX:1:1101:7638:1221 1:N:0:NTCCGC
@ST-E00126:176:HL7H5CCXX:1:1101:8572:1221 1:N:0:NTCCGC
R2
@ST-E00126:176:HL7H5CCXX:1:1101:7638:1221 2:N:0:NTCCGC
@ST-E00126:176:HL7H5CCXX:1:1101:8572:1221 2:N:0:NTCCGC

 

  • 参考基因组（）是个什么玩意? 
>Chr1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
CTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCC
TAAACCCTAACCCTAAACCCTAACCCTAAACCCTAAACCCTAAACCCTAA
ACCCTAAACCCTAAACAGCTGACAGTACGATAGATCCACGCGAGAGGAAC
CGGAGAGACAACGGGATCCAGGCGCCAGCGACGGATCCGGGCGAGAGGGG
AGTGGAGATCATGGATCCGTGCGGGAGGGGAAGAAGTCGCCGAATCCGAC
......

chr1一共有865419行, 每一行50个碱基, 最后一行23个碱基, 一共43270923个碱基, 约为43Mb.

该染色体的头部 尾部 以及中间有少量的NNNN组成的gap, 是没有组装出来的区域.

 

分籼、粳稻两个亚种, 一共12对染色体, 基因组大小:430Mb,  预测有32000～56000个基因,

 

  • BWA 比对的原理, 以及之后生成的 SAM 文件该如何解读, SAM格式是如何存储比对信息的?

 

 

额外参考：

[Perl] Getopt 函数来接收用户参数的使用

awk命令

posted @ 2016-06-20 18:39  Life·Intelligence  阅读(6036)  评论(0)  编辑  收藏  举报
TOP