小男孩‘自慰网亚洲一区二区,亚洲一级在线播放毛片,亚洲中文字幕av每天更新,黄aⅴ永久免费无码,91成人午夜在线精品,色网站免费在线观看,亚洲欧洲wwwww在线观看

分享

項(xiàng)目一:使用二代測(cè)序數(shù)據(jù)進(jìn)行基因組組裝(局部組裝)

 生物_醫(yī)藥_科研 2019-02-13

項(xiàng)目數(shù)據(jù):

  • 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)

工具:

策略:

  • 將測(cè)序的二代reads使用BWA比對(duì)到參考基因組,分成不同的窗口,按窗口進(jìn)行局部組裝,然后合并。

預(yù)備知識(shí):

  • 能用熟練使用 Perl 和 shell 寫腳本
  • 會(huì)熟練使用 PBS 提交任務(wù)
  • BWA使用方法
  • IGV使用方法
  • SOAPdenovo使用方法

具體步驟:

1.前期準(zhǔn)備

復(fù)制代碼
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
復(fù)制代碼

2.比對(duì)

/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.可視化的前處理

復(fù)制代碼
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
復(fù)制代碼

4.按窗口分類reads

寫perl腳本,提取SAM中的reads名稱,去fastQ里提取原始reads,按窗口分類,為下一步的局部組裝做準(zhǔn)備。

 

5.局部組裝

 

 

 

局部組裝的問(wèn)題:

已經(jīng)有兩批人沒(méi)組出來(lái)了,局部組裝大多不可能組裝出完整的100K窗口,因?yàn)槎蛄衦eads太短,重復(fù)序列太多,重復(fù)序列會(huì)導(dǎo)致連接中斷,一個(gè)窗口會(huì)出現(xiàn)很多片段,而且也沒(méi)有方法將其繼續(xù)連接起來(lái),所以他們都半途而廢了。

后續(xù)可能會(huì)遇到的情況,必須借助后期的分析手段,將諸多片段連接成完整的序列。

杜發(fā)的文章,完全是在無(wú)參考基因組的情況下,denovo組裝,利用多種手段,才將零碎的序列組裝成完整的基因組。

 

其他:

PCRfree,就是DNA樣品不是通過(guò)PCR進(jìn)行擴(kuò)增的,防止PCR中的錯(cuò)誤造成影響.

map,就是確定基因在染色體中的位置.

眾多軟件都可以利用 fastq.gz 文件, less也可以查看此類型的文件.

 

問(wèn)題:

  • fastq文件里都有些什么?  參見(jiàn) fastQ格式

簡(jiǎn)介:fastQ是二代測(cè)序下機(jī)數(shù)據(jù)的格式, 它存儲(chǔ)了核酸序列和相應(yīng)質(zhì)量評(píng)價(jià)的信息,該格式包含四行:

第一行: @HWUSI-EAS100R:6:73:941:1973#0/1 ; 以@開頭,后面是reads的ID以及其他信息,例如上例中 HWUSI-EAS100R代表Illmina設(shè)備名稱,6代表flowcell中的第六個(gè)lane,73代表第六個(gè)lane中的第73個(gè) tile,941:1973代表該read在該tile中的x:y坐標(biāo)信息;#0,若為多樣本的混合作為輸入樣本,則該標(biāo)志代表樣本的編號(hào),用來(lái)區(qū)分個(gè)樣本中的reads;/1代表paired end中的前一個(gè)read。

第二行: GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTT ; read序列

第三行: +HWUSI-EAS100R:6:73:941:1973#0/1 ; 以“+”開頭,跟隨者該read的名稱(一般于@后面的內(nèi)容相同),但有時(shí)可以省略,但“+”一定不能省。

第四行: !''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC6 ; 以“+”開頭,跟隨者該read的名稱(一般于@后面的內(nèi)容相同),但有時(shí)可以省略,但“+”一定不能省。

本項(xiàng)目中的原始reads格式如下: (每條read均為150 bp)

復(fù)制代碼
@ST-E00126:176:HL7H5CCXX:1:1101:23368:6319 1:N:0:GTCCGC
TATCGCTTGCTCAAACGCTGGGCAACTAGTCTCTAGTCTTGGTCGAGTGTGTGGTGGACTTGGTCGTGGACATGCTTGGTTCTTAGATCGTGTTTTGTATTCAGGGTTGCTTGTACCCTTTCTTTCTTGCACTGTCCATATTCTAATGCA
+
AAAAFKKFKKKKKKKKKKKKKKKKKKKKKKKFFKKKKKKKKKKKKFAFKK,,FKK,A<F<KAFKKKFKKKFKKKKKKKFKF<,,,AKKKKFK,FFKKKKKKFA<KK7<,<<,,7<AFFFKKKAFFKKKKKKK77FFA,,<FKKKKAAFAF
復(fù)制代碼

補(bǔ)充:雙端測(cè)序時(shí),fastq文件中,R1 和 R2 兩個(gè)文件中的reads是如何排列的?

復(fù)制代碼
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
復(fù)制代碼

 

  • 參考基因組()是個(gè)什么玩意?
復(fù)制代碼
>Chr1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
CTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCTAAACCC
TAAACCCTAACCCTAAACCCTAACCCTAAACCCTAAACCCTAAACCCTAA
ACCCTAAACCCTAAACAGCTGACAGTACGATAGATCCACGCGAGAGGAAC
CGGAGAGACAACGGGATCCAGGCGCCAGCGACGGATCCGGGCGAGAGGGG
AGTGGAGATCATGGATCCGTGCGGGAGGGGAAGAAGTCGCCGAATCCGAC
......
復(fù)制代碼

chr1一共有865419行, 每一行50個(gè)堿基, 最后一行23個(gè)堿基, 一共43270923個(gè)堿基, 約為43Mb.

該染色體的頭部 尾部 以及中間有少量的NNNN組成的gap, 是沒(méi)有組裝出來(lái)的區(qū)域.

 

分秈、粳稻兩個(gè)亞種, 一共12對(duì)染色體, 基因組大小:430Mb,  預(yù)測(cè)有32000~56000個(gè)基因,

 

  • BWA 比對(duì)的原理, 以及之后生成的 SAM 文件該如何解讀, SAM格式是如何存儲(chǔ)比對(duì)信息的?

 

 

額外參考:

[Perl] Getopt 函數(shù)來(lái)接收用戶參數(shù)的使用

awk命令

    本站是提供個(gè)人知識(shí)管理的網(wǎng)絡(luò)存儲(chǔ)空間,所有內(nèi)容均由用戶發(fā)布,不代表本站觀點(diǎn)。請(qǐng)注意甄別內(nèi)容中的聯(lián)系方式、誘導(dǎo)購(gòu)買等信息,謹(jǐn)防詐騙。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容,請(qǐng)點(diǎn)擊一鍵舉報(bào)。
    轉(zhuǎn)藏 分享 獻(xiàn)花(0

    0條評(píng)論

    發(fā)表

    請(qǐng)遵守用戶 評(píng)論公約

    類似文章 更多