RNA-seq之比对（我还没做完）

昨天的质控跑完了，今天继续，虽然我还没做完。没做完。没做完。。。
fastqc过后，生成报告中，发现有adaptor的存在，所以需要去掉adaptor

1. 去除adaptor

1.1 去除adaptor测试

#将SRR1039508_1.fastq.gz  ./SRR1039508_2.fastq.gz的adaptor 去掉，放在  ../clean的路径下，输入文件为./SRR1039508_1.fastq.gz  ./SRR1039508_2.fastq.gz
trim_galore --phred33 -q 25 -e 0.1 --length 36 --stringency 3 --paired -o ../clean  ./SRR1039508_1.fastq.gz  ./SRR1039508_2.fastq.gz

1.2 批量去除fastq.gz 中的adaptor

ls /home/vip11/project/rna/raw_data/*_1.fastq.gz >1
ls /home/vip11/project/rna/raw_data/*_2.fastq.gz >2
paste 1 2 >config
cat >qc_clean.sh #写脚本
dir='/home/vip11/project/rna/clean'
cat $1|while read id
do
 arr=(${id})
 fq1=${arr[0]}
 fq2=${arr[1]}
trim_galore --phred33 -q 25 -e 0.1 --length 36 --stringency 3 --paired -o $dir $fq1 $fq2
done
:wq #保存退出
nohup bash qc_clean.sh config & #后台运行脚本

2. 建立索引

下载建立好的索引或是自己构建，此处用已经下载好的索引

3. 批量比对

3.1 批量比对前进行测试

#批量取出每个id的前1000个
mkdir project/rna/test && cd project/rna/test 
ls $project/rna/clean/*gz |while read id;do (zcat ${id}|head -1000>  $(basename ${id} ".gz"));done
#用id=id=SRR1039508的文件进行测试
hisat2 -p 2 -x /home/qmcui/database/reference/index/hisat/hg38/genome -1 /home/vip11/project/rna/clean/SRR1039508_1_val_1.fq  -2 /home/vip11/project/rna/clean/SRR1039508_2_val_2.fq -S SRR1039508.hisat.sam
#批量测试
hisat2 -p 2 -x /home/qmcui/database/reference/index/hisat/hg38/genome -1 $/home/vip11/project/rna/test/{id}_1_val_1.fq  -2 $/home/vip11/project/rna/test/{id}_2_val_2.fq -S $/home/vip11/project/rna/alignment{id}.hisat.sam

3.2 批量比对的代码

cd project/rna/clean

cat >alignment.sh
ls *gz|cut -d"_" -f 1 |sort -u |while read id;do
ls -lh ${id}_1_val_1.fq.gz   ${id}_2_val_2.fq.gz 
hisat2 -p 2 -x /home/qmcui/database/reference/index/hisat/hg38/genome -1 ${id}_1_val_1.fq.gz   -2 ${id}_2_val_2.fq.gz  -S ${id}.hisat.sam 1>/home/vip11/project/rna/alignment/${id}.log 2>&1
done

nohup bash alignment.sh &

4. sam 转bam

ls *.sam|while read id ;do (samtools sort -O bam -@ 2 -o $(basename ${id} ".sam").bam   ${id});done
rm *.sam 
#为bam建立索引
ls *.bam |xargs -i samtools index {}
#reads的比对情况
ls *.bam |xargs -i samtools flagstat -@ 2  {}  >
ls *.bam |while read id ;do ( samtools flagstat -@ 1 $id >  $(basename ${id} ".bam").flagstat  );done
source deactivate

写在最后：

1. sra 转为 fastq.gz 后，需要用fastqc生成报告，在当前文件夹菜面进行
   ls *gz |xargs fastqc -t 10 #调用10个线程同时处理
   ls *.sam|head #head 前面输出的前10个文件
   ls *.sam |xgars head #head 前面输出的文件的前10个内容
2. 查看占用多少内存
   du -h
   构建索引：不同软件的构建方法
   https://www.jianshu.com/p/071c1757ded1

批量杀任务的其中一种方法
 vip11是用户名；python是提交任务中的搜索关键词；用awk输出第二列的PID;循环用kill执行杀掉任务的命令
ps -ef|grep vip11|grep python|awk '{print $2}'|while read id ;do kill $id;done