Structural Variant Calling for Nanopore Sequencing (edited)

gene_x 0 like s 398 view s

Tags: variation, pipeline

  1. install mambaforge https://conda-forge.org/miniforge/ (recommended)
    1. #download Mambaforge-24.9.2-0-Linux-x86_64.sh from website
    2. chmod +x Mambaforge-24.9.2-0-Linux-x86_64.sh
    3. ./Mambaforge-24.9.2-0-Linux-x86_64.sh
    5. To activate this environment, use:
    6.     micromamba activate /home/jhuang/mambaforge
    7. Or to execute a single command in this environment, use:
    8.     micromamba run -p /home/jhuang/mambaforge mycommand
    9. installation finished.
    11. Do you wish to update your shell profile to automatically initialize conda?
    12. This will activate conda on startup and change the command prompt when activated.
    13. If you'd prefer that conda's base environment not be activated on startup,
    14.   run the following command when conda is activated:
    16. conda config --set auto_activate_base false
    18. You can undo this by running `conda init --reverse $SHELL`? [yes|no]
    19. [no] >>> yes
    20. no change     /home/jhuang/mambaforge/condabin/conda
    21. no change     /home/jhuang/mambaforge/bin/conda
    22. no change     /home/jhuang/mambaforge/bin/conda-env
    23. no change     /home/jhuang/mambaforge/bin/activate
    24. no change     /home/jhuang/mambaforge/bin/deactivate
    25. no change     /home/jhuang/mambaforge/etc/profile.d/conda.sh
    26. no change     /home/jhuang/mambaforge/etc/fish/conf.d/conda.fish
    27. no change     /home/jhuang/mambaforge/shell/condabin/Conda.psm1
    28. no change     /home/jhuang/mambaforge/shell/condabin/conda-hook.ps1
    29. no change     /home/jhuang/mambaforge/lib/python3.12/site-packages/xontrib/conda.xsh
    30. no change     /home/jhuang/mambaforge/etc/profile.d/conda.csh
    31. modified      /home/jhuang/.bashrc
    32. ==> For changes to take effect, close and re-open your current shell. <==
    33. no change     /home/jhuang/mambaforge/condabin/conda
    34. no change     /home/jhuang/mambaforge/bin/conda
    35. no change     /home/jhuang/mambaforge/bin/conda-env
    36. no change     /home/jhuang/mambaforge/bin/activate
    37. no change     /home/jhuang/mambaforge/bin/deactivate
    38. no change     /home/jhuang/mambaforge/etc/profile.d/conda.sh
    39. no change     /home/jhuang/mambaforge/etc/fish/conf.d/conda.fish
    40. no change     /home/jhuang/mambaforge/shell/condabin/Conda.psm1
    41. no change     /home/jhuang/mambaforge/shell/condabin/conda-hook.ps1
    42. no change     /home/jhuang/mambaforge/lib/python3.12/site-packages/xontrib/conda.xsh
    43. no change     /home/jhuang/mambaforge/etc/profile.d/conda.csh
    44. no change     /home/jhuang/.bashrc
    45. No action taken.
    46. WARNING conda.common.path.windows:_path_to(100): cygpath is not available, fallback to manual path conversion
    47. WARNING conda.common.path.windows:_path_to(100): cygpath is not available, fallback to manual path conversion
    48. Added mamba to /home/jhuang/.bashrc
    49. ==> For changes to take effect, close and re-open your current shell. <==
    50. Thank you for installing Mambaforge!
    52. Close your terminal window and open a new one, or run:
    53. #source ~/mambaforge/bin/activate
    54. conda --version
    55. mamba --version
    57. https://github.com/conda-forge/miniforge/releases
    58. Note
    60.     * After installation, please make sure that you do not have the Anaconda default channels configured.
    61.         conda config --show channels
    62.         conda config --remove channels defaults
    63.         conda config --add channels conda-forge
    64.         conda config --show channels
    65.         conda config --set channel_priority strict
    66.         #conda clean --all
    67.         conda config --remove channels biobakery
    69.     * !!!!Do not install anything into the base environment as this might break your installation. See here for details.!!!!
    71. # --Deprecated method: mamba installing on conda--
    72. #conda install -n base --override-channels -c conda-forge mamba 'python_abi=*=*cp*'
    73. #    * Note that installing mamba into any other environment than base is not supported.
    74. #
    75. #conda activate base
    76. #conda install conda
    77. #conda uninstall mamba
    78. #conda install mamba

