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

gene_x 0 like s 880 view s

Tags: variation, pipeline

  1. call variant calling using snippy
    1. mv snippy_HDRNA_01 snippy_CP133676
    3. mv snippy_HDRNA_03 snippy_CP133677
    5. cp -r snippy_HDRNA_06 snippy_CP133678
    6. mv snippy_HDRNA_06 snippy_CP133679
    8. cp -r snippy_HDRNA_07 snippy_CP133680
    9. mv snippy_HDRNA_07 snippy_CP133681
    11. cp -r snippy_HDRNA_08 snippy_CP133682
    12. mv snippy_HDRNA_08 snippy_CP133683
    14. cp -r snippy_HDRNA_12 snippy_CP133684
    15. cp -r snippy_HDRNA_12 snippy_CP133685
    16. cp -r snippy_HDRNA_12 snippy_CP133686
    17. mv snippy_HDRNA_12 snippy_CP133687
    19. cp -r snippy_HDRNA_16 snippy_CP133688
    20. cp -r snippy_HDRNA_16 snippy_CP133689
    21. cp -r snippy_HDRNA_16 snippy_CP133690
    22. cp -r snippy_HDRNA_16 snippy_CP133691
    23. mv snippy_HDRNA_16 snippy_CP133692
    25. cp -r snippy_HDRNA_17 snippy_CP133693
    26. cp -r snippy_HDRNA_17 snippy_CP133694
    27. mv snippy_HDRNA_17 snippy_CP133695
    29. cp -r snippy_HDRNA_19 snippy_CP133696
    30. cp -r snippy_HDRNA_19 snippy_CP133697
    31. cp -r snippy_HDRNA_19 snippy_CP133698
    32. mv snippy_HDRNA_19 snippy_CP133699
    34. cp -r snippy_HDRNA_20 snippy_CP133700
    35. mv snippy_HDRNA_20 snippy_CP133701
    37.     #git clone https://github.com/huang/bacto
    38.     #mv bacto/* ./
    39.     #rm -rf bacto
    41. for sample in snippy_CP133678 snippy_CP133679 snippy_CP133680 snippy_CP133681 snippy_CP133682 snippy_CP133683 snippy_CP133684 snippy_CP133685 snippy_CP133686 snippy_CP133687 snippy_CP133688 snippy_CP133689 snippy_CP133690 snippy_CP133691 snippy_CP133692 snippy_CP133693 snippy_CP133694 snippy_CP133695 snippy_CP133696 snippy_CP133697 snippy_CP133698 snippy_CP133699 snippy_CP133700 snippy_CP133701; do
    42.     cd ${sample};
    43.     ln -s ~/Tools/bacto/db/ .;
    44.     ln -s ~/Tools/bacto/envs/ .;
    45.     ln -s ~/Tools/bacto/local/ .;
    46.     cp ~/Tools/bacto/Snakefile .;
    47.     cp ~/Tools/bacto/bacto-0.1.json .;
    48.     cp ~/Tools/bacto/cluster.json .;
    49.     cd ..;
    50. done
    51. #prepare raw_data and bacto-0.1.json
    52. #set "reference": "db/CP133***.gb" in bacto-0.1.json
    54. conda activate bengal3_ac3
    55. #cd snippy_CP133676
    56. #/home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    57. #cd snippy_CP133677
    58. #/home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    59. cd snippy_CP133678
    60. (bengal3_ac3) /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    61. cd ../snippy_CP133679
    62. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    63. cd ../snippy_CP133680
    64. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    65. cd ../snippy_CP133681
    66. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    67. cd ../snippy_CP133682
    68. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    70. cd snippy_CP133683
    71. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    72. cd ../snippy_CP133684
    73. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    74. cd ../snippy_CP133685
    75. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    76. cd ../snippy_CP133686
    77. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    78. cd ../snippy_CP133687
    79. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    80. cd ../snippy_CP133688
    81. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    82. cd ../snippy_CP133689
    83. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    84. cd ../snippy_CP133690
    85. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    86. cd ../snippy_CP133691
    87. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    88. cd ../snippy_CP133692
    89. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    90. cd ../snippy_CP133693
    91. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    92. cd ../snippy_CP133694
    93. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    94. cd ../snippy_CP133695
    95. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    96. cd ../snippy_CP133696
    97. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    98. cd ../snippy_CP133697
    99. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    100. cd ../snippy_CP133698
    101. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    102. cd ../snippy_CP133699
    103. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    104. cd ../snippy_CP133700
    105. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    106. cd ../snippy_CP133701
    107. /home/jhuang/miniconda3/envs/snakemake_4_3_1/bin/snakemake --printshellcmds
    108. cd ..

