跟着Genome Research学数据分析：msmc2软件分析种群历史动态

用户7010445

发布于 2023-08-23 10:39:42

5020

发布于 2023-08-23 10:39:42

论文

A butterfly pan-genome reveals that a large amount of structural variation underlies the evolution of chromatin accessibility

https://genome.cshlp.org/content/early/2022/09/15/gr.276839.122

A butterfly pan-genome reveals that a large amount of structural variation underlies the evolution of chromatin accessibility.pdf

论文的代码基本都公开了，很好的学习材料，今天的推文我们学习一下其中关于分析种群历史动态的代码

参考链接 https://github.com/StevenVB12/Genomics/blob/master/PAN_SV_chromatin_Genomics/PSMC/Run_PSMC.sh

这里我用到的是酵母的数据，分析这个种群历史动态是把同一个物种的二代测序数据比对到自己的参考基因组，然后去分析，为什么是这样做暂时还搞不明白，今天的推文只是跑通代码，还有好多不明白的地方需要去看

这里比对是bwa samtools是sam bam 格式转换

picard做重复序列标记

bcftools做变异检测，论文里提供的代码是使用samtools做的变异检测，但是我这边一直报错，可能是samtools的版本和论文中使用的版本不一致

bwa index N44.genome.fa
bwa mem -t 12 -M N44.genome.fa SRR4074394_1.fastq.gz SRR4074394_2.fastq.gz -o SRR4074394.sam

samtools view -f 0x02 -q 20 -b SRR4074394.sam -o SRR4074394.bam

samtools sort SRR4074394.bam -o SRR4074394.sorted.bam -@ 4

picard MarkDuplicates I=SRR4074394.sorted.bam O=SRR4074394.sorted.md.bam Remove_Duplicates=true ASSUME_SORTED=true M=SRR4074394.sorted_dup_metrics.txt
samtools index SRR4074394.sorted.md.bam

这以上都是一些比较常规的代码

~/biotools/software.package/bcftools-1.17/bcftools mpileup -q 20 -Q 20 -C 50 --fasta-ref N44.genome.fa -r chrI SRR4074394.sorted.md.bam | bcftools call -c -V indels | ~/biotools/software.package/msmc-tools-master/bamCaller.py 100 mask.bed.gz | gzip -c SRR4074394.chrI.vcf.gz

这一步的代码用到了bamCaller.py这个脚本来自于

https://github.com/stschiff/msmc-tools

这里生成了mask.bed.gz 和 SRR4074394.chrI.vcf.gz文件，这里是只有一条染色体

接下来的代码生成了 SRR4074394.chrI.SCAF.txt 文件

~/biotools/software.package/msmc-tools-master/generate_multihetsep.py --mask=mask.bed.gz SRR4074394.chrI.vcf.gz > SRR4074394.chrI.SCAF.txt

然后用这个SRR4074394.chrI.SCAF.txt 文件去做分析

msmc2_Linux -t 4 -o SRR4074394 SRR4074394.chrI.SCAF.txt

分析完得到了一个 SRR4074394.final.txt 文件

内容是

time_index left_time_boundary right_time_boundary lambda
0 0 5.9236e-07 28.6566
1 5.9236e-07 1.28527e-06 28.6566
2 1.28527e-06 2.09581e-06 1484.13
3 2.09581e-06 3.04393e-06 19546.7
4 3.04393e-06 4.153e-06 80653.5
5 4.153e-06 5.45033e-06 197381
6 5.45033e-06 6.96788e-06 218423
7 6.96788e-06 8.74303e-06 139242
8 8.74303e-06 1.08195e-05 111799
9 1.08195e-05 1.32485e-05 99175.2
10 1.32485e-05 1.60898e-05 90119.1
11 1.60898e-05 1.94134e-05 86319.4
12 1.94134e-05 2.33012e-05 87526
13 2.33012e-05 2.78489e-05 93032.5
14 2.78489e-05 3.31686e-05 103957
15 3.31686e-05 3.93913e-05 124967
16 3.93913e-05 4.66704e-05 167097
17 4.66704e-05 5.5185e-05 166593
18 5.5185e-05 6.5145e-05 16108.1
19 6.5145e-05 7.67958e-05 1004.82
20 7.67958e-05 9.04242e-05 113.209
21 9.04242e-05 0.000106366 28.6566
22 0.000106366 0.000125014 28.6566
23 0.000125014 0.000146828 28.6566
24 0.000146828 0.000172344 28.6566
25 0.000172344 0.000202192 28.6566
26 0.000202192 0.000237106 28.6566
27 0.000237106 0.000277947 28.6566
28 0.000277947 0.000325721 28.6566
29 0.000325721 0.000381604 28.6566
30 0.000381604 0.000446974 28.6566
31 0.000446974 inf 28.6566

这个是用来作图的数据吗？暂时还搞不明白

论文里还写到

We scaled the PSMC estimates using a generation time of 0.25 yr
and a mutation rate of 2 × 10−9 as estimated for H. melpomene
(i.e., spontaneous Heliconius mutation rate corrected for selective
constraint

这里提到的突变率怎么获得呢？

还有好多问题需要解决

还找到了一些其他参考资料

GitHub - stschiff/msmc-tools: Tools and Utilities for msmc and msmc2 https://github.com/stschiff/msmc-tools

GitHub - stschiff/msmc: Implementation of the multiple sequential markovian coalescent https://github.com/stschiff/msmc