2: install required Tools on the mamba env

  1.     * Sniffles2: Detect structural variants, including transposons, from long-read alignments.
  2.     * RepeatModeler2: Identify and classify transposons de novo.
  3.     * RepeatMasker: Annotate known transposable elements using transposon libraries.
  4.     * SVIM: An alternative structural variant caller optimized for long-read sequencing, if needed.
  5.     * SURVIVOR: Consolidate structural variants across samples for comparative analysis.
  7.     mamba deactivate
  8.     # Create a new conda environment
  9.     mamba create -n transposon_long python=3.6 -y
  11.     # Activate the environment
  12.     mamba activate transposon_long
  14.     mamba install -c bioconda sniffles
  15.     mamba install -c bioconda repeatmodeler repeatmasker
  17.     # configure repeatmasker database
  18.     mamba info --envs
  19.     cd /home/jhuang/mambaforge/envs/transposon_long/share/RepeatMasker
  21.     #mamba install python=3.6
  22.     mamba install -c bioconda svim
  23.     mamba install -c bioconda survivor

  1. Test the installed tools

    1. # Check versions
    2. sniffles --version
    3. RepeatModeler -h
    4. RepeatMasker -h
    5. svim --help
    6. SURVIVOR --help
    7. mamba install -c conda-forge perl r

  2. Data Preparation

    1. Raw Signal Data: Nanopore devices generate electrical signal data as DNA passes through the nanopore.
    2. Basecalling: Tools like Guppy or Dorado are used to convert raw signals into nucleotide sequences (FASTQ files).

  3. Preprocessing

    1. Quality Filtering: Remove low-quality reads using tools like Filtlong or NanoFilt.
    2. Adapter Trimming: Identify and remove sequencing adapters with tools like Porechop.

  4. (Optional) Variant Calling for SNP and Indel Detection:

    1. Tools like Medaka, Longshot, or Nanopolish analyze the aligned reads to identify SNPs and small indels.

  5. Alignment and Structural Variant Calling: Tools such as Sniffles or SVIM detect large insertions, deletions, and other structural variants. 使用长读长测序工具如 SVIM 或 Sniffles 检测结构变异(e.g. 散在性重复序列)。

    1.   #NOTE that the ./batch1_depth25/trycycler_WT/reads.fastq and F24A430001437_BACctmoD/BGI_result/Separate/${sample}/1.Cleandata/${sample}.filtered_reads.fq.gz are the same!
    2.   ./4/1.Cleandata/4.filtered_reads.fq.gz
    3.   ./3/1.Cleandata/3.filtered_reads.fq.gz
    4.   ./2/1.Cleandata/2.filtered_reads.fq.gz
    5.   ./8/1.Cleandata/8.filtered_reads.fq.gz
    6.   ./5/1.Cleandata/5.filtered_reads.fq.gz
    7.   ./WT/1.Cleandata/WT.filtered_reads.fq.gz
    8.   ./9/1.Cleandata/9.filtered_reads.fq.gz
    9.   ./10/1.Cleandata/10.filtered_reads.fq.gz
    10.   ./7/1.Cleandata/7.filtered_reads.fq.gz
    11.   ./1/1.Cleandata/1.filtered_reads.fq.gz
    13.   # -- Alignment and Detect structural variants in each sample using SVIM (failed due to the strange output from SVIM!)
    14.   #mamba install -c bioconda ngmlr
    15.   mamba install -c bioconda svim
    16.   for sample in WT 1 2 3 4 5 7 8 9 10; do
    17.       svim reads --aligner ngmlr --nanopore svim_reads_ngmlr_${sample} F24A430001437_BACctmoD/BGI_result/Separate/${sample}/1.Cleandata/${sample}.filtered_reads.fq.gz CP020463.fasta  --cores 10;
    18.   done
    19.   for sample in WT 1 2 3 4 5 7 8 9 10; do
    20.   for sample in 1; do
    21.       #INS,INV,DUP:TANDEM,DUP:INT,BND
    22.       svim reads svim_reads_minimap2_${sample} F24A430001437_BACctmoD/BGI_result/Separate/${sample}/1.Cleandata/${sample}.filtered_reads.fq.gz CP020463.fasta --aligner minimap2 --nanopore --cores 20 --types INS --min_sv_size 100 --sequence_allele --insertion_sequences --read_names;
    23.   done
    24.   #svim alignment svim_alignment_minmap2_1_re 1.sorted.bam CP020463_.fasta --types INS --sequence_alleles --insertion_sequences --read_names
    26.   # -- Results1: Detect structural variants using Minamap2+Sniffles2:
    28.   Minimap2: A commonly used aligner for nanopore sequencing data.
    29.       Align Long Reads to the WT Reference using Minimap2
    31.   for sample in WT 1 2 3 4 5 7 8 9 10; do
    32.       minimap2 --MD -t 60 -ax map-ont CP020463.fasta ./batch1_depth25/trycycler_${sample}/reads.fastq | samtools sort -o ${sample}.sorted.bam
    33.       samtools index ${sample}.sorted.bam
    34.   done
    36.   #sniffles -m WT.sorted.bam -v WT.vcf -s 10 -l 50 -t 60
    37.   #  -s 20: Requires at least 20 reads to support an SV for reporting. --> 10
    38.   #  -l 50: Reports SVs that are at least 50 base pairs long.
    39.   #  -t 4: Uses 4 threads for faster processing. --> 60
    41.   for sample in WT 1 2 3 4 5 7 8 9 10; do
    42.       sniffles -m ${sample}.sorted.bam -v ${sample}.vcf -s 10 -l 50 -t 60
    43.       #QUAL < 20 ||
    44.       bcftools filter -e "INFO/SVTYPE != 'INS'" ${sample}.vcf > ${sample}_filtered.vcf
    45.   done
    47.   # -- Results2: Detect structural variants using NGMLR+Sniffles2
    49.   for sample in WT 1 2 3 4 5 7 8 9 10; do
    50.       #ERROR: No MD string detected! Check bam file! Otherwise generate using e.g. samtools. --> No results!
    51.       #sniffles -m svim_reads_minimap2_${sample}/${sample}.filtered_reads.fq.minimap2.coordsorted.bam -v #sniffles_minimap2_${sample}.vcf -s 10 -l 50 -t 60
    52.       bcftools filter -e "INFO/SVTYPE != 'INS'" sniffles_minimap2_${sample}.vcf > sniffles_minimap2_${sample}_filtered.vcf
    53.       #Using
    54.       sniffles -m svim_reads_ngmlr_${sample}/${sample}.filtered_reads.fq.ngmlr.coordsorted.bam -v sniffles_ngmlr_${sample}.vcf -s 10 -l 50 -t 60
    55.       bcftools filter -e "INFO/SVTYPE != 'INS'" sniffles_ngmlr_${sample}.vcf > sniffles_ngmlr_${sample}_filtered.vcf
    56.   done
    58.   # -- Compare the results1 and results2, and check them each position in IGV!
    60.   #minimap2+sniffles2
    61.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 WT_filtered.vcf | grep -v "##"
    62.   POS
    63.   1855752
    64.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 1_filtered.vcf | grep -v "##"
    65.   POS
    66.   529416
    67.   1855752
    68.   2422820
    69.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 2_filtered.vcf | grep -v "##"
    70.   POS
    71.   529416
    72.   1855752
    73.   2422820
    74.   2424590
    75.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 3_filtered.vcf | grep -v "##"
    76.   POS
    77.   529416
    78.   529416
    79.   529418
    80.   1855752
    81.   2422820
    82.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 4_filtered.vcf | grep -v "##"
    83.   POS
    84.   55682
    85.   529416
    86.   1855752
    87.   2422820
    88.   2424590
    89.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 5_filtered.vcf | grep -v "##"
    90.   POS
    91.   529416
    92.   1855752
    93.   2422820
    94.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 7_filtered.vcf | grep -v "##"
    95.   POS
    96.   518217
    97.   1855752
    98.   2424590
    99.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 8_filtered.vcf | grep -v "##"
    100.   POS
    101.   529416
    102.   1855752
    103.   2422820
    104.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 9_filtered.vcf | grep -v "##"
    105.   POS
    106.   529416
    107.   1855752
    108.   2422820
    109.   2424590
    110.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 10_filtered.vcf | grep -v "##"
    111.   POS
    112.   529416
    113.   1855752
    114.   2422818
    115.   2424590
    117.   #ngmlr+sniffles2
    118.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 sniffles_ngmlr_WT_filtered.vcf | grep -v "##"
    119.   POS
    120.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 sniffles_ngmlr_1_filtered.vcf | grep -v "##"
    121.   POS
    122.   529419
    123.   2422819
    124.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 sniffles_ngmlr_2_filtered.vcf | grep -v "##"
    125.   POS
    126.   529418
    127.   2422820
    128.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 sniffles_ngmlr_3_filtered.vcf | grep -v "##"
    129.   POS
    130.   529418
    131.   2422820
    132.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 sniffles_ngmlr_4_filtered.vcf | grep -v "##"
    133.   POS
    134.   529419
    135.   2422820
    136.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 sniffles_ngmlr_5_filtered.vcf | grep -v "##"
    137.   POS
    138.   529419
    139.   2422820
    140.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 sniffles_ngmlr_7_filtered.vcf | grep -v "##"
    141.   POS
    142.   518219
    143.   2422820
    144.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 sniffles_ngmlr_8_filtered.vcf | grep -v "##"
    145.   POS
    146.   529419
    147.   2422820
    148.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 sniffles_ngmlr_9_filtered.vcf | grep -v "##"
    149.   POS
    150.   529419
    151.   2422820
    152.   (base) jhuang@WS-2290C:/mnt/md1/DATA_md1/Data_Patricia_Transposon$ cut -d$'\t' -f2 sniffles_ngmlr_10_filtered.vcf | grep -v "##"
    153.   POS
    154.   529418
    155.   2422820
    157.   #~/Tools/csv2xls-0.4/csv_to_xls.py sniffles_ngmlr_WT_filtered.vcf sniffles_ngmlr_1_filtered.vcf sniffles_ngmlr_2_filtered.vcf sniffles_ngmlr_3_filtered.vcf sniffles_ngmlr_4_filtered.vcf sniffles_ngmlr_5_filtered.vcf sniffles_ngmlr_7_filtered.vcf sniffles_ngmlr_8_filtered.vcf sniffles_ngmlr_9_filtered.vcf sniffles_ngmlr_10_filtered.vcf -d$'\t' -o putative_transposons2.xls
    159.   # -- Filtering low-complexity insertions using RepeatMasker (TODO: how to use RepeatModeler to generate own lib?)
    161.   python vcf_to_fasta.py variants.vcf variants.fasta
    162.   #python filter_low_complexity.py variants.fasta filtered_variants.fasta retained_variants.fasta
    163.   #Using RepeatMasker to filter the low-complexity fasta, the used h5 lib is
    164.   /home/jhuang/mambaforge/envs/transposon_long/share/RepeatMasker/Libraries/Dfam.h5    #1.9G
    165.   python /home/jhuang/mambaforge/envs/transposon_long/share/RepeatMasker/famdb.py -i /home/jhuang/mambaforge/envs/transposon_long/share/RepeatMasker/Libraries/Dfam.h5 names 'bacteria' | head
    166.   Exact Matches
    167.   =============
    168.   2 bacteria (blast name), Bacteria <bacteria> (scientific name), eubacteria (genbank common name), Monera <bacteria> (in-part), Procaryotae <bacteria> (in-part), Prokaryota <bacteria> (in-part), Prokaryotae <bacteria> (in-part), prokaryote <bacteria> (in-part), prokaryotes <bacteria> (in-part)
    169.   Non-exact Matches
    170.   =================
    171.   1783272 Terrabacteria group (scientific name)
    172.   91061 Bacilli (scientific name), Bacilli Ludwig et al. 2010 (authority), Bacillus/Lactobacillus/Streptococcus group (synonym), Firmibacteria (synonym), Firmibacteria Murray 1988 (authority)
    173.   1239 Bacillaeota (synonym), Bacillaeota Oren et al. 2015 (authority), Bacillota (synonym), Bacillus/Clostridium group (synonym), clostridial firmicutes (synonym), Clostridium group firmicutes (synonym), Firmacutes (synonym), firmicutes (blast name), Firmicutes (scientific name), Firmicutes corrig. Gibbons and Murray 1978 (authority), Low G+C firmicutes (synonym), low G+C Gram-positive bacteria (common name), low GC Gram+ (common name)
    175.   Summary of Classes within Firmicutes:
    176.     * Bacilli (includes many common pathogenic and non-pathogenic Gram-positive bacteria, taxid=91061)
    177.         * Bacillus (e.g., Bacillus subtilis, Bacillus anthracis)
    178.         * Staphylococcus (e.g., Staphylococcus aureus, Staphylococcus epidermidis)
    179.         * Streptococcus (e.g., Streptococcus pneumoniae, Streptococcus pyogenes)
    180.         * Listeria (e.g., Listeria monocytogenes)
    181.     * Clostridia (includes many anaerobic species like Clostridium and Clostridioides)
    182.     * Erysipelotrichia (intestinal bacteria, some pathogenic)
    183.     * Tissierellia (less-studied, veterinary relevance)
    184.     * Mollicutes (cell wall-less, includes Mycoplasma species)
    185.     * Negativicutes (includes some Gram-negative, anaerobic species)
    187.   RepeatMasker -species Bacilli -pa 4 -xsmall variants.fasta
    188.   python extract_unmasked_seq.py variants.fasta.masked unmasked_variants.fasta
    190.   #bcftools filter -i 'QUAL>30 && INFO/SVLEN>100' variants.vcf -o filtered.vcf
    191.   #
    192.   #bcftools view -i 'SVTYPE="INS"' variants.vcf | bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO\n' > insertions.txt
    193.   #mamba install -c bioconda vcf2fasta
    194.   #vcf2fasta variants.vcf -o insertions.fasta
    195.   #grep "SEQS" variants.vcf | awk '{ print $1, $2, $4, $5, $8 }' > insertions.txt
    196.   #python3 filtering_low_complexity.py
    197.   #
    198.   #vcftools --vcf input.vcf --recode --out filtered_output --minSVLEN 100
    199.   #bcftools filter -e 'INFO/SEQS ~ "^(G+|C+|T+|A+){4,}"' variants.vcf -o filtered.vcf
    201.   # -- calculate the percentage of reads
    203.   To calculate the percentage of reads that contain the insertion from the VCF entry, use the INFO and FORMAT fields provided in the VCF record.
    204.   Step 1: Extract Relevant Information
    206.   In the provided VCF entry:
    208.       RE (Reads Evidence): 733 – the total number of reads supporting the insertion.
    209.       GT (Genotype): 1/1 – this indicates a homozygous insertion, meaning all reads covering this region are expected to have the insertion.
    210.       AF (Allele Frequency): 1 – a 100% allele frequency, indicating that every read in this sample supports the insertion.
    211.       DR (Depth Reference): 0 – the number of reads supporting the reference allele.
    212.       DV (Depth Variant): 733 – the number of reads supporting the variant allele (insertion).
    214.   Step 2: Calculate Percentage of Reads Supporting the Insertion
    216.   Using the formula:
    217.   Percentage of reads with insertion=(DVDR+DV)×100
    218.   Percentage of reads with insertion=(DR+DVDV​)×100
    220.   Substitute the values:
    221.   Percentage=(7330+733)×100=100%
    222.   Percentage=(0+733733​)×100=100%
    223.   Conclusion
    225.   Based on the VCF record, 100% of the reads support the insertion, indicating that the insertion is fully present in the sample (homozygous insertion). This is consistent with the AF=1 and GT=1/1 fields.

  6. (failed) using own scripts direct analyze the bam-file via cigarString (failed due to too many short insertions!)

    1. transposons.fasta is a file containing the transposon sequences in FASTA format.
    2. python your_script.py input.bam reference.fasta transposons.fasta
    3. #Transposon_Sequence    Insertion_Frequency
    4. #Tn5                    10
    5. #Tn10                   5
    6. #Unknown                3
    8. python putative_transposons_with_counts.py mapping_WT.sorted.bam CP020463.fasta
    10. rule trim_short_reads:
    11.     input:
    12.         "/data/short-reads.fq.gz"
    13.     output:
    14.         "/data/trimmed-short-reads.fasta"
    15.     shell:
    16.         "python3 trim_by_tag_length.py /data/short-reads.fq.gz 10 > /data/trimmed-short-reads.fasta"
    18. rule trim_long_reads:
    19.     input:
    20.         "/data/long-reads.fq.gz"
    21.     output:
    22.         "/data/trimmed-long-reads.fasta"
    23.     shell:
    24.         "python3 trim_by_tag_length.py /data/long-reads.fq.gz 92 > /data/trimmed-long-reads.fasta"
    26. rule install_bwa:
    27.     output:
    28.         "bwa-mem2-2.0pre2_x64-linux/bwa-mem2"
    29.     shell:
    30.         "curl -L https://github.com/bwa-mem2/bwa-mem2/releases/download/v2.0pre2/bwa-mem2-2.0pre2_x64-linux.tar.bz2 | tar jxf -"
    32. rule map_short_reads:
    33.     input:
    34.         "bwa-mem2-2.0pre2_x64-linux/bwa-mem2",
    35.         "/data/reference.fasta",
    36.         "/data/trimmed-short-reads.fasta"
    37.     output:
    38.         "/data/mapping.sam"
    39.     shell:
    40.         """
    41.         bwa-mem2-2.0pre2_x64-linux/bwa-mem2 index /data/reference.fasta
    42.         bwa-mem2-2.0pre2_x64-linux/bwa-mem2 mem /data/reference.fasta /data/trimmed-short-reads.fasta > /data/mapping.sam
    43.         """
    45. rule map_long_reads:
    46.     input:
    47.         "/data/reference.fasta",
    48.         "/data/trimmed-long-reads.fasta"
    49.     output:
    50.         "/data/mapping.bam"
    51.     conda:
    52.         "minimap2.yml"
    53.     shell:
    54.         """
    55.         minimap2 -x map-ont -d reference /data/reference.fasta > /dev/null 2>&1
    56.         minimap2 -c -a -o /data/mapping.nonunique.sam -N 1 -x map-ont reference /data/trimmed-long-reads.fasta
    57.         samtools view -bq 1 /data/mapping.nonunique.sam > /data/mapping.bam
    58.         """
    60. rule convert_sam_to_bam:
    61.     input:
    62.         "/data/mapping.sam"
    63.     output:
    64.         "/data/mapping.bam",
    65.     conda:
    66.         "samtools.yml"
    67.     shell:
    68.         "samtools view -S -b /data/mapping.sam > /data/mapping.bam"
    70. rule get_unmapped_reads:
    71.     input:
    72.         "/data/mapping.bam"
    73.     output:
    74.         "/data/mapping.sorted.bam"
    75.     conda:
    76.         "samtools.yml"
    77.     shell:
    78.         """
    79. #        samtools view -f 4 /data/mapping.bam > /data/unmapped.sam
    80. #        samtools view -S -b /data/unmapped.sam > /data/unmapped.bam
    81. #        samtools bam2fq /data/unmapped.bam | seqtk seq -A - > /data/unmapped.fa
    82.         samtools sort  /data/mapping.bam -o /data/mapping.sorted.bam
    83.         samtools index /data/mapping.sorted.bam
    84.         """
    86. rule create_insertion_plot:
    87.     input:
    88.         "/data/mapping.sorted.bam"
    89.     output:
    90.         "/data/summary-stats.tsv"
    91.     shell:
    92.         """
    93.         python3 ~/Scripts/sam_to_insert_plot.py /data/mapping.sorted.bam /data/reference.fasta > /data/summary-stats.tsv
    94.         """

  7. Polishing of assembly: Use tools like Medaka to refine variant calls by leveraging consensus sequences derived from nanopore data.

    1.   mamba install -c bioconda medaka
    2.   medaka-consensus -i aligned_reads.bam -r reference.fasta -o polished_output -t 4

  8. Compare Insertions Across Samples

    1. Merge Variants Across Samples: Use SURVIVOR to merge and compare the detected insertions in all samples against the WT:
    3. SURVIVOR merge input_vcfs.txt 1000 1 1 1 0 30 merged.vcf
    5.     Input: List of VCF files from Sniffles2.
    6.     Output: A consolidated VCF file with shared and unique variants.
    8. Filter WT Insertions:
    10.     Identify transposons present only in samples 1–9 by subtracting WT variants using bcftools:
    12.         bcftools isec WT.vcf merged.vcf -p comparison_results

  9. Validate and Visualize

    1. Visualize with IGV: Use IGV to inspect insertion sites in the alignment and confirm quality.
    3. igv.sh
    5. Validate Findings:
    6.     Perform PCR or additional sequencing for key transposon insertion sites to confirm results.

  10. Alternatives to TEPID for Long-Read Data

    1. If you’re looking for transposon-specific tools for long reads:
    3.     REPET: A robust transposon annotation tool compatible with assembled genomes.
    4.     EDTA (Extensive de novo TE Annotator):
    5.         A pipeline to identify, classify, and annotate transposons.
    6.         Works directly on your assembled genomes.
    8.         perl EDTA.pl --genome WT.fasta --type all

  11. The WT.vcf file in the pipeline is generated by detecting structural variants (SVs) in the wild-type (WT) genome aligned against itself or using it as a baseline reference. Here’s how you can generate the WT.vcf:

    1. Steps to Generate WT.vcf
    2. 1. Align WT Reads to the WT Reference Genome
    4. The goal here is to create an alignment of the WT sequencing data to the WT reference genome to detect any self-contained structural variations, such as native insertions, deletions, or duplications.
    6. Command using Minimap2:
    8. minimap2 -ax map-ont WT.fasta WT_reads.fastq | samtools sort -o WT.sorted.bam
    10. Index the BAM file:
    12. samtools index WT.sorted.bam
    14. 2. Detect Structural Variants with Sniffles2
    16. Run Sniffles2 on the WT alignment to call structural variants:
    18. sniffles --input WT.sorted.bam --vcf WT.vcf
    20. This step identifies:
    22.     Native transposons and insertions present in the WT genome.
    23.     Other structural variants that are part of the reference genome or sequencing artifacts.
    25. Key parameters to consider:
    27.     --min_support: Adjust based on your WT sequencing coverage.
    28.     --max_distance: Define proximity for merging variants.
    29.     --min_length: Set a minimum SV size (e.g., >50 bp for transposons).

  12. Clean and Filter the WT.vcf, Variant Filtering: Remove low-confidence variants based on read depth, quality scores, or allele frequency.

    1. To ensure the WT.vcf only includes relevant transposons or SVs:
    3.     Use bcftools or similar tools to filter out low-confidence variants:
    5.     bcftools filter -e "QUAL < 20 || INFO/SVTYPE != 'INS'" WT.vcf > WT_filtered.vcf
    6.     bcftools filter -e "QUAL < 1 || INFO/SVTYPE != 'INS'" 1_.vcf > 1_filtered_.vcf

  13. NOTE that in this pipeline, the WT.fasta (reference genome) is typically a high-quality genome sequence from a database or a well-annotated version of your species' genome. It is not assembled from the WT.fastq sequencing reads in this context. Here's why:

    1. Why Use a Reference Genome (WT.fasta) from a Database?
    3.     Higher Quality and Completeness:
    4.         Database references (e.g., NCBI, Ensembl) are typically well-assembled, highly polished, and annotated. They serve as a reliable baseline for variant detection.
    6.     Consistency:
    7.         Using a standard reference ensures consistent comparisons across your WT and samples (1–9). Variants detected will be relative to this reference, not influenced by possible assembly errors.
    9.     Saves Time:
    10.         Assembling a reference genome from WT reads requires significant computational effort. Using an existing reference streamlines the analysis.
    12. Alternative: Assembling WT from FASTQ
    14. If you don’t have a high-quality reference genome (WT.fasta) and must rely on your WT FASTQ reads:
    16.     Assemble the genome from your WT.fastq:
    17.         Use long-read assemblers like Flye, Canu, or Shasta to create a draft genome.
    19.     flye --nano-raw WT.fastq --out-dir WT_assembly --genome-size <size>
    21.     Polish the assembly using tools like Racon (with the same reads) or Medaka for higher accuracy.
    22.     Use the assembled and polished genome as your WT.fasta reference for further steps.
    24. Key Takeaways:
    26.     If you have access to a reliable, high-quality reference genome, use it as the WT.fasta.
    27.     Only assemble WT.fasta from raw reads (WT.fastq) if no database reference is available for your organism.

  14. Annotate Transposable Elements: Tools like ANNOVAR or SnpEff provide functional insights into the detected variants.

    1. # -- (successful!) MANUALLY Search for all found insertion sequences at https://tncentral.ncc.unesp.br/blast/ !
    2. # Or use the program available at https://github.com/danillo-alvarenga/tncomp_finder if there are numerous matches.
    3. #https://tncentral.ncc.unesp.br/report/te/Tn551-Y13600.1
    5. # -- (failed!) try TEPID for annotation
    6. mamba install tepid=0.10 -c bioconda
    7. #(tepid_env)
    8. for sample in WT 1 2 3 4 5 7 8 9 10; do
    9.     tepid-map-se -x CP020463 -p 10 -n ${sample}_tepid -q  ../batch1_depth25/trycycler_${sample}/reads.fastq;
    10.     tepid-discover -k -i -p 10 -n ${sample}_tepid -c ${sample}_tepid --se;
    11. done
    13. tepid-discover -k -i -p 10 -n 1_tepid -c 1.sorted.bam --se;
    15. tepid-refine [-h] [--version] [-k] [-i INSERTIONS] [-d DELETIONS]
    16.             [-p PROC] -t TE -n NAME -c CONC -s SPLIT -a AL
    18. # -- (failed!) try EDTA for annotation
    19. https://github.com/oushujun/EDTA
    20. (transposon_long) mamba install -c conda-forge -c bioconda edta
    21. mamba install bioconda::rmblast  # cd RepeatMasker; ./configure
    22. EDTA.pl --genome CP020463.fasta --species others --threads 40
    24. For general-purpose TE annotation: EDTA, RepeatMasker, or RepeatScout are your best options.
    25. For de novo repeat identification: RepeatScout is highly effective.
    26. For LTR retrotransposons: Use LTR_retriever.
    27. For bacterial-specific annotations: Transposome, TEfinder, and ISfinder can be useful.

  15. Validation: Cross-validate with short-read sequencing data if available.

like unlike

点赞本文的读者

还没有人对此文章表态

本文有评论

没有评论

看文章，发评论，不要沉默

最受欢迎文章

  1. Why Do Significant Gene Lists Change After Adding Additional Conditions in Differential Gene Expression Analysis?
  2. Motif Discovery in Biological Sequences: A Comparison of MEME and HOMER
  3. Calling peaks using findPeaks of HOMER
  4. Updating Human Gene Identifiers using Ensembl BioMart: A Step-by-Step Guide
  5. Setup conda environments
  6. Should the inputs for GSVA be normalized or raw?
  7. pheatmap vs heatmap.2
  8. File format for single channel analysis of Agilent microarray data with Limma?
  9. 洛那法尼(lonafarnib)：从抗癌症到抗病毒的多功能药物
  10. RNAseq running with umi_tools

最新文章

最多评论文章

推荐相似文章

Call and merge SNP and Indel results from using two variant calling methods

Scaffolding and finishing an assembly with a reference genome

Snakefile

Analysis of the RNA binding protein motifs for RNA Seq

© 2023 XGenes.com Impressum