2, using spandx calling variants (almost the same results to the one from viral-ngs!)

  1.     mkdir ~/miniconda3/envs/spandx/share/snpeff-5.1-2/data/PP810610
  2.     cp PP810610.1.gb  ~/miniconda3/envs/spandx/share/snpeff-5.1-2/data/PP810610/genes.gbk
  3.     vim ~/miniconda3/envs/spandx/share/snpeff-5.1-2/snpEff.config
  4.     /home/jhuang/miniconda3/envs/spandx/bin/snpEff build  PP810610     -d
  5.     gzip hCoV229E_Rluc_trimmed_P_1.fastq hCoV229E_Rluc_trimmed_P_2.fastq p10_DMSO_trimmed_P_1.fastq p10_DMSO_trimmed_P_2.fastq p10_K22_trimmed_P_1.fastq p10_K22_trimmed_P_2.fastq p10_K7523_trimmed_P_1.fastq p10_K7523_trimmed_P_2.fastq
  6.     ln -s /home/jhuang/Tools/spandx/ spandx
  7.     (spandx) nextflow run spandx/main.nf --fastq "trimmed/*_P_{1,2}.fastq.gz" --ref PP810610.fasta --annotation --database PP810610 -resume

3, summarize all SNPs and Indels from the snippy result directory.

  1.     #Output: snippy_CP133676/snippy/summary_snps_indels.csv
  2.     python3 summarize_snippy_res.py snippy_CP133676/snippy
  4.     #-- post-processing of Outputs_CP133676 produced by SPANDx (temporary not necessary) --
  5.     #cd Outputs_CP133676/Phylogeny_and_annotation
  6.     #awk '{if($3!=$7) print}' < All_SNPs_indels_annotated.txt > All_SNPs_indels_annotated_.txt
  7.     #cut -d$'\t' -f1-7 All_SNPs_indels_annotated_.txt > f1_7
  8.     #grep -v "/" f1_7 > f1_7_
  9.     #grep -v "\." f1_7_ > f1_7__
  10.     #grep -v "*" f1_7__ > f1_7___
  11.     #grep -v "INDEL" f1_7___ > f1_7_SNPs
  12.     #cut -d$'\t' -f2-2 f1_7_SNPs > f1_7_SNPs_f2
  13.     #grep -wFf <(awk 'NR>0 {print $1}' f1_7_SNPs_f2) All_SNPs_indels_annotated_.txt > All_SNPs_annotated.txt
  14.     #grep -v "SNP" f1_7___ > f1_7_indels
  15.     #cut -d$'\t' -f2-2 f1_7_indels > f1_7_indels_f2
  16.     #grep -wFf <(awk 'NR>0 {print $1}' f1_7_indels_f2) All_SNPs_indels_annotated_.txt > All_indels_annotated.txt

4, merge the following two files summary_snps_indels.csv (192) and All_SNPs_indels_annotated.txt (248) --> merged_variants_CP133676.csv (94)

  1.     python3 merge_snps_indels.py snippy_CP133676/snippy/summary_snps_indels.csv Outputs_CP133676/Phylogeny_and_annotation/All_SNPs_indels_annotated.txt merged_variants_CP133676.csv
  2.     #check if the number of the output file is correct?
  3.     comm -12 <(cut -d, -f2 snippy_CP133676/snippy/summary_snps_indels.csv | sort | uniq) <(cut -f2 Outputs_CP133676/Phylogeny_and_annotation/All_SNPs_indels_annotated.txt | sort | uniq) | wc -l
  4.     comm -12 <(cut -d, -f2 snippy_CP133676/snippy/summary_snps_indels.csv | sort | uniq) <(cut -f2 Outputs_CP133676/Phylogeny_and_annotation/All_SNPs_indels_annotated.txt | sort | uniq)

5, (optional, since it should not happed) filter rows without change in snippy_CP133676/snippy/summary_snps_indels.csv (94)

  1.     awk -F, '
  2.     NR == 1 {print; next}  # Print the header line
  3.     {
  4.         ref = $3;  # Reference base (assuming REF is in the 3rd column)
  5.         same = 1;  # Flag to check if all bases are the same as reference
  6.         samples = "";  # Initialize variable to hold sample data
  7.         for (i = 6; i <= NF - 8; i++) {  # Loop through all sample columns, adjusting for the number of fixed columns
  8.             samples = samples $i " ";  # Collect sample data
  9.             if ($i != ref) {
  10.                 same = 0;
  11.             }
  12.         }
  13.         if (!same) {
  14.             print $0;  # Print the entire line if not all bases are the same as the reference
  15.             #print "Samples: " samples;  # Print all sample data for checking
  16.         }
  17.     }
  18.     ' merged_variants_CP133676.csv > merged_variants_CP133676_.csv
  20.     #Explanation:
  21.     #    -F, sets the field separator to a comma.
  22.     #    NR == 1 {print; next} prints the header line and skips to the next line.
  23.     #    ref = $3 sets the reference base (assumed to be in the 3rd column).
  24.     #    same = 1 initializes a flag to check if all sample bases are the same as the reference.
  25.     #    samples = ""; initializes a string to collect sample data.
  26.     #    for (i = 6; i <= NF - 8; i++) loops through the sample columns. This assumes the first 5 columns are fixed (CHROM, POS, REF, ALT, TYPE), and the last 8 columns are fixed (Effect, Impact, Functional_Class, Codon_change, Protein_and_nucleotide_change, Amino_Acid_Length, Gene_name, Biotype).
  27.     #    samples = samples $i " "; collects sample data.
  28.     #    if ($i != ref) { same = 0; } checks if any sample base is different from the reference base. If found, it sets same to 0.
  29.     #    if (!same) { print $0; print "Samples: " samples; } prints the entire line and the sample data if not all sample bases are the same as the reference.

