环境是以前的ATAC
mkdir pfcchow
cd pfcchow #新文件夹
cat SRR_Acc_List.txt |while read id; do (nohup prefetch $id &); done
转换格式
mv */* ./ #把子文件夹里的所有文件移动到当前
ls ./*sra |while read id; do (nohup fastq-dump --split-3 --gzip --skip-technical --clip ${id} &); done
QC
mkdir ./QCclean
fastqc -t 8 -o ./QCclean *.fastq.gz
批量QC
永远在报错，暂时没管，理论上应该
multiqc ./
#看着像是python2.x不兼容，后来重新创了个py3.11的环境
conda create --name py3.11 python=3.11
conda activate py3.11
conda install multiqc
报错 ModuleNotFoundError: No module named 'typing_extensions'
pip install typing-extensions解决了
trimmomatic修剪
**双端
for filename in SRR25056478 SRR25056479 SRR25056480 SRR25056481 SRR25056482 SRR25056483 #真没学过linux不知道这句应该怎么改，直接提取有重复，只会写单端
do
nohup trimmomatic PE -threads 10 ./${filename}_1.fastq.gz ./${filename}_2.fastq.gz ./${filename}_1.clean.fq.gz ./${filename}_1.drop.fq.gz ./${filename}_2.clean.fq.gz ./${filename}_2.drop.fq.gz HEADCROP:15 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:16 &
done
**单端
for file in *.fastq.gz;
do
nohuptrimmomatic SE -threads 10 ./"$file" ./"${file%.fastq.gz}.clean.fastq.gz"HEADCROP:15 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:24 >"${file%.fastq.gz}.log" 2>&1 &
done
参数：
PE: 双端模式。需要两个输入文件，正向测序序列和反向测序序列：sample_R1.fastq.gz和sample_R2.fastq.gz；
以及四个输出文件：sample_1.clean.fq.gz和sample_1.drop.fq.gz；sample_2.clean.fq.gz和sample_2.drop.fq.gz
-threads: 线程数
HEADCROP: 从 reads 的开头切掉指定数量的碱基，即从fastqc中看的
LEADING: 从 reads 的开头切除质量值低于阈值的碱基
TRAILING: 从 reads 的末尾开始切除质量值低于阈值的碱基
SLIDINGWINDOW: 从 reads 的5' 端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗
MINLEN: 如果经过剪切后 reads 的长度低于阈值则丢弃这条reads，即上面的计算方法：（n-m）×40%（40%为什么我没查）
想写成脚本但是报错了不知道为什么，用上面的跑了，脚本👇
#!/bin/bash
for file in /home/data/t210424/yes/data/pfcchow/*.fastq.gz; do
file_name=$(basename "$file" .fastq.gz)
filename=$(echo $file_name | cut -d'_' -f1)
nohup trimmomatic PE -threads 10 /home/data/t210424/yes/data/pfcchow/${filename}_1.fastq.gz /home/data/t210424/yes/data/pfcchow/${filename}_2.fastq.gz /home/data/t210424/yes/data/pfcchow/trim/${filename}_1.clean.fq.gz /home/data/t210424/yes/data/pfcchow/trim/${filename}_1.drop.fq.gz /home/data/t210424/yes/data/pfcchow/trim/${filename}_2.clean.fq.gz /home/data/t210424/yes/data/pfcchow/trim/${filename}_2.drop.fq.gz HEADCROP:15 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:16 > /home/data/t210424/yes/data/pfcchow/trim/"${filename}.log" &
done
之前用trimglore但是QC以后更差了，原因不知道，代码记录一下
mkdir ./clean
cd clean
ls ~/yes/data/pfcchow/*_1.fastq.gz >1
ls ~/yes/data/pfcchow/*_2.fastq.gz >2
paste 1 2 > config
cat >qc.sh
bin_trim_galore=trim_galore #这一步本是多此一举，但是很多时候并不想将软件写入环境变量，就可是在这里使用绝对路径，之后引用变量就可以了。
dir='/data1/lixianlong/YST/clean' #将结果导向的位置。
cat $1 |while read id
do
arr=(${id}) # $id表示$1的每一行，这里是config的每一行，以空格键分割。
fq1=${arr[0]} #空格键分割的第一个元素
fq2=${arr[1]} #空格键分割的第二个元素，曾健明说直接拿来用就好，为什么这样他也不知道。余顺太发现这个其实属于《LINUX命令行与shell脚本》6.7 数组变量的内容。
nohup $bin_trim_galore -q 25 --phred33--length 36 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
$ bash qc.shconfig #config就是qc.sh脚本的$1
再次QC并调整
mkdir -p ./trim && mv *.fq.gz ./trim/
mkdir ./QCclean
fastqc -t 8 -o ./QCclean *.clean.fq.gz
hisat2索引+比对
#下载 H. sapiens, UCSC hg38 genome
wget https://cloud.biohpc.swmed.edu/index.php/s/hg38/download
#下载 H. sapiens, UCSC hg38 genome_tran
wget https://cloud.biohpc.swmed.edu/index.php/s/hg38_tran/download
#下载 M. musculus, UCSC mm10 genome
wget https://cloud.biohpc.swmed.edu/index.php/s/mm10/download
mapping.sh
#!/bin/bash
dir=/home/data/t210424/yes/data/pfcchow/trim/bwa_sam/
hisat2_mm10_index='/home/data/refdir/database/server/reference/index/hisat/mm10/genome'
for id in SRR25056478 SRR25056479 SRR25056480 SRR25056481 SRR25056482 SRR25056483
do
fq1=${id}_1.clean.fq.gz
fq2=${id}_2.clean.fq.gz
nohup hisat2 -p 10 -x $hisat2_mm10_index -1 $fq1 -2 $fq2 -S $dir/${id}.sam &
done
报错原因好像是length不能短于20，重新trim了一次，QC倒是差别不大但是能mapping了
bash mapping.sh
sam格式转bam格式
ls *.sam|while read id; do (samtools sort -O bam -@ 5 -o $(basename $id ".sam").bam $id);done
#sam文件转成bam文件，文件变的更小，sort是很耗时间的。-@表示线程数量
$ ls *bam|xargs -i samtools index {}
#为bam文件建立索引。{}类似于$id，就是ls *bam的每一个结果会通过xargs传给后面的参数，这个参数就是{}。
原文摘过来的，这个没做，我的igv在Windows下看的
#排好序并构建好索引，之后就可以使用IGV看了。
$ samtools view SRR1039508.bam |less -S # bam文件的查看使用samtools view，因为bam文件是按染色体sort过的，
$ samtools tview 08.bam #曾健明说这个用不着，平时都用IGV，因为IGV是可视化的，所以跳转很方便。
$ / #会出现goto，随便输入一个位置chr1:65891003看下，通过这种方式可快速到达指定位置。
$ samtools tview --reference /data1/lixianlong/YST/reference/genome/hg19/hg19.fa 08.bam #可以看到比对的情况，输入 chr1:220142425再试下。
$ samtools flagstat 08.bam # 也可以对bam文件进行统计，得到比对率之类的信息。samtools flagstat不能对sam进行统计，只能对bam统计。
counts
gtf="/home/data/refdir/database/annotation/gencode/v34/gencode.v34.annotation.gtf.gz"
gtf="/home/data/refdir/database/server/reference/gtf/gencode/gencode.vM23.annotation.gtf.gz"
featureCounts -T 5 -p -t exon -g gene_id -a $gtf -o all.id.txt *.bam 1>counts.id.log 2>&1 &
# -T线程，-p双端测序，-t exon根据外显子来count，-g gene_id 表示count之后取基因名，得到的 all.id.txt文件就是表达矩阵，这个featureCounts有很多参数可调。$ source deactivate
-a 指定注释文件
-o 指定结果输出目录及文件名
-p 能用在paired-end的情况中，会统计fragment而不统计read
-t 指定feature的类型，默认是exon，当然gtf里面还有gene、CDS或者直接以feature命名的分类方式。
#其它参数：
-f参数该参数设置后统计的是feature层面（默认是exon）的参数，如果不设置则是直接统计meta-feature参数（即一个gene中的多个exon）
