I'm trying to call SNPs from RNA-seq data.
The workflow I'm trying to follow is: tophat2 --> remove duplicates --> local realignment and base quality recalibration by GATK --> call SNPs by samtools mpileup or GATK
These are some specific questions:
1. Should I filter the alignment from tophat2 to get only the uniquely mapped reads (reads that could only be mapped to single location on the genome)? If yes, how should I filter them?
I know there is a --multihits (-g) parameter in tophat that can be set to 1 to report only one hit per mapped read. But this is not the strict definition of "unique mapping". I saw people saying that in the alignment sam file, if MAPQ=50, that means this read can only be mapped to a single position. Could I, for example, run tophat2 with g=20 (default), and then select only the aligned reads with MAPQ=50?
2. GATK cannot accept N in cigar. How should I deal with the aligned reads with N in the cigar? GATK has two options: the first is to remove these reads; the second is to allow N in cigar. The former strategy may lose a lot of information, while the latter is at risk (I don't know what the risk is exactly).
3. To use samtools mpileup, how should I set the -C (parameter for adjusting mapQ) parameter? I know it recommends 50 for bwa aligned reads, but what it is for tophat2 alignment?
4. If you have experience in this kind of data analysis, can you see any other potential problems in the pipeline or the parameters of each program? Do you have any recommendation?
By the way, I know there are a lot of caveats of using RNA-seq data to call SNPs. The major purpose of this project is to measure expression, but calling SNPs is to try to get more use of the data (or the money).
Thanks a lot for any suggestions or advice!
The workflow I'm trying to follow is: tophat2 --> remove duplicates --> local realignment and base quality recalibration by GATK --> call SNPs by samtools mpileup or GATK
These are some specific questions:
1. Should I filter the alignment from tophat2 to get only the uniquely mapped reads (reads that could only be mapped to single location on the genome)? If yes, how should I filter them?
I know there is a --multihits (-g) parameter in tophat that can be set to 1 to report only one hit per mapped read. But this is not the strict definition of "unique mapping". I saw people saying that in the alignment sam file, if MAPQ=50, that means this read can only be mapped to a single position. Could I, for example, run tophat2 with g=20 (default), and then select only the aligned reads with MAPQ=50?
2. GATK cannot accept N in cigar. How should I deal with the aligned reads with N in the cigar? GATK has two options: the first is to remove these reads; the second is to allow N in cigar. The former strategy may lose a lot of information, while the latter is at risk (I don't know what the risk is exactly).
3. To use samtools mpileup, how should I set the -C (parameter for adjusting mapQ) parameter? I know it recommends 50 for bwa aligned reads, but what it is for tophat2 alignment?
4. If you have experience in this kind of data analysis, can you see any other potential problems in the pipeline or the parameters of each program? Do you have any recommendation?
By the way, I know there are a lot of caveats of using RNA-seq data to call SNPs. The major purpose of this project is to measure expression, but calling SNPs is to try to get more use of the data (or the money).
Thanks a lot for any suggestions or advice!