6, improve the header

  1.     sed -i '1s/_trimmed_P//g' merged_variants_CP133676.csv

7, check the REF and K1 have the same base and delete those records with difference.

  1.     cut -f3 -d',' merged_variants_CP133676.csv > f3
  2.     cut -f6 -d',' merged_variants_CP133676.csv > f6
  3.     diff f3 f6
  4.     awk -F, '$3 == $6 || NR==1' merged_variants_CP133676.csv > filtered_merged_variants_CP133676.csv #(93)
  5.     cut -f3 -d',' filtered_merged_variants_CP133676.csv > f3
  6.     cut -f6 -d',' filtered_merged_variants_CP133676.csv > f6
  7.     diff f3 f6

8, MANUALLY REMOVE the column f6 in filtered_merged_variants_CP133676.csv, and rename CHROM to HDRNA_01_K01 in the header, summarize chr and plasmids SNPs of a sample together to a single list, save as an Excel-file.

9, code of summarize_snippy_res.py

  1.     import pandas as pd
  2.     import glob
  3.     import argparse
  4.     import os
  6.     #python3 summarize_snps_indels.py snippy_HDRNA_01/snippy
  8.     #The following record for 2365295 is wrong, since I am sure in the HDRNA_01_K010, it should be a 'G', since in HDRNA_01_K010.csv
  9.     #CP133676,2365295,snp,A,G,G:178 A:0
  10.     #
  11.     #The current output is as follows:
  12.     #CP133676,2365295,A,snp,A,A,A,A,A,A,A,A,A,A,None,,,,,,None,None
  13.     #CP133676,2365295,A,snp,A,A,A,A,A,A,A,A,A,A,nan,,,,,,nan,nan
  14.     #grep -v "None,,,,,,None,None" summary_snps_indels.csv > summary_snps_indels_.csv
  15.     #BUG: CP133676,2365295,A,snp,A,A,A,A,A,A,A,A,A,A,nan,,,,,,nan,nan
  17.     import pandas as pd
  18.     import glob
  19.     import argparse
  20.     import os
  22.     def main(base_directory):
  23.         # List of isolate identifiers
  24.         isolates = [f"HDRNA_01_K0{i}" for i in range(1, 11)]
  25.         expected_columns = ["CHROM", "POS", "REF", "ALT", "TYPE", "EFFECT", "LOCUS_TAG", "GENE", "PRODUCT"]
  27.         # Find all CSV files in the directory and its subdirectories
  28.         csv_files = glob.glob(os.path.join(base_directory, '**', '*.csv'), recursive=True)
  30.         # Create an empty DataFrame to store the summary
  31.         summary_df = pd.DataFrame()
  33.         # Iterate over each CSV file
  34.         for file_path in csv_files:
  35.             # Extract isolate identifier from the file name
  36.             isolate = os.path.basename(file_path).replace('.csv', '')
  37.             df = pd.read_csv(file_path)
  39.             # Ensure all expected columns are present, adding missing ones as empty columns
  40.             for col in expected_columns:
  41.                 if col not in df.columns:
  42.                     df[col] = None
  44.             # Extract relevant columns
  45.             df = df[expected_columns]
  47.             # Ensure consistent data types
  48.             df = df.astype({"CHROM": str, "POS": int, "REF": str, "ALT": str, "TYPE": str, "EFFECT": str, "LOCUS_TAG": str, "GENE": str, "PRODUCT": str})
  50.             # Add the isolate column with the ALT allele
  51.             df[isolate] = df["ALT"]
  53.             # Drop columns that are not needed for the summary
  54.             df = df.drop(["ALT"], axis=1)
  56.             if summary_df.empty:
  57.                 summary_df = df
  58.             else:
  59.                 summary_df = pd.merge(summary_df, df, on=["CHROM", "POS", "REF", "TYPE", "EFFECT", "LOCUS_TAG", "GENE", "PRODUCT"], how="outer")
  61.         # Fill missing values with the REF allele for each isolate column
  62.         for isolate in isolates:
  63.             if isolate in summary_df.columns:
  64.                 summary_df[isolate] = summary_df[isolate].fillna(summary_df["REF"])
  65.             else:
  66.                 summary_df[isolate] = summary_df["REF"]
  68.         # Rename columns to match the required format
  69.         summary_df = summary_df.rename(columns={
  70.             "CHROM": "CHROM",
  71.             "POS": "POS",
  72.             "REF": "REF",
  73.             "TYPE": "TYPE",
  74.             "EFFECT": "Effect",
  75.             "LOCUS_TAG": "Gene_name",
  76.             "GENE": "Biotype",
  77.             "PRODUCT": "Product"
  78.         })
  80.         # Replace any remaining None or NaN values in the non-isolate columns with empty strings
  81.         summary_df = summary_df.fillna("")
  83.         # Add empty columns for Impact, Functional_Class, Codon_change, Protein_and_nucleotide_change, Amino_Acid_Length
  84.         summary_df["Impact"] = ""
  85.         summary_df["Functional_Class"] = ""
  86.         summary_df["Codon_change"] = ""
  87.         summary_df["Protein_and_nucleotide_change"] = ""
  88.         summary_df["Amino_Acid_Length"] = ""
  90.         # Reorder columns
  91.         cols = ["CHROM", "POS", "REF", "TYPE"] + isolates + ["Effect", "Impact", "Functional_Class", "Codon_change", "Protein_and_nucleotide_change", "Amino_Acid_Length", "Gene_name", "Biotype"]
  92.         summary_df = summary_df.reindex(columns=cols)
  94.         # Remove duplicate rows
  95.         summary_df = summary_df.drop_duplicates()
  97.         # Save the summary DataFrame to a CSV file
  98.         output_file = os.path.join(base_directory, "summary_snps_indels.csv")
  99.         summary_df.to_csv(output_file, index=False)
  101.         print("Summary CSV file created successfully at:", output_file)
  103.     if __name__ == "__main__":
  104.         parser = argparse.ArgumentParser(description="Summarize SNPs and Indels from CSV files.")
  105.         parser.add_argument("directory", type=str, help="Base directory containing the CSV files in subdirectories.")
  106.         args = parser.parse_args()
  108.         main(args.directory)

