search

Hollenbach-lab / PING

An R-based bioinformatic pipeline to determine killer-cell immunoglobulin-like receptor (KIR) copy number and high-resolution genotypes from short-read sequencing data.
MIT License
8 stars 6 forks source link

PING for WES and WGS data #8

Open xuxiaochen1209 opened 2 years ago

xuxiaochen1209 commented 2 years ago

Thanks a lot for your hardworking on PING. I was trying to use my WES and WGS data to do KIR genotyping, but it returned with partial results. manualCopyNumberFrame.csv and finalAlleleCalls.csv are not found in the output directory. Both my WES and WGS data are over 100X (yes, the WGS data is also 100X and I think the coverage should not be a problem). And I found your paper published in 2016 also used WES or WGS data to validate this pipeline. So WES or WGS should be OK, right?

I attach part of my rstudio console for your reference:

/usr/bin/bowtie2 -x Resources/extractor_resources/reference/output -5 3 -3 7 -L 20 -i S,1,0.5 --score-min L,0,-0.187 -I 75 -X 1000 --no-unal -p 16 -1 /home/dell/XX/Sequanta/IP02WGS/C22002690LD02-AO-F-P15-2_R1.fastq.gz -2 /home/dell/XX/Sequanta/IP02WGS/C22002690LD02-AO-F-P15-2_R2.fastq.gz -S /home/dell/XX/PING/KIR_output_IP02WGS/extractedFastq/C22002690LD02-AO-F-P15-2_R.sam --al-conc-gz /home/dell/XX/PING/KIR_output_IP02WGS/extractedFastq/C22002690LD02-AO-F-P15-2_RKIR%.fastq.gz Use of uninitialized value $bt2_args[18] in join or string at /usr/bin/bowtie2 line 423. Use of uninitialized value $bt2_args[19] in join or string at /usr/bin/bowtie2 line 423. Use of uninitialized value in exists at /usr/bin/bowtie2 line 81. Use of uninitialized value in exists at /usr/bin/bowtie2 line 81. Use of uninitialized value $bt2_args[18] in join or string at /usr/bin/bowtie2 line 459. Use of uninitialized value $bt2_args[19] in join or string at /usr/bin/bowtie2 line 459. 1058339564 reads; of these: 1058339564 (100.00%) were paired; of these: 1057271393 (99.90%) aligned concordantly 0 times 16585 (0.00%) aligned concordantly exactly 1 time 1051586 (0.10%) aligned concordantly >1 times

1057271393 pairs aligned concordantly 0 times; of these:
  71 (0.00%) aligned discordantly 1 time
----
1057271322 pairs aligned 0 times concordantly or discordantly; of these:
  2114542644 mates make up the pairs; of these:
    2110232716 (99.80%) aligned 0 times
    102006 (0.00%) aligned exactly 1 time
    4207922 (0.20%) aligned >1 times

0.30% overall alignment rate

Successfully extracted KIR reads for C22002690LD02-AO-F-P15-2_R Cleaning up alignment files.

----- PING2_extractor is complete. Extracted reads are deposited in /home/dell/XX/PING/KIR_output_IP02WGS/extractedFastq -----

PING2 gene content and copy number --------------------------------------

source('Resources/genotype_alignment_functions.R')

Writing allele SNP resource files. KIR3DP1 KIR2DS5 KIR2DL3 KIR2DP1 KIR2DS3 KIR2DS2 KIR2DL4 KIR3DL3 KIR3DL1 KIR3DS1 KIR2DL2 KIR3DL2 KIR2DS4 KIR2DL1 KIR2DS1 KIR2DL5 Finished. Allele SNPs written to KIR_output_IP02WGS//alleleFiles as [LOCUS]_alleleSNPs.csv> source('Resources/alleleSetup_functions.R')

cat('\n\n----- Moving to PING gene content and copy determination -----')

----- Moving to PING gene content and copy determination ----->

sampleList <- ping_copy.graph(sampleList=sampleList,threads=threads,resultsDirectory=outDir$path,forceRun=F,onlyKFF=F,fullAlign = F) # set forceRun=T if you want to force alignments Current working directory: /home/dell/XX/PING Max read count set at: 50000

Reading in the KFF probelist file: /home/dell/XX/PING/Resources/gc_resources/probelist_2021_01_24.csv

Processing C22002690LD02-AO-F-P15-2_R

Counting KFF primer matches. Normalizing KFF primer matches. Determining KIR locus presence/absence

Finished with presence/absence determination, moving to copy number determination.

Current used memory: 1253724832

Performing bowtie2 alignment for this sample.

/usr/bin/bowtie2 -x /home/dell/XX/PING/Resources/gc_resources/filled_kir_reference/KIR_compact_filled -5 0 -3 6 -N 0 --end-to-end -p 16 --score-min "L,-2,-0.08" -I 75 -X 1000 -1 /home/dell/XX/PING/KIR_output_IP02WGS/extractedFastq/C22002690LD02-AO-F-P15-2_R_KIR_1.fastq.gz -2 /home/dell/XX/PING/KIR_output_IP02WGS/extractedFastq/C22002690LD02-AO-F-P15-2_R_KIR_2.fastq.gz --no-unal -a --np 1 --mp 2,2 --rdg 1,1 --rfg 1,1 -S /home/dell/XX/PING/KIR_output_IP02WGS/gc_bam_files/C22002690LD02-AO-F-P15-2_R.sam Use of uninitialized value $bt2_args[26] in join or string at /usr/bin/bowtie2 line 423. Use of uninitialized value $bt2_args[27] in join or string at /usr/bin/bowtie2 line 423. Use of uninitialized value in exists at /usr/bin/bowtie2 line 81. Use of uninitialized value in exists at /usr/bin/bowtie2 line 81. Use of uninitialized value $bt2_args[26] in join or string at /usr/bin/bowtie2 line 459. Use of uninitialized value $bt2_args[27] in join or string at /usr/bin/bowtie2 line 459. 1068171 reads; of these: 1068171 (100.00%) were paired; of these: 91112 (8.53%) aligned concordantly 0 times 31906 (2.99%) aligned concordantly exactly 1 time 945153 (88.48%) aligned concordantly >1 times

