different result for population and individual

[Elevensyl]

@stschiff Hi stschiff！
I recently used msmc2 to calculate historical effective population size. I used two individuals from each of five populations of sheep for analysis, but the results for individual sheep and the populations were different. I don't know what went wrong, I would be grateful if you could answer.
here is my script:
step1:

cat chr.list | while read i; do samtools depth -r $i sample_name.sorted.markdup.bam | awk '{sum += $3} END {print sum / NR}' > sample_name.depth_$i.txt; done && echo "sample_name done" || echo "sample_name fail"

step2

sample="sample_name"
ref="~/workplace/sheep_genomic/GCF_016772045.2_ARS-UI_Ramb_v3.0_genomic.fna"
depth_dir="~/workplace/demographic_inference/msmc2/depth_all"
mkdir -p vcf_mask_${sample}
cat chr.name | while read chr; do
    echo "Processing ${sample} ${chr}..."

    depth=$(cat ${depth_dir}/${sample}.depth_${chr}.txt)

    bcftools mpileup -C 50 -Ou -r ${chr} -f ${ref} ${sample}.sorted.markdup.bam --threads 8 | \
    bcftools call -c -V indels --threads 8 | \
  ~/software/msmc-tools/bamCaller.py ${depth} vcf_mask_${sample}/${sample}_${chr}.mask.bed.gz | \
    bgzip -c > vcf_mask_${sample}/${sample}_${chr}.vcf.gz

    if [ $? -eq 0 ]; then
        echo "${sample} ${chr} done"
    else
        echo "${sample} ${chr} failed"
    fi

done

echo "All done for ${sample}!"

step3 for population

for i in $(cat chr.name)
do
~/software/msmc-tools/generate_multihetsep.py --chr ${i} \
    --mask ${mask_input}CB_${i}.mask.bed.gz \
    --mask ${mask_input}CC_${i}.mask.bed.gz \
    --mask ${mask_input}COK1_${i}.mask.bed.gz \
    --mask ${mask_input}COK4_${i}.mask.bed.gz \
    --mask ${mask_input}1408_${i}.mask.bed.gz \
    --mask ${mask_input}1409_${i}.mask.bed.gz \
    --mask ${mask_input}165M_${i}.mask.bed.gz \
    --mask ${mask_input}166M_${i}.mask.bed.gz \
    --mask ${mask_input}ERR12389563_${i}.mask.bed.gz \
    --mask ${mask_input}ERR12389564_${i}.mask.bed.gz \
    ${vcf_input}CB_${i}.vcf.gz \
    ${vcf_input}CC_${i}.vcf.gz \
    ${vcf_input}COK1_${i}.vcf.gz \
    ${vcf_input}COK4_${i}.vcf.gz \
    ${vcf_input}1408_${i}.vcf.gz \
    ${vcf_input}1409_${i}.vcf.gz \
    ${vcf_input}165M_${i}.vcf.gz \
    ${vcf_input}166M_${i}.vcf.gz \
    ${vcf_input}ERR12389563_${i}.vcf.gz \
    ${vcf_input}ERR12389564_${i}.vcf.gz \
    > ${output_dir}input_${i}.msmcInput.txt && echo "${i} done" || echo "${i} fail" 
done

step3 for individual

for i in $(cat chr.name)
do
~/software/msmc-tools/generate_multihetsep.py --chr ${i} \
    --mask ${mask_input}sample_name_${i}.mask.bed.gz \
    ${vcf_input}sample_name_${i}.vcf.gz \
    > ${output_dir}sample_name_${i}.msmcInput.txt && echo "sample_name ${i} done" || echo "sample_name ${i} fail" 
done

step4 for population

