This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Retrieve reads overlapping indels

I am doing allelic-specific RNA-seq analysis. Given a list of variants, is there an efficient program to separate reads at each variant position based on whether the variant is present? (especially for variant=indels since they accounted for 10% of my variant list). This looks quite obvious when visualized by IGV but I did not find a program doing it.

(Note: preferably I will need the output to be reads instead of other summarized statistics

Thanks

variants reads indel
• 1,108 views • link
updated by Medhat Medhat • written by yliueagle yliueagle

1 answer

Extract intervals as a bed file (do not forget to add padding like 30 bp).
Later, you need to use something like that:

 function extract_fastq {
   myFile=$1
   while read -r line; do a=$(echo $line | awk '{print $1":"$2"-"$3}'); samtools view  "${myFile}".bam  "$a" | awk '{print $1"\t"$10}' |  sort -u -k1,1 |  awk '{print ">"$1"\n"$2}' | bgzip > "${a/:/_}".fq.gz; done < "${myFile}".bed
 }

To use it:

extract_fastq your_bed_file
• 0 views • link
written by Medhat Medhat

Log in to answer this question.

More posts like this