项目数据:
工具:
策略:
预备知识:
具体步骤:
1.前期准备
NP=`cat $PBS_NODEFILE | wc -l`NN=`cat $PBS_NODEFILE | sort -u | wc -l`cd $PBS_O_WORKDIRexport LD_LIBRARY_PATH=\$LD_LIBRARY_PATH:/public/software/htslib-1.3/libdatesamtoolsdir=/public/software/samtools-1.3/binbwadir=/public/software/bwa-0.7.12-inteldir=/public/scripts/liyan/scripts2016sample=Y255out=/public/home/zxli/Project/Project_Sliced_Assemblyref=/public/genome/rice/msu/version_7.0/all.dir/all.chrs.con.fastafq1=/public/home/zxli/data/rice/Y255_PCRfree_.TCCGC_L005_R1_001.fastq.gzfq2=/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.samdate
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:GTCCGCTATCGCTTGCTCAAACGCTGGGCAACTAGTCTCTAGTCTTGGTCGAGTGTGTGGTGGACTTGGTCGTGGACATGCTTGGTTCTTAGATCGTGTTTTGTATTCAGGGTTGCTTGTACCCTTTCTTTCTTGCACTGTCCATATTCTAATGCA+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:NTCCGCR2@ST-E00126:176:HL7H5CCXX:1:1101:7638:1221 2:N:0:NTCCGC@ST-E00126:176:HL7H5CCXX:1:1101:8572:1221 2:N:0:NTCCGC
>Chr1NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAACCCTAAACCCTAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACAGCTGACAGTACGATAGATCCACGCGAGAGGAACCGGAGAGACAACGGGATCCAGGCGCCAGCGACGGATCCGGGCGAGAGGGGAGTGGAGATCATGGATCCGTGCGGGAGGGGAAGAAGTCGCCGAATCCGAC
......
chr1一共有865419行, 每一行50个碱基, 最后一行23个碱基, 一共43270923个碱基, 约为43Mb.
该染色体的头部 尾部 以及中间有少量的NNNN组成的gap, 是没有组装出来的区域.
分籼、粳稻两个亚种, 一共12对染色体, 基因组大小:430Mb, 预测有32000~56000个基因,
额外参考:
联系客服