论文
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
推文记录的是自己的学习笔记，很可能存在错误，请大家批判着看
