H-pylori: Generate BAM for 1000+ samples by using 4 different reference

[lxwgcool]

Three different purposes

• Generate BAM for 1000+ samples by using 4 different reference
• Generate Mapped (unique) BAM
• Generate Unmapped BAM

[lxwgcool]

Working Directory

• Data Root

  /scratch/lix33/lxwg/Project/H_pylori/Processed

• BAM

  
  /scratch/lix33/lxwg/Project/H_pylori/Processed/BAM

Contains 3 different types of BAM, including Raw, Mapped (Unique) and Unmapped

Raw BAM /scratch/lix33/lxwg/Project/H_pylori/Processed/BAM/*/Raw

Mapped BAM (Unique) /scratch/lix33/lxwg/Project/H_pylori/Processed/BAM/*/Mapped

UnMapped BAM /scratch/lix33/lxwg/Project/H_pylori/Processed/BAM/*/UnMapped

* Variant Calling results

/scratch/lix33/lxwg/Project/H_pylori/Processed/CSV

* Log

/scratch/lix33/lxwg/Project/H_pylori/Processed/Log (1) The calling command line and related history are recorded here

[lxwgcool]

The method of getting unique mapped BAM

• CMD

  samtools view -@ ${threads} -q 10 -F 4 -F 256 -h ${input.bam} | grep -v -E -e '\bXA:Z:' -e '\bSA:Z:' | samtools view -b - > ${output.bam}

• Details

  • -q 10: only consider the MAPQ >= 10
  • -F 4: do not consider unmapped reads
  • -F 256: do not consider not primary alignment
  • -h: keep header
  • '\bXA:Z:': do not consider alternative hits
  • '\bSA:Z:': do not consider other canonical alignments in a chimeric alignment

• Reference

  • MAPQ (1) https://genome.sph.umich.edu/wiki/Mapping_Quality_Scores (2) https://www.researchgate.net/figure/number-of-alignments-per-sample-High-quality-MapQ-10-alignments-per-sample-are-shown_fig3_312960209
  • Flags (1) https://broadinstitute.github.io/picard/explain-flags.html (2) https://www.biostars.org/p/150063/ (3) https://www.samformat.info/sam-format-alignment-tags
  • Other people's strategy (1) https://gist.github.com/sbamin/27acf13f2a28161efbf89a273559bca4

[lxwgcool]

The method of getting unmapped reads

• Command line

  samtools view -h -b -f 4 ${input.bam} > ${output.bam}

• Details

  • -f 4: only consider read unmapped

[lxwgcool]

Output Explain

• 3 different types of BAM
  • Raw: original BAM
  • Mapped: unique mapped BAM
  • Unmapped: unmapped BAM

• In the Mapped(unique)/Unmapped output directory, each sample has three outputs, including
  • sorted BAM
  • BAM index
  • gene.list: contains all genes (reads name) contained by current BAM.

[lxwgcool]

New Features and Fixed Bugs

New Features

• Append index into Reads Name
• Re-generate an updated ffa.gz file
• Use the updated reads ffa file for mapping

Fixed bug

• Solve the issue caused by "$" in command line to get the correct gene list.

[lxwgcool]

New Function: find disordered Gene

1: Easy case
(1) Normal, discored (1), normal

2: Midean Level case
(1) front goes to back (304)
    * still pick it out)

3: Hard level case
(1) Number of reads from back to front
    * not report all
    * only find if there is discord in these number of reads
(2) order goes to reverse order 
    * just ignore it
(3) number of reads together but the order arrangment is pretty complex.
(4) the case of Chimeric

[lxwgcool]

Big Progress: new function of calculating Methylation Position

1. Use The raw position in gene fna file and the position in gffv3 file to collect all methylations that belong to current gene
2. Align Gene back to reference
3. Find the real mapping position of the methylation.

Hey Difei,

Everything is all set.

The results have been uploaded into: 
/scratch/lix33/lxwg/Project/H_pylori/Processed/Methylation

Please check the files below:

The related github issue is: 
https://github.com/NCI-CGR/IlluminaSequencingAnalysis/issues/57

The latest commit is: 
https://github.com/NCI-CGR/IlluminaSequencingAnalysis/commit/f49dc3cf0531e8477567029643b34bcd2a9d0d08

The column “LiftoverAlignedPosRef(Methylation)” in CSV file is what you are looking for.

I am happy this job could be done in 2022.

Have a great new year
Best
Xin

[lxwgcool]

Reference

1: Explain CIGAR Mask https://samtools.sourceforge.net/samtools/bam/PDefines/PDefines.html

2: How to do offset for different CIGAR code https://sourceforge.net/p/samtools/mailman/message/29373646/

[lxwgcool]

New function: Add three more info into the CSV file

• m6A/m4C from the GFF3
• The base in the reference at the same position like the snp in VCF file
• The third is the strand info from GFF3

[lxwgcool]

Modify the logic: Add mask to orphan Methylation

• Orphan Methylation

  • Cannot find the related base in reference

• Set the position to its nearest reference

• Add mask "(*)"

  • example

• The related reference base will be set as "*"