昨天的质控跑完了，今天继续，虽然我还没做完。没做完。没做完。。。
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
写在最后：
- sra 转为 fastq.gz 后，需要用fastqc生成报告，在当前文件夹菜面进行
ls *gz |xargs fastqc -t 10#调用10个线程同时处理
ls *.sam|head#head 前面输出的前10个文件
ls *.sam |xgars head#head 前面输出的文件的前10个内容 - 查看占用多少内存
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
