RNAseq-踩坑02 -- Aligment 比对率低

RNAseq的第二步就是Aligment, 一般都是有参考基因组的比对，据同事介绍，比对这一步，比对率再80%以上才算是正常的，很多样本的比对结果都能达到90%以上。

然而，然而，然而，对于第一次跑流程的我，第一个项目，居然不到50% ！！！

QC 统计， 均在90%以上，说明测序数据质量尚可

QC

遇到问题：Aligment 统计，比对率 在 28%-47%之间。

Aligment

解决思路：找到比对不上的序列，进行blast 比对，查看比对上的是什么。

具体排查：

1. Aligment 下会有 比对上的fq ID; (A)
   QC 中有 所有fq ID （B）
   QC 中不含Aligment_fq_ID 的剩下所有ID 为 比对不上的ID （B-A）

2. 根据相应软件(blast)，查看这些比对不上的ID 到底是什么原因导致的。

简单介绍用到的软件

1. samtools

如今序列比对已成为各种生物学分析中不可缺少的重要环节，通过将未知的基因片段与已知具体信息的基因或基因组进行比较，并分析其中的相同部分与差异部分，就可以得到该基因片段SNP位点、所属物种以及可能具有的生物学功能等重要信息。sam与bam是两种最常用的比对结果输出文件格式，（如转录组STAR分析软件输出的比对结果为.bam文件等）
bam文件格式是sam文件的二进制格式，占用的存储空间更小，更利于节省存储资源，而且bam文件的计算处理也更快，但二进制无法直接查看，这就需要借助于工具查看了。

• Samtools 软件就是用于处理sam与bam格式的工具软件，能够实现二进制查看、格式转换、排序及合并等功能，结合sam格式中的flag、tag等信息，还可以完成比对结果的统计汇总。同时利用linux中的grep、awk等操作命令，还可以大大扩展samtools的使用范围与功能

基本命令
samtools view [options] <输入bam文件>

• 主要参数:
  -b：以bam格式输出
  -u：以未压缩的bam格式输出，一般与linux命令配合使用
  -h：在sam输出中包含header信息
  -H：只输出header信息
  -f [INT]：只输出在比对flag中包含该整数的序列信息
  -F [INT]：跳过比对flag中含有该整数的序列
  -o [file]：标准输出结果文件

bam/sam文件 每一列的内容

• 1. read名称，通常包括测序平台等信息
• 2. SAM标记（Flag），没有mapping的标记为“ * ”
• 3. chromosome
• 4. 比对上的位置，注意是从1开始计数。
• 5. MAPQ（mapping quality，描述比对的质量，数字越大，特异性越高，说明该read比对到参考基因组上的位置越唯一）
• 6. CIGAR字串，记录插入，删除，错配以及splice junctions（后剪切拼接的接头）
• 7. mate名称，记录mate pair信息
• 8. mate的位置
• 9. 模板的长度
• 10. read序列
• 11. read质量
• 12. 程序用标记

fastq文件格式
FASTQ文件中每个序列通常有四行：

• 序列标识以及相关的描述信息，以‘@’开头；
• 第二行是序列
• 第三行以‘+’开头，后面是序列标示符、描述信息，或者什么也不加
• 第四行，是质量信息，和第二行的序列相对应，每一个序列都有一个质量评分，根据评分体系的不同，每个字符的含义表示的数字也不相同。

fasta文件格式
在生物信息学中，FASTA格式是一种基于文本的、用于表示核苷酸序列或氨基酸序列的格式。
FASTA文件以 序列表示 和 序列 作为一个基本单元，各行记录信息如下：

• 第一行是由大于号">"开头的任意文字说明。用于序列标记，为了保证后续分析软件能够区分每条序列，单个序列的标识必须具有唯一性；

• 第二行开始为序列本身。只允许使用既定的核苷酸或氨基酸编码符号。通常核苷酸符号大小写均可，而氨基酸常用大写字母。使用时应注意有些程序对大小写有明确要求。文件每行的字母一般不应超过80个字符。

  
  fasta

• 将FASTA序列与质量数据放到一起，就变成了FASTQ

• 去掉FASTQ文件的后两行=FASTA

2. blast

• blastn: 核酸序列比对
• blastp: 蛋白质序列比对
• blastx: 核酸翻译成蛋白质序列，再与蛋白质序列库比对
• tblastn:核酸序列库翻译成蛋白质序列库，再与蛋白质序列比对
• tblastx:核酸与核酸序列库都翻译为蛋白质序列，再在蛋白质水平上比对
• makeblastdb:建立自定义的比对序列库

使用方法
安装路径\bin\blastn -query 比对序列文件名.fasta -db database -out 比对结果文件名 -evalue le-5 -outfmt 0(或其他比对结果参数)

• 可以看到，blastn 要求的是fasta格式的数据。

实现具体步骤 (eg. pLVX-ShRNA-1)

1. 提取Aligment之后（bam文件）匹配上的fq ID
   先读取Align.sorted.bam 文件,提取bam第一列（f1）fq ID列，并保存到Paired_ID_list.txt文件里

• samtool view + bam文件
  $ samtools view Alignment/pLVX-ShRNA-1/Align.sorted.bam | cut -f1 >pLVX-ShRNA-1_unpaired_blast/Paired_ID_list.txt
  Paired_ID_list

