搜索
查看: 4736|回复: 1

从基因组的fasta序列里面模拟测序的reads

[复制链接]

633

主题

1189

帖子

4054

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4054
发表于 2017-1-19 22:27:20 | 显示全部楼层 |阅读模式
节选自我的博客:http://www.bio-info-trainee.com/1081.html(模拟Y染色体测序判断,并比对到X染色体上面,看同源性)
下载参考基因组,或者一条染色体均可:
然后用下面的程序来模拟测序文件:

全部代码如下:
前提是你用我上面的perl代码自己生成了read1.fa和read2.fa这两个模拟测序文件!
[Shell] 纯文本查看 复制代码
cd tmp/chrX_Y/hg19/
wget [url=http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrX.fa.gz;]http://hgdownload.cse.ucsc.edu/g ... mosomes/chrX.fa.gz;[/url]
wget [url=http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrY.fa.gz;]http://hgdownload.cse.ucsc.edu/g ... mosomes/chrY.fa.gz;[/url]
gunzip chrX.fa.gz
gunzip chrY.fa.gz
~/biosoft/bwa/bwa-0.7.15/bwa index chrX.faperl simulate.pl chrY.fa ##[color=#333333][font=Consolas, &quot][size=14.4px]这个perl脚本在 [img]http://www.bio-info-trainee.com/wp-content/uploads/2015/10/tmp.png[/img] [/size][/font][/color]
~/biosoft/bwa/bwa-0.7.15/bwa mem -t 5 -M chrX.fa read*.fa >read.sam
samtools view -bS read.sam >read.bam
samtools flagstat read.bam
samtools sort -@ 5 -o read.sorted.bam  read.bam
samtools view -h -F4  -q 5 read.sorted.bam |samtools view -bS |samtools rmdup -  read.filter.rmdup.bam
samtools index read.filter.rmdup.bam
samtools mpileup -ugf ~/tmp/chrX_Y/hg19/chrX.fa  read.filter.rmdup.bam  |bcftools call -vmO z -o read.bcftools.vcf.gz




那么我们就可以通过自己的数据处理能力来探索一下X和Y染色体的同源区有多少,是哪里的问题!
首先下载X,Y染色体的fasta序列,在UCSC上面下载即可。
然后把X染色体构建bwa的索引
接着模拟一个Y染色体的测序数据,模拟的程序很简单,模拟Y染色体的测序片段(PE100,insert400)
最后把模拟测序数据比对到X染色体的参考,统计一下比对结果即可!
我自己看sam文件也发现真的同源性好高呀,总共就模拟了380万reads,就有120万是百分百比对上了。
所以对女性个体来说,测序判断比对到Y染色体是再正常不过的了。如果要判断性别,必须要找那些X,Y差异性区段!
对男性来说,更是如此!
本次测试涉及到的文件如下:

1.png

[Perl] 纯文本查看 复制代码
while(<>){
chomp;
$chrY.=uc $_;
}
$j=0;
open FH_L,">read1.fa";
open FH_R,">read2.fa";
foreach (1..4){
        for ($i=600;$i<(length($chrY)-600);$i = $i+50+int(rand(10))){
                $up = substr($chrY,$i,100);
                $down=substr($chrY,$i+400,100);
                next unless $up=~/[ATCG]/;
                next unless $down=~/[ATCG]/;
                $down=reverse $down;
                $down=~tr/ATCG/TAGC/;
                $j++;
                print FH_L ">read_$j/1\n";
                print FH_L "$up\n";
                print FH_R ">read_$j/2\n";
                print FH_R "$down\n";
        }
}
close FH_L;
close FH_R;













上一篇:线性回归中“回归”的含义
下一篇:TCGA的metadata简介
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-1-20 09:12:12 | 显示全部楼层
赞 原来还可以这么搞
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-6-1 14:04 , Processed in 0.024581 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.