only unique mapped paired-end reads
1
0
Entering edit mode
15 months ago
Ric ▴ 350

Hi, I would like to only unique mapped paired-end reads. Reads which are not anymore in pair I would like to remove. Are the below steps correctly to do it?

  1. aligned pair-end data with bowtie2.
  2. Removed duplicates with MarkDuplicates
  3. kept aligned reads only with mapping score above 40 with samtools
  4. fixed read information with
  5. filtered only proper paired with samtools .

usage: sh unique_reads_dedup_pbs.sh /path/to/FQdata

for r1 in $(find $1 -name "*.bam"); do

#cat <<eof qsub="" &lt;<eof<="" p="">

!/bin/bash -l

PBS -N $r1

PBS -l walltime=43:00:00

PBS -j oe

PBS -l mem=190G

PBS -l ncpus=2

PBS -m bea

cd \$PBS_O_WORKDIR

conda activate picard picard -Xmx490g MarkDuplicates ASSUME_SORT_ORDER=coordinate REMOVE_DUPLICATES=true INPUT=$r1 OUTPUT=${r1}.dedup.bam METRICS_FILE=${r1}.dedup.metrics.txt conda deactivate

conda activate hisat2 samtools index ${r1}.dedup.bam conda deactivate

conda activate hisat2 samtools view -Sb -q 40 ${r1}.dedup.bam > ${r1}.dedup.q40.bam samtools index ${r1}.dedup.q40.bam conda deactivate

conda activate picard TMP_DIR=pwd/tmp mkdir \$TMP_DIR echo \$TMP_DIR picard -Xmx190g FixMateInformation \ I=${r1} \ O=${r1}.dedup.q40.FixMate.bam \ ASSUME_SORTED=false \ SORT_ORDER=coordinate conda deactivate

conda activate hisat2 samtools view -@ 12 -hb -f 0x2 -F 2316 ${r1}.dedup.q40.FixMate.bam > ${r1}.dedup.q40.FixMate.clean.bam samtools index ${r1}.dedup.q40.FixMate.clean.bam conda deactivate

EOF

done

Thank you in advance

alignment next-gen sequence • 310 views
ADD COMMENT
0
Entering edit mode
15 months ago
gb ★ 1.9k

Think you need to add -F 2 to samtools view. Here you can find which options you have https://broadinstitute.github.io/picard/explain-flags.html

ADD COMMENT
0
Entering edit mode

In the first or second samtools view command?

ADD REPLY
0
Entering edit mode

Think second is better, you are already doing -F 2316. You could change this to 2317 or 2319

ADD REPLY
0
Entering edit mode

Thank you, 2317 looks good and now I can also remove -f 0x2

ADD REPLY

Login before adding your answer.

Traffic: 2330 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6