>>> import pandas as pd
>>> import numpy as np
>>> id_list = pd.read_table("Paired_ID_list.txt",sep = "\t",header=None)
>>> id_list.head()
                                         0
0   A00202:784:HKJKKDSXY:3:1138:3369:12305
1  A00202:784:HKJKKDSXY:3:1123:15953:30968
2   A00202:784:HKJKKDSXY:3:2241:6894:18818
3   A00202:784:HKJKKDSXY:3:1138:3369:12305
4   A00202:784:HKJKKDSXY:3:2241:6894:18818
>>> unilist = np.unique(id_list.values)
>>> print("matchID unique number = ", len(unilist))
matchID unique number =  11704064

2. 提取QC之后的所有的fq ID

   
   trim.out
   before QC
   after QC

$ less ./QC/pLVX-ShRNA-1/pLVX-ShRNA-1_val_1.fq.gz
$ less ./QC/pLVX-ShRNA-1/pLVX-ShRNA-1_val_2.fq.gz

为了测试，这里只提取前1000条fastq数据，

$ zcat QC/pLVX-ShRNA-1/pLVX-ShRNA-1_val_1.fq.gz | head -n 1000 > pLVX-ShRNA-1_unpaired_blast/pLVX-ShRNA-1_val_1_tmp
$ zcat QC/pLVX-ShRNA-1/pLVX-ShRNA-1_val_2.fq.gz | head -n 1000 > pLVX-ShRNA-1_unpaired_blast/pLVX-ShRNA-1_val_2_tmp

检查一下行数

$ wc -l pLVX-ShRNA-1_val_1_tmp 
1000 pLVX-ShRNA-1_val_1_tmp

下面就是进行ID匹配了。

image.png

pLVX-ShRNA-1_val_1_tmp中仍旧是fastq的标准4行的格式，我们需要提取其中的seqID，然后检查该seqID是否在上面的Paired_ID_list.txt当中，我们需要提取不在的部分fastq序列进行分析。

• 提取fastq的seqID
  方法1：根据fastq标准格式的规律，提取seqID

fastqidlst = []
with open("pLVX-ShRNA-1_val_1_tmp",'r') as fastq:
    for line in fastq:
        if line.startswith('@'):
            fastqid = line.strip().split()[0][1:]
            fastqidlst.append(fastqid)

方法2：借助于 python Bio::SeqIO

from Bio import SeqIO
fastqidlst = []
for record in SeqIO.parse("pLVX-ShRNA-1_val_1_tmp",format="fastq"):
        fastqidlst.append(record.id)

但是，我们需要的是unpaired seqID对应的fastq序列。

seqIO可以一步到位，筛选seqID, 并保存seqID对应的fastq文件

## 根据 id 来获取（过滤）fastq
fq1 = open("pLVX-ShRNA-1_val_1_tmp", "r")
fq2 = open("pLVX-ShRNA-1_val_2_tmp", "r")
OUT1=open("pLVX-ShRNA-1_val_1_tmp_unpaired1","w")
OUT2=open("pLVX-ShRNA-1_val_2_tmp_unpaired2","w")

for record in SeqIO.parse(fq1,format="fastq"):
    if record.id not in unilist:
        OUT1.write(record.format('fastq'))

for record in SeqIO.parse(fq2,format="fastq"):
    if record.id not in unilist:
        OUT2.write(record.format('fastq'))

查看fastq文件
less -S pLVX-ShRNA-1_val_2_tmp_unpaired2
pLVX-ShRNA-1_val_2_tmp有 1000条数据，居然有496条没匹配上。

wc -l pLVX-ShRNA-1_val_2_tmp_unpaired2 
496 pLVX-ShRNA-1_val_2_tmp_unpaired2

3. 利用blastn比对unpaired fastq
   blastn输入的是fasta格式，首先需要将fastq转成fasta
   awk '{if(NR%4 == 1){print ">" substr($0, 2)}}{if(NR%4 == 2){print}}' pLVX-ShRNA-1_val_2_tmp_unpaired2 > pLVX-ShRNA-1_val_2_tmp_unpaired2.fasta
   检查一下，看是不是fasta格式的文件

$ less -S pLVX-ShRNA-1_val_2_tmp_unpaired2.fasta
>A00202:784:HKJKKDSXY:3:1101:5412:1000 2:N:0:GACTAGTA+TCTTTCCC
TGTCGCTCCATCAAGCTTTCGCTCATTGTGGAAAATTCCTTACTGCTGCCTCCCGTAGGAGTCTGGGCCGTATCTCAGTCCCAGTGTGGCCGTACAGCCTCTCGGCCCGGCTAAACATCATCGTCTTGGTAGGCCATTACCCTAC

blastn进行核酸序列比对

/home/Software/ncbi-blast-2.11.0+/bin/blastn -task blastn -query pLVX-ShRNA-1_val_2_tmp_unpaired2.fasta \
-db /home/database/nt_blast/nt -show_gis -evalue 1e-3 -max_target_seqs 5 -num_threads 6 -out unpaired2_blast.nt.aln \
-outfmt "6 qseqid qlen sseqid slen qstart qend sstart send length pident nident bitscore evalue staxids stitle"

bastnresult

发现都是支原体污染，问了一下养过细胞的小伙伴，说是细胞被支原体污染是件很常见的事情啊啊啊啊啊啊啊啊啊啊啊啊啊