10, code of merge_snps_indels.py

  1.     import pandas as pd
  2.     import argparse
  3.     import os
  5.     def merge_files(input_file1, input_file2, output_file):
  6.         # Read the input files
  7.         df1 = pd.read_csv(input_file1)
  8.         df2 = pd.read_csv(input_file2, sep='\t')
  10.         # Merge the dataframes on the 'POS' column, keeping only the rows that have common 'POS' values
  11.         merged_df = pd.merge(df2, df1[['POS']], on='POS', how='inner')
  13.         # Remove duplicate rows
  14.         merged_df.drop_duplicates(inplace=True)
  16.         # Save the merged dataframe to the output file
  17.         merged_df.to_csv(output_file, index=False)
  19.         print("Merged file created successfully at:", output_file)
  21.     if __name__ == "__main__":
  22.         parser = argparse.ArgumentParser(description="Merge two SNP and Indel files based on the 'POS' column.")
  23.         parser.add_argument("input_file1", type=str, help="Path to the first input file (summary_snps_indels.csv).")
  24.         parser.add_argument("input_file2", type=str, help="Path to the second input file (All_SNPs_indels_annotated.txt).")
  25.         parser.add_argument("output_file", type=str, help="Path to the output file.")
  26.         args = parser.parse_args()
  28.         merge_files(args.input_file1, args.input_file2, args.output_file)

like unlike

点赞本文的读者

还没有人对此文章表态

本文有评论

没有评论

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

最受欢迎文章

  1. Motif Discovery in Biological Sequences: A Comparison of MEME and HOMER
  2. Why Do Significant Gene Lists Change After Adding Additional Conditions in Differential Gene Expression Analysis?
  3. Calling peaks using findPeaks of HOMER
  4. Updating Human Gene Identifiers using Ensembl BioMart: A Step-by-Step Guide
  5. pheatmap vs heatmap.2
  6. Should the inputs for GSVA be normalized or raw?
  7. Setup conda environments
  8. PiCRUST2 Pipeline for Functional Prediction and Pathway Analysis in Metagenomics
  9. Kraken2 Installation and Usage Guide
  10. File format for single channel analysis of Agilent microarray data with Limma?

最新文章

最多评论文章

推荐相似文章

Variant calling (inter-host + intra-host) for Data_Pietschmann_229ECoronavirus_Mutations_2025 (via docker own_viral_ngs) v2

Structural Variant Calling for Nanopore Sequencing (edited)

Scaffolding and finishing an assembly with a reference genome

Snakefile

© 2023 XGenes.com Impressum