~/software/msmc2-master/build/release/msmc2 --skipAmbiguous -I 0,1 -i 25 -t 4 -o pop/kurdi input_*.msmcInput.txt && echo "kurdi done" || echo "kurdi fail"
~/software/msmc2-master/build/release/msmc2 --skipAmbiguous -I 2,3 -i 25 -t 4 -o pop/cok input_*.msmcInput.txt && echo "cok done" || echo "cok fail"
~/software/msmc2-master/build/release/msmc2 --skipAmbiguous -I 4,5 -i 25 -t 4 -o pop/asiatic input_*.msmcInput.txt && echo "asiatic done" || echo "asiatic fail"
~/software/msmc2-master/build/release/msmc2 --skipAmbiguous -I 6,7 -i 25 -t 4 -o pop/ermani input_*.msmcInput.txt && echo "ermani done" || echo "ermani fail"
~/software/msmc2-master/build/release/msmc2 --skipAmbiguous -I 8,9 -i 25 -t 4 -o pop/cyprus input_*.msmcInput.txt && echo "cyprus done" || echo "cyprus fail"

step4 for individual

for i in {CB,CC,COK1,COK4,1408,1409,165M,166M,ERR12389563,ERR12389564}
do
~/software/msmc2-master/build/release/msmc2 --skipAmbiguous -i 25 -t 4 -o msmc2/${i} \
${i}_NC_056054.1.msmcInput.txt ${i}_NC_056055.1.msmcInput.txt ${i}_NC_056056.1.msmcInput.txt ${i}_NC_056057.1.msmcInput.txt ${i}_NC_056058.1.msmcInput.txt \
${i}_NC_056059.1.msmcInput.txt ${i}_NC_056060.1.msmcInput.txt ${i}_NC_056061.1.msmcInput.txt ${i}_NC_056062.1.msmcInput.txt ${i}_NC_056063.1.msmcInput.txt \
${i}_NC_056064.1.msmcInput.txt ${i}_NC_056065.1.msmcInput.txt ${i}_NC_056066.1.msmcInput.txt ${i}_NC_056067.1.msmcInput.txt ${i}_NC_056068.1.msmcInput.txt \
${i}_NC_056069.1.msmcInput.txt ${i}_NC_056070.1.msmcInput.txt ${i}_NC_056071.1.msmcInput.txt ${i}_NC_056072.1.msmcInput.txt ${i}_NC_056073.1.msmcInput.txt \
${i}_NC_056074.1.msmcInput.txt ${i}_NC_056075.1.msmcInput.txt ${i}_NC_056076.1.msmcInput.txt ${i}_NC_056077.1.msmcInput.txt ${i}_NC_056078.1.msmcInput.txt ${i}_NC_056079.1.msmcInput.txt
done

[stschiff]

Hmm, I'm not quite sure I can follow this. Both plots show 5 lines. But your "individual" run actually consists of 10 runs over 10 individuals, so there should be 10 curves, right? For the "population", you choose 5 sets of haplotypes, with -I 0,1 corresponding to the two haplotypes of the first diploid sample, -I 2,3 to the two haplotypes of the second individual, and so on. So how many samples do you have, and how do you divide them up into populations?

I think you have 10 samples, which would mean 20 haplotypes.

[Elevensyl]

Hmm, I'm not quite sure I can follow this. Both plots show 5 lines. But your "individual" run actually consists of 10 runs over 10 individuals, so there should be 10 curves, right? For the "population", you choose 5 sets of haplotypes, with -I 0,1 corresponding to the two haplotypes of the first diploid sample, -I 2,3 to the two haplotypes of the second individual, and so on. So how many samples do you have, and how do you divide them up into populations?

I think you have 10 samples, which would mean 20 haplotypes.

I'm sorry that I didn't make my meaning clear. You understood it correctly. There are indeed 10 individuals. It's just that when I drew the graph, I only chose five to represent the population because the other five lines have the same trend.

I understand what you mean. My -I parameter setting is wrong. I have 10 individuals, so the -I parameter should be 0, 1, 2, 3...19, but I set it with 0,1,2,....9.