最近在学习当下很火的环状DNA(eccDNA)的建库技术以及生信分析,特此记录。
No.1 去低质量&去接头
这一步需要注意的是低质量碱基的quality cutoff, 2018年NC的paper中他们-q参数取的10,而我们自己通常喜欢取20.
接头序列对于很多入门新手来说不容易搞清楚,因此推荐使用trim_galory软件,因为它会自动识别adapter序列,参考代码如下:
trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 --paired $rawdatadir/${i}/${i}_1.fq.gz $rawdatadir/${i}/${i}_2.fq.gz --gzip -o $trimdir
No.2 比对 Alignment
这一步一定要注意,由于后续Circle-Map需要用到好几种bam文件进行realign,所以一定按照下述步骤进行:
(1) fastq比对生成bam (bam0)
time bwa mem -q -t 20 hg38.fa sample.trimed.1.fastq.gz sample.trimed.2.fastq.gz | samtools view -bS > sample.bam
(2) 对初始sam进行按Read名排序(bam1)
time samtools sort -n -o sample_nameSorted.bam sample.bam
sample.nameSort.bam这个文件是Circle-Map用来提取环形DNA候选reads的文件,并用来估计先验分布的参数的。
(3) 对初始sam进行按最左边映射坐标排序 (bam2)
time samtools sort -o sample_LeftReadSorted.bam sample.bam
time samtools index sample_LeftReadSorted.bam
Cirle-Map使用sample_LeftReadSorted.bam这个文件来计算基因组覆盖率指标。
(4) 使用Circle-Map提取环形DNA circular_structural_reads的bam(bam3)。主要是split_reads和discordant_reads
#首先是Circle-Map的安装
conda install -c bioconda Circle-Map
#这里需要安装biopython=1.7.7版本,原本的biopython1.78在使用时报错了
pip install biopython=1.7.7
##PS: pip和pip3是一样的用法,区别是如果python的版本是2,那么只能使用pip,而如果是python3,则都可以使用,虽然是这么说,但是我的python3.9.6的版本下使用pip3安装失败,没有找到这个命令,而pip使用没有问题
#都装好后,输入bam1进行reads提取,产生bam3
time Circle-Map ReadExtractor -i sample_nameSorted.bam -o sample_circular_candidate.bam
#Circle-Map ReadExtractor 这一步会产生无法找到bam1的报错,即使使用了bam1的绝对路径,但是这个报错不会影响正常的bam3的输出。
#接下来对bam3进行按坐标排序并建索引
time samtools sort -o sample_circular_candidate.sorted.bam sample_circular_candidate.bam
time samtools index sample_circular_candidate.sorted.bam
No.3 环形DNA(eccDNA)检测——生信输出
总结这一步所需要的前面生成的3个bam,以及其输入顺序为 -i bam3, -qbam bam1. -sbam bam2
Circle-Map Realign -t 4 -i sample_circular_candidate.sorted.bam -qbam sample_nameSorted.bam -sbam sample_LeftReadSorted.bam -fasta hg38.fa -o sample_unknown_circle.bed