跟着NAR学数据分析：利用二代测序数据估计拟南芥基因组大小

论文

Pushing the limits of HiFi assemblies reveals
centromere diversity between two Arabidopsis
thaliana genomes

https://academic.oup.com/nar/article/50/21/12309/6858746?login=false

2个拟南芥NAR.pdf

代码链接

https://github.com/frabanal/A.thaliana_CLR_vs_HiFi/tree/v1.0

下载二代测序数据

~/biotools/kingfisher/bin/kingfisher get -r ERR8666067 -m ena-ftp

对测序数据进行过滤

cutadapt -j 16 -q 20,15 -b TruSeq1=AGATCGGAAGAGC -b Nextera1=CTGTCTCTTATACACATCT -b Nextera1rc=AGATGTGTATAAGAGACAG -B TruSeq2=AGATCGGAAGAGC -B Nextera2=CTGTCTCTTATACACATCT -B Nextera2rc=AGATGTGTATAAGAGACAG --trim-n --minimum-length 75 -o ERR8666067_cutadapt_R1.fastq.gz --paired-output ERR8666067_cutadapt_R2.fastq.gz ERR8666067_1.fastq.gz ERR8666067_2.fastq.gz

这里怎么确定具体用到的是那种接头序列一直没有搞明白

bwa比对

bwa index at_ChrC_ChrM_phiX74.fa
bwa mem -t 16 -R "@RG\tID:ERR8666067\tSM:ERR8666067" at_ChrC_ChrM_phiX74.fa ERR8666067_cutadapt_R1.fastq.gz ERR8666067_cutadapt_R2.fastq.gz -o ERR8666067.sam
samtools sort -@ 16 -n -O BAM -o ERR8666067.sorted.bam ERR8666067.sam

samtools 基本统计

samtools stats ERR8666067.sorted.bam

提取没有比对上的数据

r1没有比对上r2比对上

samtools view -@ 16 -b -f 4 -F 264 ERR8666067.sorted.bam -O BAM -o ERR8666067.unmap_map.bam

r1比对上，r2没有比对上

samtools view -@ 16 -b -f 8 -F 260 ERR8666067.sorted.bam -O BAM -o ERR8666067.map_unmap.bam

r1 r2 都没有比对上

samtools view -@ 16 -b -f 12 -F 256 ERR8666067.sorted.bam -O BAM -o ERR8666067.unmap_unmap.bam

这里有一步没有看懂是什么意思

echo "Merging both combinations of SE-mapped..."
rm $outputTMP/$ACC.$ref.TMP1_unmapped.bam
$SAMTOOLS merge $outputTMP/$ACC.$ref.TMP1_unmapped.bam $outputTMP/$ACC.$ref.unmap_map.bam $outputTMP/$ACC.$ref.map_unmap.bam

echo "Discarding supplementary alignments..."
$SAMTOOLS view -b -F 2048 $outputTMP/$ACC.$ref.TMP1_unmapped.bam > $outputTMP/$ACC.$ref.TMP2_unmapped.bam
$SAMTOOLS sort -n -o $outputTMP/$ACC.$ref.single_unmapped.bam $outputTMP/$ACC.$ref.TMP2_unmapped.bam

Discarding supplementary alignments...是什么意思暂时没有搞明白

后续只用r1和r2都没有比对上的reads

bam转fastq

bamToFastq -i ERR8666067.unmap_unmap.bam -fq ERR8666067.unmap_unmap.R1.fastq -fq2 ERR8666067.unmap_unmap.R2.fastq

估计基因组大小

cat ERR8666067.unmap_unmap.*.fastq | jellyfish count /dev/fd/0 -C -o ERR8666067.21mer -m 21 -t 16 -s 5G
jellyfish histo -h 3000000 -o ERR8666067.21mer.histo ERR8666067.21mer

findGSE的链接

https://github.com/schneebergerlab/findGSE/tree/master

library(findGSE)
findGSE(histo="ERR8666067.21mer.histo",size=21,outdir="ERR8666067_21mer")

findGSE initialized...
   Info: histo file provided as ERR8666067.21mer.histo
   Info: size k set as 21
   Info: output folder set as ERR8666067_21mer
   Info: expected coverage of homozygous k-mers set as 0
   Info: het observed set as false ==> heterozygous fitting not asked. 

Iterative fitting process for sample ERR8666067.21mer.histo started... 
    Size 21 at itr 1
    Info: min_valid_pos:  33 
    Info: signal error border:  33 
    Info: hom_xfit_left  for hom fitting at itr  1 :  89 
    Info: hom_xfit_right for hom fitting at itr  1 :  129 
       Fitting has to be repeated: sizek 21 at itr 1
    Info: hom_xfit_left  for hom fitting at itr  1 :  89 
    Info: hom_xfit_right for hom fitting at itr  1 :  129 
    Size 21 at itr 2
    Info: hom_xfit_left  for hom fitting at itr  2 :  111 
    Info: hom_xfit_right for hom fitting at itr  2 :  435 
    Size 21 at itr 3
    Info: hom_xfit_left  for hom fitting at itr  3 :  234 
    Info: hom_xfit_right for hom fitting at itr  3 :  541 
    Warning: data does not follow assumed distribution anymore at itr 3, fitting stopped.
Iterative fitting done.

Genome size estimate for ERR8666067.21mer.histo: 142872866 bp.

论文中的数据是 we estimated the genome size to be 143.12 Mb

本次的重复结果比论文中的小一点

image.png image.png

推文记录的是自己的学习笔记，很可能存在错误，请大家批判着看