Variant calling for Data_Huang_Human_herpesvirus_3 using snippy+spandx+viralngs

Leave a reply

  1. Input files

    mkdir raw_data; cd raw_data;

# Note that the names must be ending with fastq.gz
ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb180/01_VZV_20S_S1_R1_001.fastq.gz VZV_20S_R1.fastq.gz
ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb180/01_VZV_20S_S1_R2_001.fastq.gz VZV_20S_R2.fastq.gz
ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb181/02_VZV_20c_S2_R1_001.fastq.gz VZV_20c_R1.fastq.gz
ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb181/02_VZV_20c_S2_R2_001.fastq.gz VZV_20c_R2.fastq.gz

ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb182/03_VZV_60S_S3_R1_001.fastq.gz VZV_60S_R1.fastq.gz
ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb182/03_VZV_60S_S3_R2_001.fastq.gz VZV_60S_R2.fastq.gz
ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb183/04_VZV_60c_S4_R1_001.fastq.gz VZV_60c_R1.fastq.gz
ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb183/04_VZV_60c_S4_R2_001.fastq.gz VZV_60c_R2.fastq.gz

ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb184/05_VZV_1451S_S5_R1_001.fastq.gz VZV_1451S_R1.fastq.gz
ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb184/05_VZV_1451S_S5_R2_001.fastq.gz VZV_1451S_R2.fastq.gz
ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb185/06_Pcc1_1451_S6_R1_001.fastq.gz Pcc1_1451_R1.fastq.gz
ln -s ../VZV/241121_VH00358_117_AAGFF7FM5_Dongdong/wb185/06_Pcc1_1451_S6_R2_001.fastq.gz Pcc1_1451_R2.fastq.gz

ln -s ../VZV/2024_081_wb_dongdong/wb190/PCC1_VZV_20_1_S36_R1_001.fastq.gz PCC1_VZV_20_1_R1.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb190/PCC1_VZV_20_1_S36_R2_001.fastq.gz PCC1_VZV_20_1_R2.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb191/PCC1_VZV_20_2_S37_R1_001.fastq.gz PCC1_VZV_20_2_R1.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb191/PCC1_VZV_20_2_S37_R2_001.fastq.gz PCC1_VZV_20_2_R2.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb192/PCC1_VZV_20_5_S38_R1_001.fastq.gz PCC1_VZV_20_5_R1.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb192/PCC1_VZV_20_5_S38_R2_001.fastq.gz PCC1_VZV_20_5_R2.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb193/PCC1_VZV_60_1_S39_R1_001.fastq.gz PCC1_VZV_60_1_R1.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb193/PCC1_VZV_60_1_S39_R2_001.fastq.gz PCC1_VZV_60_1_R2.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb194/PCC1_VZV_60_4_S40_R1_001.fastq.gz PCC1_VZV_60_4_R1.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb194/PCC1_VZV_60_4_S40_R2_001.fastq.gz PCC1_VZV_60_4_R2.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb195/PCC1_VZV_60_6_S41_R1_001.fastq.gz PCC1_VZV_60_6_R1.fastq.gz
ln -s ../VZV/2024_081_wb_dongdong/wb195/PCC1_VZV_60_6_S41_R2_001.fastq.gz PCC1_VZV_60_6_R2.fastq.gz

  2. Call variant calling using snippy

    ln -s ~/Tools/bacto/db/ .;
ln -s ~/Tools/bacto/envs/ .;
ln -s ~/Tools/bacto/local/ .;
cp ~/Tools/bacto/Snakefile .;
cp ~/Tools/bacto/bacto-0.1.json .;
cp ~/Tools/bacto/cluster.json .;

#download CU459141.gb from GenBank
mv ~/Downloads/sequence\(1\).gb db/NC_001348.gb = X04370.1
mv ~/Downloads/sequence\(2\).gb db/AB097932.gb
#X04370
#setting the following in bacto-0.1.json
    "fastqc": false,
    "taxonomic_classifier": false,
    "assembly": true,
    "typing_ariba": false,
    "typing_mlst": true,
    "pangenome": true,
    "variants_calling": true,
    "phylogeny_fasttree": true,
    "phylogeny_raxml": true,
    "recombination": false, (due to gubbins-error set false)
        "genus": "Varicella-zoster virus",
        "kingdom": "Viruses",
        "species": "Varicella-zoster virus"(in both prokka and mykrobe)
    "reference": "db/NC_001348.gb"
