HBV分析

anneng

四医大HBV分析记录
1.使用bbmerge合并R1 R2 下面脚本的意思是使用find找到样本名称然后使用这个样本名称传递给parallel并发处理

find ../all_data/*_L001_R1_001.fastq.gz | sed 's/_L001_R1_001.fastq.gz$//' | parallel 'bbmerge.sh in1=../all_data/{}_L001_R1_001.fastq.gz in2=../all_data/{}_L001_R2_001.fastq.gz out={}.fastq  outu1={}.R1.umerged outu2={}.R2.unmerged'

发现327个样本中有几个样本 R1 和 R2 的数量不一致针对这些样本使用spades进行组装取最长的序列进行第二步
因为涉及到组装无法进行混合样品的分析把这些样本当作单样本处理
将所有的fastq转成fasta(blast只识别fasta)

parallel 'seqtk seq -a {}> {.}.fasta' ::: *.fastq

2.使用blast 对样本中的序列进行分型得到每个样本中各种分型的序列数量
构建blast数据库
从hbvdb下载的参考序列有一个类别是RF 例如 https://www.ncbi.nlm.nih.gov/nucleotide/EU871985.1?report=genbank&log$=nuclalign&blast_rank=1&RID=Z8DW1MY8016 这个序列 NCBI没有标识类型 hbvdb将其注释为了BC重组型我们当前先把这种RF的去掉

makeblastdb -in all_hbvdb_Genomes.fas -dbtype nucl

blastn -task blastn -max_target_seqs 1 -query ../0-merging-pe/100_S42.fasta -db ../hbvdb/all_hbvdb_Genomes.fas -num_threads 10 -out 100_S42.m8 -outfmt 6

nohup bash -c "find ../0-merging-pe/*.fasta | sed 's/.fasta$//' |  parallel --joblog ./logs -j40 blastn -task blastn -max_target_seqs 1 -query ../0-merging-pe/{}.fasta -db ../hbvdb/A-H/HBV_A_H.fas -out {/}.m8 -outfmt 6 " &

3.比对

nohup bash -c "find ../all_data/*_L001_R1_001.fastq.gz | sed 's/_L001_R1_001.fastq.gz$//' | parallel 'bwa mem -M AB033556_hbc_type_C.fasta {}_L001_R1_001.fastq.gz {}_L001_R2_001.fastq.gz > {/}.sam' " &

nohup parallel "samtools view -bF 4 {} > {/.}.bam" ::: ./sam/*.sam &
parallel samtools sort {} -o {.}.sorted.bam ::: *.bam

4.call

nohup parallel "lofreq indelqual {} --dindel -f ../3-mapping/AB033556_hbc_type_C.fasta -o {/.}.sorted.dindel.bam " ::: ../3-mapping/bam/*.sorted.bam &

nohup parallel "lofreq call {} --call-indels -f ../3-mapping/AB033556_hbc_type_C.fasta -o {/.}.vcf " ::: *.bam &

5.分析单倍型

find /ceph_disk2/siyida_327_sample/3-mapping/sam/ -name "*.sam" -exec basename \{} .sam \; | sed 's/.sam$//' |parallel 'java -jar clique-snv.jar -m snv-illumina -in /ceph_disk2/siyida_327_sample/3-mapping/sam/{}.sam'

anneng

spades的组装

/home/bioinfo/miniconda2/envs/assembly/bin/spades.py      -1      /ceph_disk3/hbv/HBV_illumina/106/106_S46_L001_R1_001.fastq      -2      /ceph_disk3/hbv/HBV_illumina/106/106_S46_L001_R2_001.fastq      -o      /ceph_disk3/hbv/HBV_illumina/106/spades

anneng

https://www.sciencedirect.com/science/article/pii/S1386653218300970
Frequency of hepatitis B surface antigen variants (HBsAg) in hepatitis B virus genotype B and C infected East- and Southeast Asian patients: Detection by the Elecsys HBsAg II assay

anneng

https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0172101
Ultra-deep sequencing reveals high prevalence and broad structural diversity of hepatitis B surface antigen mutations in a global population
https://github.com/spabinger/HBV_data_publication_2016_07
an MHR variant was defined as a nucleotide sequence change in the S gene region (encoding amino acids 99 to 170) with an allele frequency >5% (in both sequencing directions) and at least 3 variant reads present on the forward as well as on the reverse strand.

anneng

https://sci-hub.st/10.1159/000361076
Hepatitis B Virus Drug Resistance Tools:
One Sequence, Two Predictions
www.genafor.org/services.php

HIV-GRADE HBV

文章提到了一些工具用于分型、耐药、免疫逃逸的分析

anneng

Genetic Diversity of Hepatitis B Virus
Strains Derived Worldwide: Genotypes,
Subgenotypes, and HBsAg Subtypes

https://sci-hub.st/10.1159/000080872
对HBV进行进化树分析里面也提到血清型和基因型之间的复杂的对应关系。
涉及的软件：
DNADIST and NEIGHBOR from the Phylip program package version 3.53

PUZZLE

Bootstrap on 1,000 replicas was performed with SEQBOOT, DNADIST, NEIGHBOR, and CONSENSE from the Phylip package.

anneng

Global Occurrence of Clinically Relevant Hepatitis B Virus
viruses-12-01344-v3.pdf

从蛋白序列预测血清型

anneng

https://www.aimspress.com/article/doi/10.3934/microbiol.2020024?viewType=HTML
突变可能造成的影响这个论文做了一个总结

anneng

https://www.nature.com/articles/s41598-019-43524-9
Illumina and Nanopore methods for whole genome sequencing of hepatitis B virus (HBV)

anneng

https://www.frontiersin.org/articles/10.3389/fmicb.2020.616023/full
Comprehensive Analysis of Clinically Significant Hepatitis B Virus Mutations in Relation to Genotype, Subgenotype and Geographic Region
使用公开数据分析HBV的突变
Table_1_Comprehensive Analysis of Clinically Significant Hepatitis B Virus Mutations in Relation to Genotype, Subgenotype and Geographic Region.XLSX

这个表格的格式可以作为分析的模板
行是样本列是突变的位置或者重要图标的代号

anneng

https://www.hiv.lanl.gov/content/sequence/ENTROPY/entropy.html

香农熵计算器

anneng

https://zhanglab.ccmb.med.umich.edu/I-TASSER/.
结构预测

anneng

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7229894/
四医大肖老师提供的一个文章这个使用clone测序方法对HBV的全长进行了测序

组装：Contig-Express 和Codon Code Aligner
序列对齐：MEGAX Clustal X

anneng

Inference with viral quasispecies diversity indices: Clonal and
NGS approaches

对突变频率香农熵做了详细分析

anneng

https://www.yacinemahdid.com/shannon-entropy-from-theory-to-python/
香农熵的python实现

anneng

https://elifesciences.org/articles/61803
The haplotypes for each sample were reconstructed for each gene segment using a previously published pipeline (Cacciabue et al., 2020). In brief, FastQC (Andrews, 2010) was used for quality assurance of the NGS paired-end raw reads followed by BBtools (Bushnell, 2014), for removing and filtering adapters and low-quality reads. Bowtie2 (Langmead and Salzberg, 2012), an aligner tool to align the trimmed reads to the selected reference of the influenza strain (i.e. the inoculum), was then used. Samtools suite (Li et al., 2009) was used to sort, index, and generate depth and coverage statistics for read alignment files. Next, CliqueSNV (Knyazev, 2020) was used to infer the haplotypes and frequencies for all eight gene segments for each sample.

anneng

https://www.sciencedirect.com/science/article/pii/S004268221630037X

anneng

走一遍qap的流程
1.fqc

docker run -v /home/bioinfo/hbv_pipeline/data:/data anneng01:8090/app/fqc fqc qc hbv s1 /data/SRR6378032_1.fastq.gz --r2 /data/SRR6378032_2.fastq.gz -o /data/qc/

2.cutadapt
具体算法见
https://cutadapt.readthedocs.io/en/stable/guide.html?highlight=max-n#dealing-with-n-bases

docker run -v /home/bioinfo/hbv_pipeline/:/workplace pegi3s/cutadapt -q 1 --max-n 0 --minimum-length 10 -o /workplace/data/SRR6378032_1.cleaned.fastq.gz -p /workplace/data/SRR6378032_2.cleaned.fastq.gz /workplace/data/SRR6378032_1.fastq.gz /workplace/data/SRR6378032_2.fastq.gz

-q 按照质量值进行过滤
--max-n 按N碱基数量进行过滤
--minimum-length 按长度进行过滤
-o R1的输出
-p R2的输出

3.qap的环状参考基因组修复过程涉及了几个自己的perl和R脚本算法质量情况不明确我们采用下面的软件来替代：
https://github.com/apeltzer/CircularMapper
java -jar generator-1.93.5.jar CircularGenerator -e 20 -i ../data/demo.fasta -s "AB033556.1"
这个软件有个bug fasta中的序列id不能包括空格有空格的话就找不到这条序列导致没有进行处理处理的算法很简单就是从头部取一定的碱基数量加到尾部
本步骤要做成一个cwl的选择项只有环状基因组需要执行这个步骤
参考序列处理完毕后就可以使用新生成的参考序列进行比对（BWA、Bowtie2等）但是比对完毕后要继续使用另外一个模块(RealignSAMFile.jar)进行重新对齐

4.序列比对
docker run -v /home/bioinfo/hbv_pipeline/data/:/data anneng01:8090/library/angs_bwa:1.0.0 bwa mem -M /data/HBV_C_AB033556_150.fasta /data/SRR6377924_1.fastq.gz /data/SRR6377924_2.fastq.gz -o /data/SRR6377924.sam

//下面的命令不支持bwa mem 的结果 先不执行
java -jar RealignSAMFile.jar -e 500 -i SRR6377924.sam -r HBV_C_AB033556.fasta

过滤没有比对上的序列
docker run -v /home/bioinfo/hbv_pipeline/:/work jweinstk/samtools samtools view -bF 4 /work/mapping/SRR6377924.sam -o /work/mapping/SRR6377924.bam

5.去除PCR重复

docker run -v /home/bioinfo/hbv_pipeline/:/work quay.io/biocontainers/samtools:1.15.1--h1170115_0 samtools collate -o /work/mapping/SRR6377924.namecollate.bam /work/mapping/SRR6377924.bam

docker run -v /home/bioinfo/hbv_pipeline/:/work quay.io/biocontainers/samtools:1.15.1--h1170115_0 samtools fixmate -m /work/mapping/SRR6377924.namecollate.bam  /work/mapping/SRR6377924.fixmate.bam

docker run -v /home/bioinfo/hbv_pipeline/:/work quay.io/biocontainers/samtools:1.15.1--h1170115_0 samtools sort -o /work/mapping/SRR6377924.sorted.bam /work/mapping/SRR6377924.fixmate.bam

docker run -v /home/bioinfo/hbv_pipeline/:/work quay.io/biocontainers/samtools:1.15.1--h1170115_0 samtools markdup -r  /work/mapping/SRR6377924.sorted.bam /work/mapping/SRR6377924.sorted.rmdup.bam

6.call snp

docker run -v /home/bioinfo/hbv_pipeline/:/workplace quay.io/biocontainers/lofreq:broken---2.5.1--py38h1bd3507_2 lofreq call -f /workplace/data/HBV_C_AB033556_150.fasta -o /workplace/calling/SRR6377924.vcf /workplace/mapping/SRR6377924.sorted.rmdup.bam

7.merge R1 R2

docker run -v /home/bioinfo/hbv_pipeline/:/workhome quay.io/biocontainers/samtools:1.15.1--h1170115_0 samtools view -h -o /workhome/mapping/SRR6377924.sorted.rmdup.sam /workhome/mapping/SRR6377924.sorted.rmdup.bam

docker run -v /home/bioinfo/hbv_pipeline/:/workhome quay.io/biocontainers/samtools:1.15.1--h1170115_0 samtools fasta -1 /workhome/merging/SRR6377924_R1.fasta -2 /workhome/merging/SRR6377924_R2.fasta  -0 /dev/null -s /dev/null -n  /workhome/mapping/SRR6377924.sorted.rmdup.sam

docker run -v /home/bioinfo/hbv_pipeline/:/workplace staphb/bbtools bbmerge.sh in1=/workplace/merging/SRR6377924_R1.fasta in2=/workplace/merging/SRR6377924_R2.fasta out=/workplace/merging/SRR6377924_QS.fasta

//======有异常先不用这个软件====

wget https://github.91chi.fun//https://github.com//neufeld/pandaseq/archive/refs/tags/v2.11.tar.gz
tar xvfz v2.11.tar.gz
sudo apt-get install build-essential libtool automake zlib1g-dev libbz2-dev pkg-config
./autogen.sh && ./configure && make && sudo make install

wget https://github.91chi.fun//https://github.com//neufeld/pandaseq-sam/archive/refs/tags/v1.4.tar.gz
./autogen.sh && ./configure && make && sudo make install

//==========================

anneng

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8656693/
The Impact of HBV Quasispecies Features on Immune Status in HBsAg+/HBsAb+ Patients With HBV Genotype C Using Next-Generation Sequencing

anneng

https://virologyj.biomedcentral.com/articles/10.1186/s12985-022-01836-9
Quality evaluation of raw reads was performed with the online tool fastqc (http:// www.bioinformatics.babraham.ac.uk/projects/fastqc/), and the reads having average base calling quality score under 20 were discarded. After quality filtration and adapter removal, paired-end reads were joined with FLASH, v1.2.10 [31]. Merged preS region sequence was genotyped with HBV STAR software as reported previously [32], and corresponding preS regions of 23 reference HBV genomes from the GenBank database were used for genotyping (Accession numbers: X02763, X51970, AF090842, D00329, AB073846, AB602818, X04615, AY123041, AB014381, X65259, M32138, X85254, X75657, AB032431, X69798, AB036910, AF223965, AF160501, AB064310, AF405706, AY090454, AY090457, AY090460). The genotype of each sample was defined as the most frequent one among all 8 types from A to H.

Data preprocessing and predictors
After sequencing the quasispecies, we collected the point mutation data for 457 positions including the positions from 1 to 61 and 2820 to 3215 in and close to the preS region. We counted the frequencies of the nucleotides in each position. To describe the mutation complexity in each position, we transformed the frequency data to Shannon entropy, which is defined as H=−∑ipilogpi, ∑ipi=1 where i∈{A,C,G,T} and pi is its frequency, xlog(x)=0 when x = 0. Entropy of all the 457 nucleotide positions of preS region were used as predictors for HCC diagnosis.