I have modified R1 reads ID by adding length of poly-T stretches in R2 reads. After that I have aligned R1 reads to genome. Now I have sam files. I can use HT-seq to count number of reads aligned to features. In addition to counts of reads aligned to feature, I would get average T length. Here is example of a reads from sam files.
J00113:322:0001:30:0:GT 16 IV 1450285 22 150M * 0 0 ACCAGTAGTGTGTCTTCTCTTTGCCTTGGCAGCCCAGTTGTGAGATCTAGTCTTAGCGGATGGGTAACCACAAGAGGAACAGGTCTTCTTTTGAACATGGAAAGAACGACGACCAC
ATCTGTTACACAAGGTGTGAGATTTGATTCTCCG 7-<-FJFJJFJJFFJJ<JJFJAF7JFFAAFJFF<FAJAJFFA-7JJ7F<AAFJJ<JFFFJJFJJFJF7<7-F77<-7JF7A-AJFJFJJJJJF<FJFJJFA7FF<--7-<--F-F77J<JFJF-<<<FFJJJJFF<JFJJJJFJA-AA-A AS:i:-23 X
30 in sequence ID field is the length of polyT stretch in R2 files which did not use for alignment. Thanks in advance for helping me with idea about extracting average T length and counts associated with features.
•
written
by
biplab
1 answer
Have htseq-count (or featureCounts, which is much faster) output a BAM/SAM file with the alignments labeled as to which gene (if any) they were assigned (the -o option). You then parse that to get the auxiliary tag it appends. That holds the gene to which the read was assigned and you can appropriately increment a value in a hash after parsing the read ID.
Log in to answer this question.