conda activate bengal3_ac3
(bengal3_ac3) /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds

  3. Summarize all SNPs and Indels from the snippy result directory.

    #Output: snippy/summary_snps_indels.csv
# IMPORTANT_ADAPT the array isolates = ["AYE-S", "AYE-Q", "AYE-WT on Tig4", "AYE-craA on Tig4", "AYE-craA-1 on Cm200", "AYE-craA-2 on Cm200"]
python3 ~/Scripts/summarize_snippy_res.py snippy
cd snippy
grep -v "None,,,,,,None,None" summary_snps_indels.csv > summary_snps_indels_.csv

  4. Using spandx calling variants (almost the same results to the one from viral-ngs!)

    mamba activate /home/jhuang/miniconda3/envs/spandx
mkdir ~/miniconda3/envs/spandx/share/snpeff-5.1-2/data/NC_001348
cp NC_001348.gb  ~/miniconda3/envs/spandx/share/snpeff-5.1-2/data/NC_001348/genes.gbk
vim ~/miniconda3/envs/spandx/share/snpeff-5.1-2/snpEff.config
/home/jhuang/miniconda3/envs/spandx/bin/snpEff build NC_001348    #-d
~/Scripts/genbank2fasta.py NC_001348.gb
mv NC_001348.gb_converted.fna NC_001348.fasta    #rename "NC_001348.1 xxxxx" to "NC_001348" in the fasta-file
ln -s /home/jhuang/Tools/spandx/ spandx
(spandx) nextflow run spandx/main.nf --fastq "trimmed/*_P_{1,2}.fastq" --ref NC_001348.fasta --annotation --database NC_001348 -resume

# Rerun SNP_matrix.sh due to the error ERROR_CHROMOSOME_NOT_FOUND in the variants annotation
cd Outputs/Master_vcf
(spandx) cp -r ../../snippy/VZV_20S/reference .
(spandx) cp ../../spandx/bin/SNP_matrix.sh ./
#Note that ${variant_genome_path}=NC_001348 in the following command, but it was not used after command replacement.
#Adapt "snpEff eff -no-downstream -no-intergenic -ud 100 -formatEff -v ${variant_genome_path} out.vcf > out.annotated.vcf" to
"/home/jhuang/miniconda3/envs/bengal3_ac3/bin/snpEff eff -no-downstream -no-intergenic -ud 100 -formatEff -c reference/snpeff.config -dataDir . ref out.vcf > out.annotated.vcf" in SNP_matrix.sh
(spandx) bash SNP_matrix.sh NC_001348 .

  5. Calling inter-host variants by merging the results from snippy+spandx (Manually!)

    # Inter-host variants（宿主间变异）:一种病毒在两个人之间有不同的基因变异，这些变异可能与宿主的免疫反应、疾病表现或病毒传播的方式相关。