91112 pairs aligned concordantly 0 times; of these:
  1 (0.00%) aligned discordantly 1 time
----
91111 pairs aligned 0 times concordantly or discordantly; of these:
  182222 mates make up the pairs; of these:
    131620 (72.23%) aligned 0 times
    2792 (1.53%) aligned exactly 1 time
    47810 (26.24%) aligned >1 times

93.84% overall alignment rate

Successfully aligned C22002690LD02-AO-F-P15-2_R to /home/dell/XX/PING/Resources/gc_resources/filled_kir_reference/KIR_compact_filled

/usr/bin/samtools view -@16 /home/dell/XX/PING/KIR_output_IP02WGS/gc_bam_files/C22002690LD02-AO-F-P15-2_R.sam -o /home/dell/XX/PING/KIR_output_IP02WGS/gc_bam_files/C22002690LD02-AO-F-P15-2_R.bam

Successfully converted /home/dell/XX/PING/KIR_output_IP02WGS/gc_bam_files/C22002690LD02-AO-F-P15-2_R.sam to /home/dell/XX/PING/KIR_output_IP02WGS/gc_bam_files/C22002690LD02-AO-F-P15-2_R.bam

Reading in /home/dell/XX/PING/KIR_output_IP02WGS/gc_bam_files/C22002690LD02-AO-F-P15-2_R.sam Setting alignment scores. Setting locus names. Setting read lengths. Removing read alignments that do not include alignment scores. Running cross-mapped read separation with readBoost.thresh (alignment score buffer) = 2 Rescued 22399 reads. Formatting reads.

Setting read start positions.

Setting read end positions. Sorting read table by locus and alignment coordinates.

Generating SNP tables for C22002690LD02-AO-F-P15-2_R: KIR2DL1 KIR2DL2 KIR2DL3 KIR2DL4 KIR2DL5 KIR2DP1 KIR2DS2 KIR2DS3-- fewer than 10 aligned reads, skipping -- KIR2DS4 KIR3DL1 KIR3DL2 KIR3DL3 KIR3DP1 KIR3DS1

----- Finished with alignment! -----

Moving on to copy number graphing. Skipping 0 samples that had fewer than 500 KIR3DL3 reads. Generating copy number graphs... [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone

ALL FINISHED!!> sampleList <- ping_copy.manual_threshold(sampleList=sampleList,resultsDirectory=outDir$path,use.threshFile = T) # this function sets copy thresholds Current working directory: /home/dell/XX/PING Running manual thresholding on: /home/dell/XX/PING/KIR_output_IP02WGS

Found kffPresenceFrame.csv in /home/dell/XX/PING/KIR_output_IP02WGS. Loading these results. Found locusCountFrame.csv in /home/dell/XX/PING/KIR_output_IP02WGS. Loading these results. Found locusRatioFrame.csv in /home/dell/XX/PING/KIR_output_IP02WGS. Loading these results. Found manualCopyThresholds.csv in /home/dell/XX/PING/KIR_output_IP02WGS/manualCopyThresholds.csv. Loading these results. Generating copy number graphs... [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone [WARNING] Deprecated: --self-contained. use --embed-resources --standalone

wesleymarin commented 2 years ago

Hello! So the PING pipeline has not been tested or validated for WES data. There are a couple roadblocks to implementing WES analysis because PING relies on intronic sequence to differentiate between similar alleles.

However, WGS is a different story, and we just made adjustments to PING to be able to handle WGS data, specifically 1000 Genomes Project data, which can be found in the wgs_snakemake branch. (the publication detailing these changes is pending)

However, since your data is much higher depth then 1000GP, I would recommend staying with the normal workflow as you have done, but there are two important considerations to take into account.

First, we have noticed that our extraction step does not work very well for WGS data. Specifically, there seems to be parts of the genome that have high sequence similarity to an intronic region of KIR3DL2, causing a huge pileup of of misaligned reads that get brought into the workflow. The effective strategy that we used with the 1000GP data was to extract from the aligned cram files (aligned to HG38), generate fastq files from the extracted reads, then send those reads through the PING workflow. The BED file used for extraction can be found here https://github.com/Hollenbach-lab/PING/blob/wgs_snakemake/Resources/1000Genomes_resources/kir_regions.bed. I don't know if you have HG38 aligned cram of bam files, but doing the extraction this way generates much cleaner data.

Second, it is really important to separate out the data by ethnicity for copy number determination. Because the allelic makup of KIR differs so much between ethnicities, having a mixed group can jumble up the copy number results and lead to incorrect genotype calls.

With those points out of the way, looking at the output you provided it looks like you are only sending a single sample through. For the copy number determination step, it needs multiple samples to differentiate between copy number groupings. If this is the only sample you have, then id suggest first doing the BAM/CRAM extraction based on the bed file i linked previously, then share the copy number graphs with me and i can help determine the copy numbers.