cp All_SNPs_indels_annotated.txt All_SNPs_indels_annotated_backup.txt
vim All_SNPs_indels_annotated.txt

  6. Calling intra-host variants using viral-ngs (http://xgenes.com/article/article-content/347/variant-calling-for-herpes-simplex-virus-1-from-patient-sample-using-capture-probe-sequencing/)

    # Intra-host variants（宿主内变异）：同一个人感染了某种病毒，但在其体内的不同细胞或器官中可能存在多个不同的病毒变异株。
mamba activate /home/jhuang/miniconda3/envs/viral-ngs4

mkdir viralngs

ln -s ~/Tools/viral-ngs/Snakefile Snakefile
ln -s  ~/Tools/viral-ngs/bin bin
cp  ~/Tools/viral-ngs/refsel.acids refsel.acids
cp  ~/Tools/viral-ngs/lastal.acids lastal.acids
cp  ~/Tools/viral-ngs/config.yaml config.yaml
cp  ~/Tools/viral-ngs/samples-runs.txt samples-runs.txt
cp  ~/Tools/viral-ngs/samples-depletion.txt samples-depletion.txt
cp  ~/Tools/viral-ngs/samples-metagenomics.txt samples-metagenomics.txt
cp  ~/Tools/viral-ngs/samples-assembly.txt samples-assembly.txt
cp  ~/Tools/viral-ngs/samples-assembly-failures.txt samples-assembly-failures.txt
mkdir data
cd data
mkdir 00_raw
cd ../..

mkdir bams
ref_fa="NC_001348.fasta";
for sample in VZV_20S VZV_20c VZV_60S VZV_60c PCC1_VZV_20_1 PCC1_VZV_20_2 PCC1_VZV_20_5 PCC1_VZV_60_1 PCC1_VZV_60_4 PCC1_VZV_60_6; do
    bwa index ${ref_fa}; \
    bwa mem -M -t 16 ${ref_fa} trimmed/${sample}_trimmed_P_1.fastq trimmed/${sample}_trimmed_P_2.fastq | samtools view -bS - > bams/${sample}_genome_alignment.bam; \
done
for sample in VZV_20S VZV_20c VZV_60S VZV_60c PCC1_VZV_20_1 PCC1_VZV_20_2 PCC1_VZV_20_5 PCC1_VZV_60_1 PCC1_VZV_60_4 PCC1_VZV_60_6; do
    picard AddOrReplaceReadGroups I=bams/${sample}_genome_alignment.bam O=viralngs/data/00_raw/${sample}.bam SORT_ORDER=coordinate CREATE_INDEX=true RGPL=illumina RGID=$sample RGSM=$sample RGLB=standard RGPU=$sample VALIDATION_STRINGENCY=LENIENT; \
done

cd viralngs
(viral-ngs4) snakemake --printshellcmds --cores 80

# -- DEBUG: If the env disappeared, reinstall the env viral-ngs4 --

# -- Running time hints --
#Note that novoalign is not installed. The used Novoalign path: /home/jhuang/Tools/novocraft_v3/novoalign; the used gatk: /usr/local/bin/gatk using /home/jhuang/Tools/GenomeAnalysisTK-3.6/GenomeAnalysisTK.jar.

#Samtools path:  #Why, the samtools in the env is v1.6?
#Novoalign path: /home/jhuang/Tools/novocraft_v3/novoalign
#GATK path: /usr/local/bin/gatk  # jar_file in the file: jar_file = '/home/jhuang/Tools/GenomeAnalysisTK-3.6/GenomeAnalysisTK.jar'

# -- in config.yaml --
#GATK_PATH: "/home/jhuang/Tools/GenomeAnalysisTK-3.6"
#NOVOALIGN_PATH: "/home/jhuang/Tools/novocraft_v3"

mamba create -n viral-ngs4 python=3.6
mamba activate viral-ngs4
mamba install blast=2.6.0 bmtagger biopython pysam pyyaml picard mvicuna pybedtools fastqc matplotlib spades last=876 -c conda-forge -c bioconda
#mafft=7.221 --> mafft since └─ mafft 7.221**  is not installable because it conflicts with any installable versions previously reported.
mamba install cd-hit cd-hit-auxtools diamond gap2seq=2.1 mafft mummer4 muscle=3.8 parallel pigz prinseq samtools=1.6 tbl2asn trimmomatic trinity unzip vphaser2 bedtools -c r -c defaults -c conda-forge -c bioconda
mamba install bwa
mamba install vphaser2=2.0

# Sovle confilict between bowtie, bowtie2 and snpeff
mamba remove bowtie
mamba install bowtie2
mamba remove snpeff
mamba install snpeff=4.1l
#which snpEff

mamba install gatk=3.6
#DEBUG if FileNotFoundError: [Errno 2] No such file or directory: '/usr/local/bin/gatk': '/usr/local/bin/gatk'
#IMPORTANT_UPDATE jar_file in the file /home/jhuang/mambaforge/envs/viral-ngs4/bin/gatk3 with "/home/jhuang/Tools/GenomeAnalysisTK-3.6/GenomeAnalysisTK.jar"
#IMPORTANT_REPLACE "sudo cp /home/jhuang/mambaforge/envs/viral-ngs4/bin/gatk3 /usr/local/bin/gatk"
#IMPORTANT_SET /home/jhuang/Tools/GenomeAnalysisTK-3.6 as GATK_PATH in config.yaml
#IMPORTANT_CHECK if it works
#        java -jar /home/jhuang/Tools/GenomeAnalysisTK-3.6/GenomeAnalysisTK.jar -T RealignerTargetCreator --help
#        /usr/local/bin/gatk -T RealignerTargetCreator --help
#IMPORTANT_NOTE that the env viral-ngs4 cannot logined from the base env due to the python3-conflict!

  7. Merge intra- and inter-host variants, comparing the variants to the alignments of the assemblies to confirm its correctness.

    cat NC_001348.fasta viralngs/data/02_assembly/VZV_20S.fasta viralngs/data/02_assembly/VZV_60S.fasta > aligned_1.fasta
mafft --clustalout aligned_1.fasta > aligned_1.aln
#~/Scripts/convert_fasta_to_clustal.py aligned_1.fasta_orig aligned_1.aln
~/Scripts/convert_clustal_to_clustal.py aligned_1.aln aligned_1_.aln
#manully delete the postion with all or '-' in aligned_1_.aln
~/Scripts/check_sequence_differences.py aligned_1_.aln
~/Scripts/check_sequence_differences.py aligned_1_.aln > aligned_1.res
grep -v " = n" aligned_1.res > aligned_1_.res

cat NC_001348.fasta viralngs/tmp/02_assembly/VZV_20S.assembly4-refined.fasta viralngs/tmp/02_assembly/VZV_60S.assembly4-refined.fasta > aligned_1.fasta
mafft --clustalout aligned_1.fasta > aligned_1.aln
~/Scripts/convert_clustal_to_clustal.py aligned_1.aln aligned_1_.aln
~/Scripts/check_sequence_differences.py aligned_1_.aln > aligned_1.res
grep -v " = n" aligned_1.res > aligned_1_.res

#Differences found at the following positions (150):
Position 8956: OP297860.1 = A, HSV1_S1-1 = A, HSV-Klinik_S2-1 = G
Position 8991: OP297860.1 = A, HSV1_S1-1 = A, HSV-Klinik_S2-1 = C
Position 8992: OP297860.1 = T, HSV1_S1-1 = C, HSV-Klinik_S2-1 = C
Position 8995: OP297860.1 = T, HSV1_S1-1 = T, HSV-Klinik_S2-1 = C
Position 9190: OP297860.1 = T, HSV1_S1-1 = A, HSV-Klinik_S2-1 = T
* Position 13659: OP297860.1 = G, HSV1_S1-1 = T, HSV-Klinik_S2-1 = G
* Position 47969: OP297860.1 = C, HSV1_S1-1 = T, HSV-Klinik_S2-1 = C
* Position 53691: OP297860.1 = G, HSV1_S1-1 = T, HSV-Klinik_S2-1 = G
* Position 55501: OP297860.1 = T, HSV1_S1-1 = C, HSV-Klinik_S2-1 = C
* Position 63248: OP297860.1 = G, HSV1_S1-1 = T, HSV-Klinik_S2-1 = G
Position 63799: OP297860.1 = T, HSV1_S1-1 = C, HSV-Klinik_S2-1 = T
* Position 64328: OP297860.1 = C, HSV1_S1-1 = A, HSV-Klinik_S2-1 = C
Position 65179: OP297860.1 = T, HSV1_S1-1 = T, HSV-Klinik_S2-1 = C
* Position 65225: OP297860.1 = G, HSV1_S1-1 = G, HSV-Klinik_S2-1 = A
* Position 95302: OP297860.1 = C, HSV1_S1-1 = A, HSV-Klinik_S2-1 = C

gunzip isnvs.annot.txt.gz
~/Scripts/filter_isnv.py isnvs.annot.txt 0.05
cut -d$'\t' filtered_isnvs.annot.txt -f1-7
chr     pos     sample  patient time    alleles iSNV_freq
OP297860        13203   HSV1_S1 HSV1_S1         T,C,A   1.0
OP297860        13203   HSV-Klinik_S2   HSV-Klinik_S2           T,C,A   1.0
OP297860        13522   HSV1_S1 HSV1_S1         G,T     1.0
OP297860        13522   HSV-Klinik_S2   HSV-Klinik_S2           G,T     0.008905554253573941
OP297860        13659   HSV1_S1 HSV1_S1         G,T     1.0
OP297860        13659   HSV-Klinik_S2   HSV-Klinik_S2           G,T     0.008383233532934131

~/Scripts/convert_clustal_to_fasta.py aligned_1_.aln aligned_1.fasta
samtools faidx aligned_1.fasta
samtools faidx aligned_1.fasta OP297860.1 > OP297860.1.fasta
samtools faidx aligned_1.fasta HSV1_S1-1 > HSV1_S1-1.fasta
samtools faidx aligned_1.fasta HSV-Klinik_S2-1 > HSV-Klinik_S2-1.fasta
seqkit seq OP297860.1.fasta -w 70 > OP297860.1_w70.fasta
diff OP297860.1_w70.fasta ../../refsel_db/refsel.fasta

  8. Consensus sequences of each and of all isolates

    cp data/02_assembly/*.fasta ./
for sample in 838_S1 840_S2 820_S3 828_S4 815_S5 834_S6 808_S7 811_S8 837_S9 768_S10 773_S11 767_S12 810_S13 814_S14 10121-16_S15 7510-15_S16 828-17_S17 8806-15_S18 9881-16_S19 8981-14_S20; do
for sample in p953-84660-tsek p938-16972-nra p942-88507-nra p943-98523-nra p944-103323-nra p947-105565-nra p948-112830-nra; do \
mv ${sample}.fasta ${sample}.fa
cat all.fa ${sample}.fa >> all.fa
done
cat RSV_dedup.fa all.fa > RSV_all.fa
mafft --adjustdirection RSV_all.fa > RSV_all.aln
snp-sites RSV_all.aln -o RSV_all_.aln

  9. Download all Human alphaherpesvirus 3 (Varicella-zoster virus) genomes

    Human alphaherpesvirus 3
acronym: HHV-3 VZV
equivalent: Human herpes virus 3

Human alphaherpesvirus 3 (Varicella-zoster virus)
    * Human herpesvirus 3 strain Dumas
    * Human herpesvirus 3 strain Oka vaccine
    * Human herpesvirus 3 VZV-32

#Taxonomy ID: 10335
esearch -db nucleotide -query "txid10335[Organism:exp]" | efetch -format fasta -email j.huang@uke.de > genome_10335_ncbi.fasta
python ~/Scripts/filter_fasta.py genome_10335_ncbi.fasta complete_genome_10335_ncbi.fasta  #2041-->165
# ---- Download related genomes from ENA ----
https://www.ebi.ac.uk/ena/browser/view/10335
#Click "Sequence" and download "Counts" (2003) and "Taxon descendants count" (2005) if there is enough time! Downloading time points is 11.03.2025.
python ~/Scripts/filter_fasta.py  ena_10335_sequence.fasta complete_genome_10335_ena_taxon_descendants_count.fasta  #2005-->153
#python ~/Scripts/filter_fasta.py ena_10335_sequence_Counts.fasta complete_genome_10335_ena_Counts.fasta  #xxx, 5.8G

  10. Run vrap

    #replace --virus to the specific taxonomy (e.g. Acinetobacter baumannii) --> change virus_user_db --> specific_bacteria_user_db
ln -s ~/Tools/vrap/ .
mamba activate /home/jhuang/miniconda3/envs/vrap

#!!!!! TODO: ignore the first parts! only take the virus genome in the vector-part!
vrap/vrap.py  -1 trimmed/VZV_20c_trimmed_P_1.fastq -2 trimmed/VZV_20c_trimmed_P_2.fastq -o vrap_VZV_20c  --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
#-1 trimmed/VZV_20S_trimmed_P_1.fastq -2 trimmed/VZV_20S_trimmed_P_2.fastq
#(vrap) vrap/vrap.py  -1 trimmed/VZV_20S_trimmed_P_1.fastq -2 trimmed/VZV_20S_trimmed_P_2.fastq -o vrap_VZV_20S --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
vrap/vrap.py  -1 trimmed/VZV_60c_trimmed_P_1.fastq -2 trimmed/VZV_60c_trimmed_P_2.fastq -o vrap_VZV_60c  --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
vrap/vrap.py  -1 trimmed/VZV_60S_trimmed_P_1.fastq -2 trimmed/VZV_60S_trimmed_P_2.fastq -o vrap_VZV_60S
vrap/vrap.py  -1 trimmed/VZV_1451S_trimmed_P_1.fastq -2 trimmed/VZV_1451S_trimmed_P_2.fastq -o vrap_VZV_1451S  --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
vrap/vrap.py  -1 trimmed/Pcc1_1451_trimmed_P_1.fastq -2 trimmed/Pcc1_1451_trimmed_P_2.fastq -o vrap_Pcc1_1451  --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
vrap/vrap.py  -1 trimmed/PCC1_VZV_20_1_trimmed_P_1.fastq -2 trimmed/PCC1_VZV_20_1_trimmed_P_2.fastq -o vrap_PCC1_VZV_20_1  --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
vrap/vrap.py  -1 trimmed/PCC1_VZV_20_2_trimmed_P_1.fastq -2 trimmed/PCC1_VZV_20_2_trimmed_P_2.fastq -o vrap_PCC1_VZV_20_2  --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
vrap/vrap.py  -1 trimmed/PCC1_VZV_20_5_trimmed_P_1.fastq -2 trimmed/PCC1_VZV_20_5_trimmed_P_2.fastq -o vrap_PCC1_VZV_20_5  --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
vrap/vrap.py  -1 trimmed/PCC1_VZV_60_1_trimmed_P_1.fastq -2 trimmed/PCC1_VZV_60_1_trimmed_P_2.fastq -o vrap_PCC1_VZV_60_1  --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
vrap/vrap.py  -1 trimmed/PCC1_VZV_60_4_trimmed_P_1.fastq -2 trimmed/PCC1_VZV_60_4_trimmed_P_2.fastq -o vrap_PCC1_VZV_60_4  --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g
vrap/vrap.py  -1 trimmed/PCC1_VZV_60_6_trimmed_P_1.fastq -2 trimmed/PCC1_VZV_60_6_trimmed_P_2.fastq -o vrap_PCC1_VZV_60_6  --bt2idx=/home/jhuang/REFs/genome  --host=/home/jhuang/REFs/genome.fa --virus=/home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/complete_genome_10335_ncbi.fasta --nt=/mnt/nvme1n1p1/blast/nt --nr=/mnt/nvme1n1p1/blast/nr  -t 100 -l 200  -g

http://xgenes.com/article/article-content/365/virus-genome-analysis-pipeline-hybrid-capture-damian-blastn-and-vrap-mapping-for-measles-ma-zhen-sample/ Draw the mapping figures on the reference, consensus reference!

  1. Using the bowtie of vrap to map the reads on ref_genome/reference.fasta (The reference refers to the closest related genome found from the list generated by vrap)

    (vrap) vrap/vrap.py  -1 trimmed/VZV_20S_trimmed_P_1.fastq -2 trimmed/VZV_20S_trimmed_P_2.fastq  -o VZV_20S_on_X04370 --host /home/jhuang/DATA/Data_Huang_Human_herpesvirus_3/X04370.fasta   -t 100 -l 200  -g
cd bowtie
mv mapped mapped.sam
samtools view -S -b mapped.sam > mapped.bam
samtools sort mapped.bam -o mapped_sorted.bam
samtools index mapped_sorted.bam
samtools view -H mapped_sorted.bam
samtools flagstat mapped_sorted.bam

  2. Show the bw on IGV

Leave a Reply

Your email address will not be published. Required fields are marked *